diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h
index 3b2fe81..e5298af 100644
--- a/include/GooseFEM/MeshQuad4.h
+++ b/include/GooseFEM/MeshQuad4.h
@@ -1,319 +1,325 @@
/**
Generate mesh with 4-noded quadrilateral elements.
\file MeshQuad4.h
\copyright Copyright 2017. Tom de Geus. All rights reserved.
\license This project is released under the GNU Public License (GPLv3).
*/
#ifndef GOOSEFEM_MESHQUAD4_H
#define GOOSEFEM_MESHQUAD4_H
#include "config.h"
namespace GooseFEM {
namespace Mesh {
namespace Quad4 {
// pre-allocation
namespace Map {
class FineLayer2Regular;
}
// Regular mesh: equi-sized elements
class Regular {
public:
Regular() = default;
Regular(size_t nelx, size_t nely, double h = 1.0);
// size
size_t nelem() const; // number of elements
size_t nnode() const; // number of nodes
size_t nne() const; // number of nodes-per-element
size_t ndim() const; // number of dimensions
size_t nelx() const; // number of elements in x-direction
size_t nely() const; // number of elements in y-direction
double h() const; // edge size
/**
Get the element type.
\return GooseFEM::Mesh::ElementType().
*/
ElementType getElementType() const;
// mesh
xt::xtensor<double, 2> coor() const; // nodal positions [nnode, ndim]
xt::xtensor<size_t, 2> conn() const; // connectivity [nelem, nne]
// boundary nodes: edges
xt::xtensor<size_t, 1> nodesBottomEdge() const;
xt::xtensor<size_t, 1> nodesTopEdge() const;
xt::xtensor<size_t, 1> nodesLeftEdge() const;
xt::xtensor<size_t, 1> nodesRightEdge() const;
// boundary nodes: edges, without corners
xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
xt::xtensor<size_t, 1> nodesRightOpenEdge() const;
// boundary nodes: corners (including aliases)
size_t nodesBottomLeftCorner() const;
size_t nodesBottomRightCorner() const;
size_t nodesTopLeftCorner() const;
size_t nodesTopRightCorner() const;
size_t nodesLeftBottomCorner() const;
size_t nodesLeftTopCorner() const;
size_t nodesRightBottomCorner() const;
size_t nodesRightTopCorner() const;
// DOF-numbers for each component of each node (sequential)
xt::xtensor<size_t, 2> dofs() const;
// DOF-numbers for the case that the periodicity if fully eliminated
xt::xtensor<size_t, 2> dofsPeriodic() const;
// periodic node pairs [:,2]: (independent, dependent)
xt::xtensor<size_t, 2> nodesPeriodic() const;
// front-bottom-left node, used as reference for periodicity
size_t nodesOrigin() const;
// element numbers as matrix
xt::xtensor<size_t, 2> elementgrid() const;
private:
double m_h; // elementary element edge-size (in all directions)
size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h")
size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h")
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
static const size_t m_nne = 4; // number of nodes-per-element
static const size_t m_ndim = 2; // number of dimensions
};
// Mesh with fine middle layer, and coarser elements towards the top and bottom
class FineLayer {
public:
FineLayer() = default;
FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1);
// Reconstruct class for given coordinates / connectivity
FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
// size
size_t nelem() const; // number of elements
size_t nnode() const; // number of nodes
size_t nne() const; // number of nodes-per-element
size_t ndim() const; // number of dimensions
size_t nelx() const; // number of elements in x-direction
size_t nely() const; // number of elements in y-direction
double h() const; // edge size
// edge size, per row of elements (in units of "h")
xt::xtensor<size_t, 1> elemrow_nhx() const;
xt::xtensor<size_t, 1> elemrow_nhy() const;
xt::xtensor<size_t, 1> elemrow_nelem() const;
// type
ElementType getElementType() const;
// mesh
xt::xtensor<double, 2> coor() const; // nodal positions [nnode, ndim]
xt::xtensor<size_t, 2> conn() const; // connectivity [nelem, nne]
// elements in the middle (fine) layer
xt::xtensor<size_t, 1> elementsMiddleLayer() const;
// extract elements along a layer
xt::xtensor<size_t, 1> elementsLayer(size_t layer) const;
// select region of elements from 'matrix' of element numbers
xt::xtensor<size_t, 1> elementgrid_ravel(
std::vector<size_t> rows_start_stop,
std::vector<size_t> cols_start_stop) const;
/**
Select region of elements from 'matrix' of element numbers around an element:
square box with edge-size ``(2 * size + 1) * h``, around ``element``.
\param element The element around which to select elements.
\param size Edge size of the square box encapsulating the selected element.
\param periodic Assume the mesh periodic.
\returns List of elements.
*/
xt::xtensor<size_t, 1> elementgrid_around_ravel(
size_t element,
size_t size,
bool periodic = true);
/**
Select region of elements from 'matrix' of element numbers around an element:
left/right from ``element`` (on the same layer).
\param element The element around which to select elements.
\param left Number of elements to select to the left.
\param right Number of elements to select to the right.
\param periodic Assume the mesh periodic.
\returns List of elements.
*/
// -
xt::xtensor<size_t, 1> elementgrid_leftright(
size_t element,
size_t left,
size_t right,
bool periodic = true);
// boundary nodes: edges
xt::xtensor<size_t, 1> nodesBottomEdge() const;
xt::xtensor<size_t, 1> nodesTopEdge() const;
xt::xtensor<size_t, 1> nodesLeftEdge() const;
xt::xtensor<size_t, 1> nodesRightEdge() const;
// boundary nodes: edges, without corners
xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
xt::xtensor<size_t, 1> nodesRightOpenEdge() const;
// boundary nodes: corners (including aliases)
size_t nodesBottomLeftCorner() const;
size_t nodesBottomRightCorner() const;
size_t nodesTopLeftCorner() const;
size_t nodesTopRightCorner() const;
size_t nodesLeftBottomCorner() const;
size_t nodesLeftTopCorner() const;
size_t nodesRightBottomCorner() const;
size_t nodesRightTopCorner() const;
// DOF-numbers for each component of each node (sequential)
xt::xtensor<size_t, 2> dofs() const;
// DOF-numbers for the case that the periodicity if fully eliminated
xt::xtensor<size_t, 2> dofsPeriodic() const;
// periodic node pairs [:,2]: (independent, dependent)
xt::xtensor<size_t, 2> nodesPeriodic() const;
// front-bottom-left node, used as reference for periodicity
size_t nodesOrigin() const;
// mapping to 'roll' periodically in the x-direction,
// returns element mapping, such that: new_elemvar = elemvar[elem_map]
xt::xtensor<size_t, 1> roll(size_t n);
private:
double m_h; // elementary element edge-size (in all directions)
double m_Lx; // mesh size in "x"
size_t m_nelem; // number of elements
size_t m_nnode; // number of nodes
static const size_t m_nne = 4; // number of nodes-per-element
static const size_t m_ndim = 2; // number of dimensions
xt::xtensor<size_t, 1> m_nelx; // number of elements in "x" (*)
xt::xtensor<size_t, 1> m_nnd; // total number of nodes in the main node layer (**)
xt::xtensor<size_t, 1> m_nhx; // element size in x-direction (*)
xt::xtensor<size_t, 1> m_nhy; // element size in y-direction (*)
xt::xtensor<int, 1> m_refine; // refine direction (-1:no refine, 0:"x" (*)
xt::xtensor<size_t, 1> m_startElem; // start element (*)
xt::xtensor<size_t, 1> m_startNode; // start node (**)
// (*) per element layer in "y"
// (**) per node layer in "y"
void init(size_t nelx, size_t nely, double h, size_t nfine = 1);
void map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular;
};
// Mesh mappings
namespace Map {
// Return "FineLayer"-class responsible for generating a connectivity
// Throws if conversion is not possible
GooseFEM::Mesh::Quad4::FineLayer FineLayer(
const xt::xtensor<double, 2>& coor,
const xt::xtensor<size_t, 2>& conn);
// Refine a regular mesh: sub-divide elements in several smaller elements
class RefineRegular {
public:
// Constructors
RefineRegular() = default;
RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny);
// return the coarse or the fine mesh objects
GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const;
GooseFEM::Mesh::Quad4::Regular getFineMesh() const;
// elements of the fine mesh per element of the coarse mesh
xt::xtensor<size_t, 2> getMap() const;
// map field
xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 1>& data) const; // scalar per el
xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
xt::xtensor<double, 4> mapToCoarse(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
// map field
xt::xtensor<double, 1> mapToFine(const xt::xtensor<double, 1>& data) const; // scalar per el
xt::xtensor<double, 2> mapToFine(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
xt::xtensor<double, 4> mapToFine(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
private:
// the meshes
GooseFEM::Mesh::Quad4::Regular m_coarse;
GooseFEM::Mesh::Quad4::Regular m_fine;
// mapping
xt::xtensor<size_t, 1> m_fine2coarse;
xt::xtensor<size_t, 1> m_fine2coarse_index;
xt::xtensor<size_t, 2> m_coarse2fine;
};
class FineLayer2Regular {
public:
// constructor
FineLayer2Regular() = default;
FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh);
// return either of the meshes
GooseFEM::Mesh::Quad4::Regular getRegularMesh() const;
GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const;
// elements of the Regular mesh per element of the FineLayer mesh
// and the fraction by which the overlap is
std::vector<std::vector<size_t>> getMap() const;
std::vector<std::vector<double>> getMapFraction() const;
- // map field
- xt::xtensor<double, 1> mapToRegular(const xt::xtensor<double, 1>& data) const; // scalar per el
- xt::xtensor<double, 2> mapToRegular(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
- xt::xtensor<double, 4> mapToRegular(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
+ /**
+ Map integration point quantities to Regular.
+
+ \tparam T The type of the data (e.g. ``double``).
+ \tparam rank Rank of the data.
+ \param arg The data.
+ \return The mapped data.
+ */
+ template <class T, size_t rank>
+ xt::xtensor<T, rank> mapToRegular(const xt::xtensor<T, rank>& arg) const;
private:
// the "FineLayer" mesh to map
GooseFEM::Mesh::Quad4::FineLayer m_finelayer;
// the new "Regular" mesh to which to map
GooseFEM::Mesh::Quad4::Regular m_regular;
// mapping
std::vector<std::vector<size_t>> m_elem_regular;
std::vector<std::vector<double>> m_frac_regular;
};
} // namespace Map
} // namespace Quad4
} // namespace Mesh
} // namespace GooseFEM
#include "MeshQuad4.hpp"
#endif
diff --git a/include/GooseFEM/MeshQuad4.hpp b/include/GooseFEM/MeshQuad4.hpp
index 71ac6cc..ac7fc7a 100644
--- a/include/GooseFEM/MeshQuad4.hpp
+++ b/include/GooseFEM/MeshQuad4.hpp
@@ -1,1646 +1,1618 @@
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_MESHQUAD4_HPP
#define GOOSEFEM_MESHQUAD4_HPP
#include "MeshQuad4.h"
namespace GooseFEM {
namespace Mesh {
namespace Quad4 {
inline Regular::Regular(size_t nelx, size_t nely, double h) : m_h(h), m_nelx(nelx), m_nely(nely)
{
GOOSEFEM_ASSERT(m_nelx >= 1ul);
GOOSEFEM_ASSERT(m_nely >= 1ul);
m_nnode = (m_nelx + 1) * (m_nely + 1);
m_nelem = m_nelx * m_nely;
}
inline size_t Regular::nelem() const
{
return m_nelem;
}
inline size_t Regular::nnode() const
{
return m_nnode;
}
inline size_t Regular::nne() const
{
return m_nne;
}
inline size_t Regular::ndim() const
{
return m_ndim;
}
inline size_t Regular::nelx() const
{
return m_nelx;
}
inline size_t Regular::nely() const
{
return m_nely;
}
inline double Regular::h() const
{
return m_h;
}
inline ElementType Regular::getElementType() const
{
return ElementType::Quad4;
}
inline xt::xtensor<double, 2> Regular::coor() const
{
xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
xt::xtensor<double, 1> x =
xt::linspace<double>(0.0, m_h * static_cast<double>(m_nelx), m_nelx + 1);
xt::xtensor<double, 1> y =
xt::linspace<double>(0.0, m_h * static_cast<double>(m_nely), m_nely + 1);
size_t inode = 0;
for (size_t iy = 0; iy < m_nely + 1; ++iy) {
for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
ret(inode, 0) = x(ix);
ret(inode, 1) = y(iy);
++inode;
}
}
return ret;
}
inline xt::xtensor<size_t, 2> Regular::conn() const
{
xt::xtensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
size_t ielem = 0;
for (size_t iy = 0; iy < m_nely; ++iy) {
for (size_t ix = 0; ix < m_nelx; ++ix) {
ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix);
ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1);
ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix);
ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1);
++ielem;
}
}
return ret;
}
inline xt::xtensor<size_t, 1> Regular::nodesBottomEdge() const
{
return xt::arange<size_t>(m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesTopEdge() const
{
return xt::arange<size_t>(m_nelx + 1) + m_nely * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesLeftEdge() const
{
return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesRightEdge() const
{
return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1) + m_nelx;
}
inline xt::xtensor<size_t, 1> Regular::nodesBottomOpenEdge() const
{
return xt::arange<size_t>(1, m_nelx);
}
inline xt::xtensor<size_t, 1> Regular::nodesTopOpenEdge() const
{
return xt::arange<size_t>(1, m_nelx) + m_nely * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesLeftOpenEdge() const
{
return xt::arange<size_t>(1, m_nely) * (m_nelx + 1);
}
inline xt::xtensor<size_t, 1> Regular::nodesRightOpenEdge() const
{
return xt::arange<size_t>(1, m_nely) * (m_nelx + 1) + m_nelx;
}
inline size_t Regular::nodesBottomLeftCorner() const
{
return 0;
}
inline size_t Regular::nodesBottomRightCorner() const
{
return m_nelx;
}
inline size_t Regular::nodesTopLeftCorner() const
{
return m_nely * (m_nelx + 1);
}
inline size_t Regular::nodesTopRightCorner() const
{
return m_nely * (m_nelx + 1) + m_nelx;
}
inline size_t Regular::nodesLeftBottomCorner() const
{
return nodesBottomLeftCorner();
}
inline size_t Regular::nodesLeftTopCorner() const
{
return nodesTopLeftCorner();
}
inline size_t Regular::nodesRightBottomCorner() const
{
return nodesBottomRightCorner();
}
inline size_t Regular::nodesRightTopCorner() const
{
return nodesTopRightCorner();
}
inline xt::xtensor<size_t, 2> Regular::nodesPeriodic() const
{
xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);
ret(0, 0) = nodesBottomLeftCorner();
ret(0, 1) = nodesBottomRightCorner();
ret(1, 0) = nodesBottomLeftCorner();
ret(1, 1) = nodesTopRightCorner();
ret(2, 0) = nodesBottomLeftCorner();
ret(2, 1) = nodesTopLeftCorner();
size_t i = 3;
xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
i += bot.size();
xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
return ret;
}
inline size_t Regular::nodesOrigin() const
{
return nodesBottomLeftCorner();
}
inline xt::xtensor<size_t, 2> Regular::dofs() const
{
return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
}
inline xt::xtensor<size_t, 2> Regular::dofsPeriodic() const
{
xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
for (size_t j = 0; j < m_ndim; ++j) {
xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
}
return GooseFEM::Mesh::renumber(ret);
}
inline xt::xtensor<size_t, 2> Regular::elementgrid() const
{
return xt::arange<size_t>(m_nelem).reshape({m_nely, m_nelx});
}
inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine)
{
this->init(nelx, nely, h, nfine);
}
inline FineLayer::FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
{
this->map(coor, conn);
}
inline void FineLayer::init(size_t nelx, size_t nely, double h, size_t nfine)
{
GOOSEFEM_ASSERT(nelx >= 1ul);
GOOSEFEM_ASSERT(nely >= 1ul);
m_h = h;
m_Lx = m_h * static_cast<double>(nelx);
// compute element size in y-direction (use symmetry, compute upper half)
// temporary variables
size_t nmin, ntot;
xt::xtensor<size_t, 1> nhx = xt::ones<size_t>({nely});
xt::xtensor<size_t, 1> nhy = xt::ones<size_t>({nely});
xt::xtensor<int, 1> refine = -1 * xt::ones<int>({nely});
// minimum height in y-direction (half of the height because of symmetry)
if (nely % 2 == 0) {
nmin = nely / 2;
}
else {
nmin = (nely + 1) / 2;
}
// minimum number of fine layers in y-direction (minimum 1, middle layer part of this half)
if (nfine % 2 == 0) {
nfine = nfine / 2 + 1;
}
else {
nfine = (nfine + 1) / 2;
}
if (nfine < 1) {
nfine = 1;
}
if (nfine > nmin) {
nfine = nmin;
}
// loop over element layers in y-direction, try to coarsen using these rules:
// (1) element size in y-direction <= distance to origin in y-direction
// (2) element size in x-direction should fit the total number of elements in x-direction
// (3) a certain number of layers have the minimum size "1" (are fine)
for (size_t iy = nfine;;) {
// initialize current size in y-direction
if (iy == nfine) {
ntot = nfine;
}
// check to stop
if (iy >= nely || ntot >= nmin) {
nely = iy;
break;
}
// rules (1,2) satisfied: coarsen in x-direction
if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) {
refine(iy) = 0;
nhy(iy) *= 2;
auto vnhy = xt::view(nhy, xt::range(iy + 1, _));
auto vnhx = xt::view(nhx, xt::range(iy, _));
vnhy *= 3;
vnhx *= 3;
}
// update the number of elements in y-direction
ntot += nhy(iy);
// proceed to next element layer in y-direction
++iy;
// check to stop
if (iy >= nely || ntot >= nmin) {
nely = iy;
break;
}
}
// symmetrize, compute full information
// allocate mesh constructor parameters
m_nhx = xt::empty<size_t>({nely * 2 - 1});
m_nhy = xt::empty<size_t>({nely * 2 - 1});
m_refine = xt::empty<int>({nely * 2 - 1});
m_nelx = xt::empty<size_t>({nely * 2 - 1});
m_nnd = xt::empty<size_t>({nely * 2});
m_startElem = xt::empty<size_t>({nely * 2 - 1});
m_startNode = xt::empty<size_t>({nely * 2});
// fill
// - lower half
for (size_t iy = 0; iy < nely; ++iy) {
m_nhx(iy) = nhx(nely - iy - 1);
m_nhy(iy) = nhy(nely - iy - 1);
m_refine(iy) = refine(nely - iy - 1);
}
// - upper half
for (size_t iy = 0; iy < nely - 1; ++iy) {
m_nhx(iy + nely) = nhx(iy + 1);
m_nhy(iy + nely) = nhy(iy + 1);
m_refine(iy + nely) = refine(iy + 1);
}
// update size
nely = m_nhx.size();
// compute the number of elements per element layer in y-direction
for (size_t iy = 0; iy < nely; ++iy) {
m_nelx(iy) = nelx / m_nhx(iy);
}
// compute the number of nodes per node layer in y-direction
for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
m_nnd(iy) = m_nelx(iy) + 1;
}
for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
m_nnd(iy + 1) = m_nelx(iy) + 1;
}
// compute mesh dimensions
// initialize
m_nnode = 0;
m_nelem = 0;
m_startNode(0) = 0;
// loop over element layers (bottom -> middle, elements become finer)
for (size_t i = 0; i < (nely - 1) / 2; ++i) {
// - store the first element of the layer
m_startElem(i) = m_nelem;
// - add the nodes of this layer
if (m_refine(i) == 0) {
m_nnode += (3 * m_nelx(i) + 1);
}
else {
m_nnode += (m_nelx(i) + 1);
}
// - add the elements of this layer
if (m_refine(i) == 0) {
m_nelem += (4 * m_nelx(i));
}
else {
m_nelem += (m_nelx(i));
}
// - store the starting node of the next layer
m_startNode(i + 1) = m_nnode;
}
// loop over element layers (middle -> top, elements become coarser)
for (size_t i = (nely - 1) / 2; i < nely; ++i) {
// - store the first element of the layer
m_startElem(i) = m_nelem;
// - add the nodes of this layer
if (m_refine(i) == 0) {
m_nnode += (5 * m_nelx(i) + 1);
}
else {
m_nnode += (m_nelx(i) + 1);
}
// - add the elements of this layer
if (m_refine(i) == 0) {
m_nelem += (4 * m_nelx(i));
}
else {
m_nelem += (m_nelx(i));
}
// - store the starting node of the next layer
m_startNode(i + 1) = m_nnode;
}
// - add the top row of nodes
m_nnode += m_nelx(nely - 1) + 1;
}
inline size_t FineLayer::nelem() const
{
return m_nelem;
}
inline size_t FineLayer::nnode() const
{
return m_nnode;
}
inline size_t FineLayer::nne() const
{
return m_nne;
}
inline size_t FineLayer::ndim() const
{
return m_ndim;
}
inline size_t FineLayer::nelx() const
{
return xt::amax(m_nelx)();
}
inline size_t FineLayer::nely() const
{
return xt::sum(m_nhy)();
}
inline double FineLayer::h() const
{
return m_h;
}
inline xt::xtensor<size_t, 1> FineLayer::elemrow_nhx() const
{
return m_nhx;
}
inline xt::xtensor<size_t, 1> FineLayer::elemrow_nhy() const
{
return m_nhy;
}
inline xt::xtensor<size_t, 1> FineLayer::elemrow_nelem() const
{
return m_nelx;
}
inline ElementType FineLayer::getElementType() const
{
return ElementType::Quad4;
}
inline xt::xtensor<double, 2> FineLayer::coor() const
{
// allocate output
xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
// current node, number of element layers
size_t inode = 0;
size_t nely = static_cast<size_t>(m_nhy.size());
// y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
// - allocate
xt::xtensor<double, 1> y = xt::empty<double>({nely + 1});
// - initialize
y(0) = 0.0;
// - compute
for (size_t iy = 1; iy < nely + 1; ++iy) {
y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h;
}
// loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes
for (size_t iy = 0;; ++iy) {
// get positions along the x- and z-axis
xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);
// add nodes of the bottom layer of this element
for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
ret(inode, 0) = x(ix);
ret(inode, 1) = y(iy);
++inode;
}
// stop at middle layer
if (iy == (nely - 1) / 2) {
break;
}
// add extra nodes of the intermediate layer, for refinement in x-direction
if (m_refine(iy) == 0) {
// - get position offset in x- and y-direction
double dx = m_h * static_cast<double>(m_nhx(iy) / 3);
double dy = m_h * static_cast<double>(m_nhy(iy) / 2);
// - add nodes of the intermediate layer
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
for (size_t j = 0; j < 2; ++j) {
ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
ret(inode, 1) = y(iy) + dy;
++inode;
}
}
}
}
// loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes
for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
// get positions along the x- and z-axis
xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);
// add extra nodes of the intermediate layer, for refinement in x-direction
if (m_refine(iy) == 0) {
// - get position offset in x- and y-direction
double dx = m_h * static_cast<double>(m_nhx(iy) / 3);
double dy = m_h * static_cast<double>(m_nhy(iy) / 2);
// - add nodes of the intermediate layer
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
for (size_t j = 0; j < 2; ++j) {
ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
ret(inode, 1) = y(iy) + dy;
++inode;
}
}
}
// add nodes of the top layer of this element
for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
ret(inode, 0) = x(ix);
ret(inode, 1) = y(iy + 1);
++inode;
}
}
return ret;
}
inline xt::xtensor<size_t, 2> FineLayer::conn() const
{
// allocate output
xt::xtensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
// current element, number of element layers, starting nodes of each node layer
size_t ielem = 0;
size_t nely = static_cast<size_t>(m_nhy.size());
size_t bot, mid, top;
// loop over all element layers
for (size_t iy = 0; iy < nely; ++iy) {
// - get: starting nodes of bottom(, middle) and top layer
bot = m_startNode(iy);
mid = m_startNode(iy) + m_nnd(iy);
top = m_startNode(iy + 1);
// - define connectivity: no coarsening/refinement
if (m_refine(iy) == -1) {
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
ret(ielem, 0) = bot + (ix);
ret(ielem, 1) = bot + (ix + 1);
ret(ielem, 2) = top + (ix + 1);
ret(ielem, 3) = top + (ix);
ielem++;
}
}
// - define connectivity: refinement along the x-direction (below the middle layer)
else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
// -- bottom element
ret(ielem, 0) = bot + (ix);
ret(ielem, 1) = bot + (ix + 1);
ret(ielem, 2) = mid + (2 * ix + 1);
ret(ielem, 3) = mid + (2 * ix);
ielem++;
// -- top-right element
ret(ielem, 0) = bot + (ix + 1);
ret(ielem, 1) = top + (3 * ix + 3);
ret(ielem, 2) = top + (3 * ix + 2);
ret(ielem, 3) = mid + (2 * ix + 1);
ielem++;
// -- top-center element
ret(ielem, 0) = mid + (2 * ix);
ret(ielem, 1) = mid + (2 * ix + 1);
ret(ielem, 2) = top + (3 * ix + 2);
ret(ielem, 3) = top + (3 * ix + 1);
ielem++;
// -- top-left element
ret(ielem, 0) = bot + (ix);
ret(ielem, 1) = mid + (2 * ix);
ret(ielem, 2) = top + (3 * ix + 1);
ret(ielem, 3) = top + (3 * ix);
ielem++;
}
}
// - define connectivity: coarsening along the x-direction (above the middle layer)
else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) {
for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
// -- lower-left element
ret(ielem, 0) = bot + (3 * ix);
ret(ielem, 1) = bot + (3 * ix + 1);
ret(ielem, 2) = mid + (2 * ix);
ret(ielem, 3) = top + (ix);
ielem++;
// -- lower-center element
ret(ielem, 0) = bot + (3 * ix + 1);
ret(ielem, 1) = bot + (3 * ix + 2);
ret(ielem, 2) = mid + (2 * ix + 1);
ret(ielem, 3) = mid + (2 * ix);
ielem++;
// -- lower-right element
ret(ielem, 0) = bot + (3 * ix + 2);
ret(ielem, 1) = bot + (3 * ix + 3);
ret(ielem, 2) = top + (ix + 1);
ret(ielem, 3) = mid + (2 * ix + 1);
ielem++;
// -- upper element
ret(ielem, 0) = mid + (2 * ix);
ret(ielem, 1) = mid + (2 * ix + 1);
ret(ielem, 2) = top + (ix + 1);
ret(ielem, 3) = top + (ix);
ielem++;
}
}
}
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::elementsMiddleLayer() const
{
size_t nely = m_nhy.size();
size_t iy = (nely - 1) / 2;
return m_startElem(iy) + xt::arange<size_t>(m_nelx(iy));
}
inline xt::xtensor<size_t, 1> FineLayer::elementsLayer(size_t iy) const
{
GOOSEFEM_ASSERT(iy < m_nelx.size());
size_t n = m_nelx(iy);
if (m_refine(iy) != -1) {
n *= 4;
}
return m_startElem(iy) + xt::arange<size_t>(n);
}
inline xt::xtensor<size_t, 1> FineLayer::elementgrid_ravel(
std::vector<size_t> start_stop_rows,
std::vector<size_t> start_stop_cols) const
{
GOOSEFEM_ASSERT(start_stop_rows.size() == 0 || start_stop_rows.size() == 2);
GOOSEFEM_ASSERT(start_stop_cols.size() == 0 || start_stop_cols.size() == 2);
std::array<size_t, 2> rows;
std::array<size_t, 2> cols;
if (start_stop_rows.size() == 2) {
std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin());
GOOSEFEM_ASSERT(rows[0] <= this->nely());
GOOSEFEM_ASSERT(rows[1] <= this->nely());
}
else {
rows[0] = 0;
rows[1] = this->nely();
}
if (start_stop_cols.size() == 2) {
std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin());
GOOSEFEM_ASSERT(cols[0] <= this->nelx());
GOOSEFEM_ASSERT(cols[1] <= this->nelx());
}
else {
cols[0] = 0;
cols[1] = this->nelx();
}
if (rows[0] == rows[1] || cols[0] == cols[1]) {
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({0});
return ret;
}
// Compute dimensions
auto H = xt::cumsum(m_nhy);
size_t yl = 0;
if (rows[0] > 0) {
yl = xt::argmax(H > rows[0])();
}
size_t yu = xt::argmax(H >= rows[1])();
size_t hx = std::max(m_nhx(yl), m_nhx(yu));
size_t xl = (cols[0] - cols[0] % hx) / hx;
size_t xu = (cols[1] - cols[1] % hx) / hx;
// Allocate output
size_t N = 0;
for (size_t i = yl; i <= yu; ++i) {
// no refinement
size_t n = (xu - xl) * hx / m_nhx(i);
// refinement
if (m_refine(i) != -1) {
n *= 4;
}
N += n;
}
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({N});
// Write output
N = 0;
for (size_t i = yl; i <= yu; ++i) {
// no refinement
size_t n = (xu - xl) * hx / m_nhx(i);
size_t h = hx;
// refinement
if (m_refine(i) != -1) {
n *= 4;
h *= 4;
}
xt::xtensor<size_t, 1> e = m_startElem(i) + xl * h / m_nhx(i) + xt::arange<size_t>(n);
xt::view(ret, xt::range(N, N + n)) = e;
N += n;
}
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::elementgrid_around_ravel(
size_t e,
size_t size,
bool periodic)
{
GOOSEFEM_WIP_ASSERT(periodic == true);
size_t iy = xt::argmin(m_startElem <= e)() - 1;
size_t nel = m_nelx(iy);
GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
auto H = xt::cumsum(m_nhy);
if (2 * size >= H(H.size() - 1)) {
return xt::arange<size_t>(this->nelem());
}
size_t hy = H(iy);
size_t l = xt::argmax(H > (hy - size - 1))();
size_t u = xt::argmax(H >= (hy + size))();
size_t lh = 0;
if (l > 0) {
lh = H(l - 1);
}
size_t uh = H(u);
size_t step = xt::amax(m_nhx)();
size_t relx = (e - m_startElem(iy)) % step;
size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
size_t dx = m_nhx(u);
size_t xl = 0;
size_t xu = nel;
if (mid >= size) {
xl = mid - size;
}
if (mid + size < nel) {
xu = mid + size + 1;
}
xl = xl - xl % dx;
xu = xu - xu % dx;
if (mid - xl < size) {
if (xl < dx) {
xl = 0;
}
else {
xl -= dx;
}
}
if (xu - mid < size) {
if (xu > nel - dx) {
xu = nel;
}
else {
xu += dx;
}
}
auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu});
auto map = this->roll(nroll);
return xt::view(map, xt::keep(ret));
}
inline xt::xtensor<size_t, 1> FineLayer::elementgrid_leftright(
size_t e,
size_t left,
size_t right,
bool periodic)
{
GOOSEFEM_WIP_ASSERT(periodic == true);
size_t iy = xt::argmin(m_startElem <= e)() - 1;
size_t nel = m_nelx(iy);
GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
size_t step = xt::amax(m_nhx)();
size_t relx = (e - m_startElem(iy)) % step;
size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
size_t dx = m_nhx(iy);
size_t xl = 0;
size_t xu = nel;
if (mid >= left) {
xl = mid - left;
}
if (mid + right < nel) {
xu = mid + right + 1;
}
xl = xl - xl % dx;
xu = xu - xu % dx;
if (mid - xl < left) {
if (xl < dx) {
xl = 0;
}
else {
xl -= dx;
}
}
if (xu - mid < right) {
if (xu > nel - dx) {
xu = nel;
}
else {
xu += dx;
}
}
auto H = xt::cumsum(m_nhy);
size_t yl = 0;
if (iy > 0) {
yl = H(iy - 1);
}
auto ret = this->elementgrid_ravel({yl, H(iy)}, {xl, xu});
auto map = this->roll(nroll);
return xt::view(map, xt::keep(ret));
}
inline xt::xtensor<size_t, 1> FineLayer::nodesBottomEdge() const
{
return m_startNode(0) + xt::arange<size_t>(m_nelx(0) + 1);
}
inline xt::xtensor<size_t, 1> FineLayer::nodesTopEdge() const
{
size_t nely = m_nhy.size();
return m_startNode(nely) + xt::arange<size_t>(m_nelx(nely - 1) + 1);
}
inline xt::xtensor<size_t, 1> FineLayer::nodesLeftEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j));
xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1));
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::nodesRightEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j)) =
xt::view(m_startNode, xt::range(i, j)) + xt::view(m_nelx, xt::range(i, j));
xt::view(ret, xt::range(k + 1, l + 1)) =
xt::view(m_startNode, xt::range(k + 1, l + 1)) + xt::view(m_nelx, xt::range(k, l));
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::nodesBottomOpenEdge() const
{
return m_startNode(0) + xt::arange<size_t>(1, m_nelx(0));
}
inline xt::xtensor<size_t, 1> FineLayer::nodesTopOpenEdge() const
{
size_t nely = m_nhy.size();
return m_startNode(nely) + xt::arange<size_t>(1, m_nelx(nely - 1));
}
inline xt::xtensor<size_t, 1> FineLayer::nodesLeftOpenEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j));
xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l));
return ret;
}
inline xt::xtensor<size_t, 1> FineLayer::nodesRightOpenEdge() const
{
size_t nely = m_nhy.size();
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
size_t i = 0;
size_t j = (nely + 1) / 2;
size_t k = (nely - 1) / 2;
size_t l = nely;
xt::view(ret, xt::range(i, j - 1)) =
xt::view(m_startNode, xt::range(i + 1, j)) + xt::view(m_nelx, xt::range(i + 1, j));
xt::view(ret, xt::range(k, l - 1)) =
xt::view(m_startNode, xt::range(k + 1, l)) + xt::view(m_nelx, xt::range(k, l - 1));
return ret;
}
inline size_t FineLayer::nodesBottomLeftCorner() const
{
return m_startNode(0);
}
inline size_t FineLayer::nodesBottomRightCorner() const
{
return m_startNode(0) + m_nelx(0);
}
inline size_t FineLayer::nodesTopLeftCorner() const
{
size_t nely = m_nhy.size();
return m_startNode(nely);
}
inline size_t FineLayer::nodesTopRightCorner() const
{
size_t nely = m_nhy.size();
return m_startNode(nely) + m_nelx(nely - 1);
}
inline size_t FineLayer::nodesLeftBottomCorner() const
{
return nodesBottomLeftCorner();
}
inline size_t FineLayer::nodesRightBottomCorner() const
{
return nodesBottomRightCorner();
}
inline size_t FineLayer::nodesLeftTopCorner() const
{
return nodesTopLeftCorner();
}
inline size_t FineLayer::nodesRightTopCorner() const
{
return nodesTopRightCorner();
}
inline xt::xtensor<size_t, 2> FineLayer::nodesPeriodic() const
{
xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);
ret(0, 0) = nodesBottomLeftCorner();
ret(0, 1) = nodesBottomRightCorner();
ret(1, 0) = nodesBottomLeftCorner();
ret(1, 1) = nodesTopRightCorner();
ret(2, 0) = nodesBottomLeftCorner();
ret(2, 1) = nodesTopLeftCorner();
size_t i = 3;
xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
i += bot.size();
xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
return ret;
}
inline size_t FineLayer::nodesOrigin() const
{
return nodesBottomLeftCorner();
}
inline xt::xtensor<size_t, 2> FineLayer::dofs() const
{
return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
}
inline xt::xtensor<size_t, 2> FineLayer::dofsPeriodic() const
{
xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
for (size_t j = 0; j < m_ndim; ++j) {
xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
}
return GooseFEM::Mesh::renumber(ret);
}
inline xt::xtensor<size_t, 1> FineLayer::roll(size_t n)
{
auto conn = this->conn();
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelem});
// loop over all element layers
for (size_t iy = 0; iy < nely; ++iy) {
// no refinement
size_t shift = n * (m_nelx(iy) / m_nelx(0));
size_t nel = m_nelx(iy);
// refinement
if (m_refine(iy) != -1) {
shift = n * (m_nelx(iy) / m_nelx(0)) * 4;
nel = m_nelx(iy) * 4;
}
// element numbers of the layer, and roll them
auto e = m_startElem(iy) + xt::arange<size_t>(nel);
xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift);
}
return ret;
}
inline void FineLayer::map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
{
GOOSEFEM_ASSERT(coor.shape(1) == 2);
GOOSEFEM_ASSERT(conn.shape(1) == 4);
GOOSEFEM_ASSERT(conn.shape(0) > 0);
GOOSEFEM_ASSERT(coor.shape(0) >= 4);
GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
if (conn.shape(0) == 1) {
this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0));
GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
return;
}
// Identify the middle layer
size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2;
size_t eleft = emid;
size_t eright = emid;
while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) {
eleft--;
}
while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) {
eright++;
}
GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0));
// Get element sizes along the middle layer
auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0);
auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1);
auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 2);
auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0);
auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1);
auto hx = xt::amin(dx)();
auto hy = xt::amin(dy)();
GOOSEFEM_CHECK(xt::allclose(hx, hy));
GOOSEFEM_CHECK(xt::allclose(dx, hx));
GOOSEFEM_CHECK(xt::allclose(dy, hy));
// Extract shape and initialise
size_t nelx = eright - eleft + 1;
size_t nely = static_cast<size_t>((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx);
this->init(nelx, nely, hx);
GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
GOOSEFEM_CHECK(xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange<size_t>(nelx))));
}
namespace Map {
inline RefineRegular::RefineRegular(
const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny)
: m_coarse(mesh)
{
m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h());
xt::xtensor<size_t, 2> elmat_coarse = m_coarse.elementgrid();
xt::xtensor<size_t, 2> elmat_fine = m_fine.elementgrid();
m_coarse2fine = xt::empty<size_t>({m_coarse.nelem(), nx * ny});
for (size_t i = 0; i < elmat_coarse.shape(0); ++i) {
for (size_t j = 0; j < elmat_coarse.shape(1); ++j) {
xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view(
elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx)));
}
}
}
inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getCoarseMesh() const
{
return m_coarse;
}
inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getFineMesh() const
{
return m_fine;
}
inline xt::xtensor<size_t, 2> RefineRegular::getMap() const
{
return m_coarse2fine;
}
inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 1>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
size_t m = m_coarse2fine.shape(0);
size_t n = m_coarse2fine.shape(1);
xt::xtensor<double, 2> ret = xt::empty<double>({m, n});
for (size_t i = 0; i < m; ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, i) = xt::view(data, xt::keep(e));
}
return ret;
}
inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 2>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
size_t m = m_coarse2fine.shape(0);
size_t n = m_coarse2fine.shape(1);
size_t N = data.shape(1);
xt::xtensor<double, 2> ret = xt::empty<double>({m, n * N});
for (size_t i = 0; i < m; ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
for (size_t q = 0; q < data.shape(1); ++q) {
xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
}
}
return ret;
}
inline xt::xtensor<double, 4> RefineRegular::mapToCoarse(const xt::xtensor<double, 4>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
size_t m = m_coarse2fine.shape(0);
size_t n = m_coarse2fine.shape(1);
size_t N = data.shape(1);
xt::xtensor<double, 4> ret = xt::empty<double>({m, n * N, data.shape(2), data.shape(3)});
for (size_t i = 0; i < m; ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
for (size_t q = 0; q < data.shape(1); ++q) {
xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
}
}
return ret;
}
inline xt::xtensor<double, 1> RefineRegular::mapToFine(const xt::xtensor<double, 1>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
xt::xtensor<double, 1> ret = xt::empty<double>({m_coarse2fine.size()});
for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, xt::keep(e)) = data(i);
}
return ret;
}
inline xt::xtensor<double, 2> RefineRegular::mapToFine(const xt::xtensor<double, 2>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
xt::xtensor<double, 2> ret = xt::empty<double>({m_coarse2fine.size(), data.shape(1)});
for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, xt::keep(e)) = xt::view(data, i);
}
return ret;
}
inline xt::xtensor<double, 4> RefineRegular::mapToFine(const xt::xtensor<double, 4>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
xt::xtensor<double, 4> ret =
xt::empty<double>({m_coarse2fine.size(), data.shape(1), data.shape(2), data.shape(3)});
for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
auto e = xt::view(m_coarse2fine, i, xt::all());
xt::view(ret, xt::keep(e)) = xt::view(data, i);
}
return ret;
}
inline FineLayer2Regular::FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh)
: m_finelayer(mesh)
{
// ------------
// Regular-mesh
// ------------
m_regular = GooseFEM::Mesh::Quad4::Regular(
xt::amax(m_finelayer.m_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h);
// -------
// mapping
// -------
// allocate mapping
m_elem_regular.resize(m_finelayer.m_nelem);
m_frac_regular.resize(m_finelayer.m_nelem);
// alias
xt::xtensor<size_t, 1> nhx = m_finelayer.m_nhx;
xt::xtensor<size_t, 1> nhy = m_finelayer.m_nhy;
xt::xtensor<size_t, 1> nelx = m_finelayer.m_nelx;
xt::xtensor<size_t, 1> start = m_finelayer.m_startElem;
// 'matrix' of element numbers of the Regular-mesh
xt::xtensor<size_t, 2> elementgrid = m_regular.elementgrid();
// cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh
xt::xtensor<size_t, 1> cum_nhy =
xt::concatenate(xt::xtuple(xt::zeros<size_t>({1}), xt::cumsum(nhy)));
// number of element layers in y-direction of the FineLayer-mesh
size_t nely = nhy.size();
// loop over layers of the FineLayer-mesh
for (size_t iy = 0; iy < nely; ++iy) {
// element numbers of the Regular-mesh along this layer of the FineLayer-mesh
auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all());
// no coarsening/refinement
// ------------------------
if (m_finelayer.m_refine(iy) == -1) {
// element numbers of the FineLayer-mesh along this layer
xt::xtensor<size_t, 1> el_old = start(iy) + xt::arange<size_t>(nelx(iy));
// loop along this layer of the FineLayer-mesh
for (size_t ix = 0; ix < nelx(iy); ++ix) {
// get the element numbers of the Regular-mesh for this element of the
// FineLayer-mesh
auto block =
xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
// write to mapping
for (auto& i : block) {
m_elem_regular[el_old(ix)].push_back(i);
m_frac_regular[el_old(ix)].push_back(1.0);
}
}
}
// refinement along the x-direction (below the middle layer)
// ---------------------------------------------------------
else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
// element numbers of the FineLayer-mesh along this layer
// rows: coarse block, columns element numbers per block
xt::xtensor<size_t, 2> el_old =
start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
// loop along this layer of the FineLayer-mesh
for (size_t ix = 0; ix < nelx(iy); ++ix) {
// get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
auto block =
xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
// bottom: wide-to-narrow
{
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(block, j, xt::range(j, nhx(iy) - j));
m_elem_regular[el_old(ix, 0)].push_back(e(0));
m_frac_regular[el_old(ix, 0)].push_back(0.5);
for (size_t k = 1; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 0)].push_back(e(k));
m_frac_regular[el_old(ix, 0)].push_back(1.0);
}
m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 0)].push_back(0.5);
}
}
// top: regular small element
{
auto e = xt::view(
block,
xt::range(nhy(iy) / 2, nhy(iy)),
xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 2)].push_back(i);
m_frac_regular[el_old(ix, 2)].push_back(1.0);
}
}
// left
{
// left-bottom: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(block, j, xt::range(0, j + 1));
for (size_t k = 0; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 3)].push_back(e(k));
m_frac_regular[el_old(ix, 3)].push_back(1.0);
}
m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 3)].push_back(0.5);
}
// left-top: regular
{
auto e = xt::view(
block,
xt::range(nhy(iy) / 2, nhy(iy)),
xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 3)].push_back(i);
m_frac_regular[el_old(ix, 3)].push_back(1.0);
}
}
}
// right
{
// right-bottom: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy)));
m_elem_regular[el_old(ix, 1)].push_back(e(0));
m_frac_regular[el_old(ix, 1)].push_back(0.5);
for (size_t k = 1; k < e.size(); ++k) {
m_elem_regular[el_old(ix, 1)].push_back(e(k));
m_frac_regular[el_old(ix, 1)].push_back(1.0);
}
}
// right-top: regular
{
auto e = xt::view(
block,
xt::range(nhy(iy) / 2, nhy(iy)),
xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 1)].push_back(i);
m_frac_regular[el_old(ix, 1)].push_back(1.0);
}
}
}
}
}
// coarsening along the x-direction (above the middle layer)
else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) {
// element numbers of the FineLayer-mesh along this layer
// rows: coarse block, columns element numbers per block
xt::xtensor<size_t, 2> el_old =
start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
// loop along this layer of the FineLayer-mesh
for (size_t ix = 0; ix < nelx(iy); ++ix) {
// get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
auto block =
xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
// top: narrow-to-wide
{
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(
block,
nhy(iy) / 2 + j,
xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1));
m_elem_regular[el_old(ix, 3)].push_back(e(0));
m_frac_regular[el_old(ix, 3)].push_back(0.5);
for (size_t k = 1; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 3)].push_back(e(k));
m_frac_regular[el_old(ix, 3)].push_back(1.0);
}
m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 3)].push_back(0.5);
}
}
// bottom: regular small element
{
auto e = xt::view(
block,
xt::range(0, nhy(iy) / 2),
xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 1)].push_back(i);
m_frac_regular[el_old(ix, 1)].push_back(1.0);
}
}
// left
{
// left-bottom: regular
{
auto e = xt::view(
block,
xt::range(0, nhy(iy) / 2),
xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 0)].push_back(i);
m_frac_regular[el_old(ix, 0)].push_back(1.0);
}
}
// left-top: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e =
xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j));
for (size_t k = 0; k < e.size() - 1; ++k) {
m_elem_regular[el_old(ix, 0)].push_back(e(k));
m_frac_regular[el_old(ix, 0)].push_back(1.0);
}
m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
m_frac_regular[el_old(ix, 0)].push_back(0.5);
}
}
// right
{
// right-bottom: regular
{
auto e = xt::view(
block,
xt::range(0, nhy(iy) / 2),
xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
for (auto& i : e) {
m_elem_regular[el_old(ix, 2)].push_back(i);
m_frac_regular[el_old(ix, 2)].push_back(1.0);
}
}
// right-top: narrow-to-wide
for (size_t j = 0; j < nhy(iy) / 2; ++j) {
auto e = xt::view(
block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy)));
m_elem_regular[el_old(ix, 2)].push_back(e(0));
m_frac_regular[el_old(ix, 2)].push_back(0.5);
for (size_t k = 1; k < e.size(); ++k) {
m_elem_regular[el_old(ix, 2)].push_back(e(k));
m_frac_regular[el_old(ix, 2)].push_back(1.0);
}
}
}
}
}
}
}
inline GooseFEM::Mesh::Quad4::Regular FineLayer2Regular::getRegularMesh() const
{
return m_regular;
}
inline GooseFEM::Mesh::Quad4::FineLayer FineLayer2Regular::getFineLayerMesh() const
{
return m_finelayer;
}
inline std::vector<std::vector<size_t>> FineLayer2Regular::getMap() const
{
return m_elem_regular;
}
inline std::vector<std::vector<double>> FineLayer2Regular::getMapFraction() const
{
return m_frac_regular;
}
-inline xt::xtensor<double, 1>
-FineLayer2Regular::mapToRegular(const xt::xtensor<double, 1>& data) const
+template <class T, size_t rank>
+inline xt::xtensor<T, rank>
+FineLayer2Regular::mapToRegular(const xt::xtensor<T, rank>& data) const
{
GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
- xt::xtensor<double, 1> ret = xt::zeros<double>({m_regular.nelem()});
+ std::array<size_t, rank> shape;
+ std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
+ shape[0] = m_regular.nelem();
- for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
- for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
- ret(m_elem_regular[e][i]) += m_frac_regular[e][i] * data(e);
- }
- }
-
- return ret;
-}
-
-inline xt::xtensor<double, 2>
-FineLayer2Regular::mapToRegular(const xt::xtensor<double, 2>& data) const
-{
- GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
-
- xt::xtensor<double, 2> ret = xt::zeros<double>({m_regular.nelem(), data.shape(1)});
-
- for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
- for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
- xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
- }
- }
-
- return ret;
-}
-
-inline xt::xtensor<double, 4>
-FineLayer2Regular::mapToRegular(const xt::xtensor<double, 4>& data) const
-{
- GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
-
- xt::xtensor<double, 4> ret =
- xt::zeros<double>({m_regular.nelem(), data.shape(1), data.shape(2), data.shape(3)});
+ xt::xtensor<T, rank> ret = xt::zeros<T>(shape);
for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
}
}
return ret;
}
} // namespace Map
} // namespace Quad4
} // namespace Mesh
} // namespace GooseFEM
#endif
diff --git a/include/GooseFEM/VectorPartitioned.hpp b/include/GooseFEM/VectorPartitioned.hpp
index dca41a3..92ba965 100644
--- a/include/GooseFEM/VectorPartitioned.hpp
+++ b/include/GooseFEM/VectorPartitioned.hpp
@@ -1,757 +1,757 @@
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_VECTORPARTITIONED_HPP
#define GOOSEFEM_VECTORPARTITIONED_HPP
#include "Mesh.h"
#include "VectorPartitioned.h"
namespace GooseFEM {
inline VectorPartitioned::VectorPartitioned(
const xt::xtensor<size_t, 2>& conn,
const xt::xtensor<size_t, 2>& dofs,
const xt::xtensor<size_t, 1>& iip)
: m_conn(conn), m_dofs(dofs), m_iip(iip)
{
m_nelem = m_conn.shape(0);
m_nne = m_conn.shape(1);
m_nnode = m_dofs.shape(0);
m_ndim = m_dofs.shape(1);
m_iiu = xt::setdiff1d(m_dofs, m_iip);
m_ndof = xt::amax(m_dofs)() + 1;
m_nnp = m_iip.size();
m_nnu = m_iiu.size();
- m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs);
+ m_part = Mesh::Reorder({m_iiu, m_iip}).apply(m_dofs);
GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)());
GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
}
inline size_t VectorPartitioned::nelem() const
{
return m_nelem;
}
inline size_t VectorPartitioned::nne() const
{
return m_nne;
}
inline size_t VectorPartitioned::nnode() const
{
return m_nnode;
}
inline size_t VectorPartitioned::ndim() const
{
return m_ndim;
}
inline size_t VectorPartitioned::ndof() const
{
return m_ndof;
}
inline size_t VectorPartitioned::nnu() const
{
return m_nnu;
}
inline size_t VectorPartitioned::nnp() const
{
return m_nnp;
}
inline xt::xtensor<size_t, 2> VectorPartitioned::dofs() const
{
return m_dofs;
}
inline xt::xtensor<size_t, 1> VectorPartitioned::iiu() const
{
return m_iiu;
}
inline xt::xtensor<size_t, 1> VectorPartitioned::iip() const
{
return m_iip;
}
inline void VectorPartitioned::copy(
const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
xt::noalias(nodevec_dest) = nodevec_src;
}
inline void VectorPartitioned::copy_u(
const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) < m_nnu) {
nodevec_dest(m, i) = nodevec_src(m, i);
}
}
}
}
inline void VectorPartitioned::copy_p(
const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) >= m_nnu) {
nodevec_dest(m, i) = nodevec_src(m, i);
}
}
}
}
inline void VectorPartitioned::asDofs(
const xt::xtensor<double, 1>& dofval_u,
const xt::xtensor<double, 1>& dofval_p,
xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
dofval.fill(0.0);
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
dofval(m_iiu(d)) = dofval_u(d);
}
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
dofval(m_iip(d)) = dofval_p(d);
}
}
inline void VectorPartitioned::asDofs(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
dofval.fill(0.0);
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
dofval(m_dofs(m, i)) = nodevec(m, i);
}
}
}
inline void VectorPartitioned::asDofs_u(
const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 1>& dofval_u) const
{
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
dofval_u(d) = dofval(m_iiu(d));
}
}
inline void VectorPartitioned::asDofs_u(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_u) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
dofval_u.fill(0.0);
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) < m_nnu) {
dofval_u(m_part(m, i)) = nodevec(m, i);
}
}
}
}
inline void VectorPartitioned::asDofs_p(
const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 1>& dofval_p) const
{
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
dofval_p(d) = dofval(m_iip(d));
}
}
inline void VectorPartitioned::asDofs_p(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_p) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
dofval_p.fill(0.0);
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) >= m_nnu) {
dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i);
}
}
}
}
inline void VectorPartitioned::asDofs(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
dofval.fill(0.0);
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
dofval(m_dofs(m_conn(e, m), i)) = elemvec(e, m, i);
}
}
}
}
inline void VectorPartitioned::asDofs_u(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_u) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
dofval_u.fill(0.0);
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m_conn(e, m), i) < m_nnu) {
dofval_u(m_part(m_conn(e, m), i)) = elemvec(e, m, i);
}
}
}
}
}
inline void VectorPartitioned::asDofs_p(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_p) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
dofval_p.fill(0.0);
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m_conn(e, m), i) >= m_nnu) {
dofval_p(m_part(m_conn(e, m), i) - m_nnu) = elemvec(e, m, i);
}
}
}
}
}
inline void VectorPartitioned::asNode(
const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 2>& nodevec) const
{
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
nodevec(m, i) = dofval(m_dofs(m, i));
}
}
}
inline void VectorPartitioned::asNode(
const xt::xtensor<double, 1>& dofval_u,
const xt::xtensor<double, 1>& dofval_p,
xt::xtensor<double, 2>& nodevec) const
{
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) < m_nnu) {
nodevec(m, i) = dofval_u(m_part(m, i));
}
else {
nodevec(m, i) = dofval_p(m_part(m, i) - m_nnu);
}
}
}
}
inline void VectorPartitioned::asNode(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
nodevec.fill(0.0);
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
nodevec(m_conn(e, m), i) = elemvec(e, m, i);
}
}
}
}
inline void VectorPartitioned::asElement(
const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 3>& elemvec) const
{
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
elemvec(e, m, i) = dofval(m_dofs(m_conn(e, m), i));
}
}
}
}
inline void VectorPartitioned::asElement(
const xt::xtensor<double, 1>& dofval_u,
const xt::xtensor<double, 1>& dofval_p,
xt::xtensor<double, 3>& elemvec) const
{
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m_conn(e, m), i) < m_nnu) {
elemvec(e, m, i) = dofval_u(m_part(m_conn(e, m), i));
}
else {
elemvec(e, m, i) = dofval_p(m_part(m_conn(e, m), i) - m_nnu);
}
}
}
}
}
inline void VectorPartitioned::asElement(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 3>& elemvec) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
#pragma omp parallel for
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
elemvec(e, m, i) = nodevec(m_conn(e, m), i);
}
}
}
}
inline void VectorPartitioned::assembleDofs(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
dofval.fill(0.0);
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
dofval(m_dofs(m, i)) += nodevec(m, i);
}
}
}
inline void VectorPartitioned::assembleDofs_u(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_u) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
dofval_u.fill(0.0);
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) < m_nnu) {
dofval_u(m_part(m, i)) += nodevec(m, i);
}
}
}
}
inline void VectorPartitioned::assembleDofs_p(
const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_p) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
dofval_p.fill(0.0);
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) >= m_nnu) {
dofval_p(m_part(m, i) - m_nnu) += nodevec(m, i);
}
}
}
}
inline void VectorPartitioned::assembleDofs(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
dofval.fill(0.0);
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
dofval(m_dofs(m_conn(e, m), i)) += elemvec(e, m, i);
}
}
}
}
inline void VectorPartitioned::assembleDofs_u(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_u) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
dofval_u.fill(0.0);
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m_conn(e, m), i) < m_nnu) {
dofval_u(m_part(m_conn(e, m), i)) += elemvec(e, m, i);
}
}
}
}
}
inline void VectorPartitioned::assembleDofs_p(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_p) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
dofval_p.fill(0.0);
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m_conn(e, m), i) >= m_nnu) {
dofval_p(m_part(m_conn(e, m), i) - m_nnu) += elemvec(e, m, i);
}
}
}
}
}
inline void VectorPartitioned::assembleNode(
const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
{
GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
xt::xtensor<double, 1> dofval = this->AssembleDofs(elemvec);
this->asNode(dofval, nodevec);
}
inline xt::xtensor<double, 1> VectorPartitioned::AsDofs(
const xt::xtensor<double, 1>& dofval_u, const xt::xtensor<double, 1>& dofval_p) const
{
xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
this->asDofs(dofval_u, dofval_p, dofval);
return dofval;
}
inline xt::xtensor<double, 1> VectorPartitioned::AsDofs(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
this->asDofs(nodevec, dofval);
return dofval;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AsDofs_u(const xt::xtensor<double, 1>& dofval) const
{
xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
this->asDofs_u(dofval, dofval_u);
return dofval_u;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AsDofs_u(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
this->asDofs_u(nodevec, dofval_u);
return dofval_u;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AsDofs_p(const xt::xtensor<double, 1>& dofval) const
{
xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
this->asDofs_p(dofval, dofval_p);
return dofval_p;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AsDofs_p(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
this->asDofs_p(nodevec, dofval_p);
return dofval_p;
}
inline xt::xtensor<double, 1> VectorPartitioned::AsDofs(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
this->asDofs(elemvec, dofval);
return dofval;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AsDofs_u(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
this->asDofs_u(elemvec, dofval_u);
return dofval_u;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AsDofs_p(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
this->asDofs_p(elemvec, dofval_p);
return dofval_p;
}
inline xt::xtensor<double, 2> VectorPartitioned::AsNode(const xt::xtensor<double, 1>& dofval) const
{
xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
this->asNode(dofval, nodevec);
return nodevec;
}
inline xt::xtensor<double, 2> VectorPartitioned::AsNode(
const xt::xtensor<double, 1>& dofval_u, const xt::xtensor<double, 1>& dofval_p) const
{
xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
this->asNode(dofval_u, dofval_p, nodevec);
return nodevec;
}
inline xt::xtensor<double, 2> VectorPartitioned::AsNode(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
this->asNode(elemvec, nodevec);
return nodevec;
}
inline xt::xtensor<double, 3>
VectorPartitioned::AsElement(const xt::xtensor<double, 1>& dofval) const
{
xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
this->asElement(dofval, elemvec);
return elemvec;
}
inline xt::xtensor<double, 3> VectorPartitioned::AsElement(
const xt::xtensor<double, 1>& dofval_u, const xt::xtensor<double, 1>& dofval_p) const
{
xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
this->asElement(dofval_u, dofval_p, elemvec);
return elemvec;
}
inline xt::xtensor<double, 3>
VectorPartitioned::AsElement(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
this->asElement(nodevec, elemvec);
return elemvec;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AssembleDofs(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
this->assembleDofs(nodevec, dofval);
return dofval;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AssembleDofs_u(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
this->assembleDofs_u(nodevec, dofval_u);
return dofval_u;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AssembleDofs_p(const xt::xtensor<double, 2>& nodevec) const
{
xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
this->assembleDofs_p(nodevec, dofval_p);
return dofval_p;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AssembleDofs(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
this->assembleDofs(elemvec, dofval);
return dofval;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AssembleDofs_u(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
this->assembleDofs_u(elemvec, dofval_u);
return dofval_u;
}
inline xt::xtensor<double, 1>
VectorPartitioned::AssembleDofs_p(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
this->assembleDofs_p(elemvec, dofval_p);
return dofval_p;
}
inline xt::xtensor<double, 2>
VectorPartitioned::AssembleNode(const xt::xtensor<double, 3>& elemvec) const
{
xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
this->assembleNode(elemvec, nodevec);
return nodevec;
}
inline xt::xtensor<double, 2> VectorPartitioned::Copy(
const xt::xtensor<double, 2>& nodevec_src, const xt::xtensor<double, 2>& nodevec_dest) const
{
xt::xtensor<double, 2> ret = nodevec_dest;
this->copy(nodevec_src, ret);
return ret;
}
inline xt::xtensor<double, 2> VectorPartitioned::Copy_u(
const xt::xtensor<double, 2>& nodevec_src, const xt::xtensor<double, 2>& nodevec_dest) const
{
xt::xtensor<double, 2> ret = nodevec_dest;
this->copy_u(nodevec_src, ret);
return ret;
}
inline xt::xtensor<double, 2> VectorPartitioned::Copy_p(
const xt::xtensor<double, 2>& nodevec_src, const xt::xtensor<double, 2>& nodevec_dest) const
{
xt::xtensor<double, 2> ret = nodevec_dest;
this->copy_p(nodevec_src, ret);
return ret;
}
inline xt::xtensor<double, 1> VectorPartitioned::AllocateDofval() const
{
xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
return dofval;
}
inline xt::xtensor<double, 2> VectorPartitioned::AllocateNodevec() const
{
xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
return nodevec;
}
inline xt::xtensor<double, 3> VectorPartitioned::AllocateElemvec() const
{
xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
return elemvec;
}
inline xt::xtensor<double, 3> VectorPartitioned::AllocateElemmat() const
{
xt::xtensor<double, 3> elemmat = xt::empty<double>({m_nelem, m_nne * m_ndim, m_nne * m_ndim});
return elemmat;
}
inline xt::xtensor<double, 1> VectorPartitioned::AllocateDofval(double val) const
{
xt::xtensor<double, 1> dofval = xt::zeros<double>({m_ndof});
dofval.fill(val);
return dofval;
}
inline xt::xtensor<double, 2> VectorPartitioned::AllocateNodevec(double val) const
{
xt::xtensor<double, 2> nodevec = xt::zeros<double>({m_nnode, m_ndim});
nodevec.fill(val);
return nodevec;
}
inline xt::xtensor<double, 3> VectorPartitioned::AllocateElemvec(double val) const
{
xt::xtensor<double, 3> elemvec = xt::zeros<double>({m_nelem, m_nne, m_ndim});
elemvec.fill(val);
return elemvec;
}
inline xt::xtensor<double, 3> VectorPartitioned::AllocateElemmat(double val) const
{
xt::xtensor<double, 3> elemmat = xt::empty<double>({m_nelem, m_nne * m_ndim, m_nne * m_ndim});
elemmat.fill(val);
return elemmat;
}
} // namespace GooseFEM
#endif
diff --git a/python/MeshQuad4.hpp b/python/MeshQuad4.hpp
index 65d8c10..b7f38f6 100644
--- a/python/MeshQuad4.hpp
+++ b/python/MeshQuad4.hpp
@@ -1,229 +1,234 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#include <GooseFEM/GooseFEM.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pyxtensor/pyxtensor.hpp>
namespace py = pybind11;
void init_MeshQuad4(py::module& m)
{
py::class_<GooseFEM::Mesh::Quad4::Regular>(m, "Regular")
.def(
py::init<size_t, size_t, double>(),
"Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
"size 'h'",
py::arg("nx"),
py::arg("ny"),
py::arg("h") = 1.)
.def("coor", &GooseFEM::Mesh::Quad4::Regular::coor)
.def("conn", &GooseFEM::Mesh::Quad4::Regular::conn)
.def("nelem", &GooseFEM::Mesh::Quad4::Regular::nelem)
.def("nnode", &GooseFEM::Mesh::Quad4::Regular::nnode)
.def("nne", &GooseFEM::Mesh::Quad4::Regular::nne)
.def("ndim", &GooseFEM::Mesh::Quad4::Regular::ndim)
.def("nelx", &GooseFEM::Mesh::Quad4::Regular::nelx)
.def("nely", &GooseFEM::Mesh::Quad4::Regular::nely)
.def("h", &GooseFEM::Mesh::Quad4::Regular::h)
.def("getElementType", &GooseFEM::Mesh::Quad4::Regular::getElementType)
.def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge)
.def("nodesTopEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge)
.def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge)
.def("nodesRightEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge)
.def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge)
.def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge)
.def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge)
.def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge)
.def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner)
.def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner)
.def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner)
.def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner)
.def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner)
.def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner)
.def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner)
.def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner)
.def("dofs", &GooseFEM::Mesh::Quad4::Regular::dofs)
.def("nodesPeriodic", &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic)
.def("nodesOrigin", &GooseFEM::Mesh::Quad4::Regular::nodesOrigin)
.def("dofsPeriodic", &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic)
.def("elementgrid", &GooseFEM::Mesh::Quad4::Regular::elementgrid)
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) {
return "<GooseFEM.Mesh.Quad4.Regular>";
});
py::class_<GooseFEM::Mesh::Quad4::FineLayer>(m, "FineLayer")
.def(
py::init<size_t, size_t, double, size_t>(),
"FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical "
"direction",
py::arg("nx"),
py::arg("ny"),
py::arg("h") = 1.,
py::arg("nfine") = 1)
.def(
py::init<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(),
"Map connectivity to generating FineLayer-object.",
py::arg("coor"),
py::arg("conn"))
.def("coor", &GooseFEM::Mesh::Quad4::FineLayer::coor)
.def("conn", &GooseFEM::Mesh::Quad4::FineLayer::conn)
.def("nelem", &GooseFEM::Mesh::Quad4::FineLayer::nelem)
.def("nnode", &GooseFEM::Mesh::Quad4::FineLayer::nnode)
.def("nne", &GooseFEM::Mesh::Quad4::FineLayer::nne)
.def("ndim", &GooseFEM::Mesh::Quad4::FineLayer::ndim)
.def("nelx", &GooseFEM::Mesh::Quad4::FineLayer::nelx)
.def("nely", &GooseFEM::Mesh::Quad4::FineLayer::nely)
.def("h", &GooseFEM::Mesh::Quad4::FineLayer::h)
.def("elemrow_nhx", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhx)
.def("elemrow_nhy", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhy)
.def("elemrow_nelem", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nelem)
.def("getElementType", &GooseFEM::Mesh::Quad4::FineLayer::getElementType)
.def("elementsMiddleLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer)
.def("elementsLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsLayer)
.def(
"elementgrid_ravel",
&GooseFEM::Mesh::Quad4::FineLayer::elementgrid_ravel,
py::arg("rows_range"),
py::arg("cols_range"))
.def(
"elementgrid_around_ravel",
&GooseFEM::Mesh::Quad4::FineLayer::elementgrid_around_ravel,
py::arg("element"),
py::arg("size"),
py::arg("periodic") = true)
.def(
"elementgrid_leftright",
&GooseFEM::Mesh::Quad4::FineLayer::elementgrid_leftright,
py::arg("element"),
py::arg("left"),
py::arg("right"),
py::arg("periodic") = true)
.def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge)
.def("nodesTopEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge)
.def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge)
.def("nodesRightEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge)
.def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge)
.def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge)
.def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge)
.def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge)
.def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner)
.def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner)
.def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner)
.def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner)
.def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner)
.def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner)
.def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner)
.def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner)
.def("dofs", &GooseFEM::Mesh::Quad4::FineLayer::dofs)
.def("nodesPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic)
.def("nodesOrigin", &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin)
.def("dofsPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic)
.def("roll", &GooseFEM::Mesh::Quad4::FineLayer::roll)
.def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer&) {
return "<GooseFEM.Mesh.Quad4.FineLayer>";
});
}
void init_MeshQuad4Map(py::module& m)
{
py::class_<GooseFEM::Mesh::Quad4::Map::RefineRegular>(m, "RefineRegular")
.def(
py::init<const GooseFEM::Mesh::Quad4::Regular&, size_t, size_t>(),
"Refine a regular mesh",
py::arg("mesh"),
py::arg("nx"),
py::arg("ny"))
.def("getCoarseMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getCoarseMesh)
.def("getFineMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getFineMesh)
.def("getMap", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getMap)
.def(
"mapToCoarse",
py::overload_cast<const xt::xtensor<double, 1>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
.def(
"mapToCoarse",
py::overload_cast<const xt::xtensor<double, 2>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
.def(
"mapToCoarse",
py::overload_cast<const xt::xtensor<double, 4>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
.def(
"mapToFine",
py::overload_cast<const xt::xtensor<double, 1>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
.def(
"mapToFine",
py::overload_cast<const xt::xtensor<double, 2>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
.def(
"mapToFine",
py::overload_cast<const xt::xtensor<double, 4>&>(
&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::RefineRegular&) {
return "<GooseFEM.Mesh.Quad4.Map.RefineRegular>";
});
py::class_<GooseFEM::Mesh::Quad4::Map::FineLayer2Regular>(m, "FineLayer2Regular")
.def(
py::init<const GooseFEM::Mesh::Quad4::FineLayer&>(),
"Map a FineLayer mesh to a Regular mesh",
py::arg("mesh"))
.def("getRegularMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getRegularMesh)
.def("getFineLayerMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getFineLayerMesh)
.def("getMap", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMap)
.def("getMapFraction", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMapFraction)
.def(
"mapToRegular",
py::overload_cast<const xt::xtensor<double, 1>&>(
- &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
+ &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 1>, py::const_))
.def(
"mapToRegular",
py::overload_cast<const xt::xtensor<double, 2>&>(
- &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
+ &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 2>, py::const_))
+
+ .def(
+ "mapToRegular",
+ py::overload_cast<const xt::xtensor<double, 3>&>(
+ &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 3>, py::const_))
.def(
"mapToRegular",
py::overload_cast<const xt::xtensor<double, 4>&>(
- &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
+ &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 4>, py::const_))
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular&) {
return "<GooseFEM.Mesh.Quad4.Map.FineLayer2Regular>";
});
}
diff --git a/test/basic/MeshQuad4.cpp b/test/basic/MeshQuad4.cpp
index ebe1f28..aabe773 100644
--- a/test/basic/MeshQuad4.cpp
+++ b/test/basic/MeshQuad4.cpp
@@ -1,674 +1,676 @@
#include <catch2/catch.hpp>
#include <xtensor/xrandom.hpp>
#include <xtensor/xmath.hpp>
#include <GooseFEM/GooseFEM.h>
TEST_CASE("GooseFEM::MeshQuad4", "MeshQuad4.h")
{
SECTION("Regular")
{
xt::xtensor<double, 2> coor_ = {{0., 0.}, {1., 0.}, {2., 0.}, {3., 0.}, {4., 0.}, {5., 0.},
{0., 1.}, {1., 1.}, {2., 1.}, {3., 1.}, {4., 1.}, {5., 1.},
{0., 2.}, {1., 2.}, {2., 2.}, {3., 2.}, {4., 2.}, {5., 2.},
{0., 3.}, {1., 3.}, {2., 3.}, {3., 3.}, {4., 3.}, {5., 3.}};
xt::xtensor<double, 2> conn_ = {{0, 1, 7, 6},
{1, 2, 8, 7},
{2, 3, 9, 8},
{3, 4, 10, 9},
{4, 5, 11, 10},
{6, 7, 13, 12},
{7, 8, 14, 13},
{8, 9, 15, 14},
{9, 10, 16, 15},
{10, 11, 17, 16},
{12, 13, 19, 18},
{13, 14, 20, 19},
{14, 15, 21, 20},
{15, 16, 22, 21},
{16, 17, 23, 22}};
size_t nelem_ = 15;
size_t nnode_ = 24;
size_t nne_ = 4;
size_t ndim_ = 2;
size_t nnodePeriodic_ = 15;
xt::xtensor<size_t, 1> nodesBottomEdge_ = {0, 1, 2, 3, 4, 5};
xt::xtensor<size_t, 1> nodesTopEdge_ = {18, 19, 20, 21, 22, 23};
xt::xtensor<size_t, 1> nodesLeftEdge_ = {0, 6, 12, 18};
xt::xtensor<size_t, 1> nodesRightEdge_ = {5, 11, 17, 23};
xt::xtensor<size_t, 1> nodesBottomOpenEdge_ = {1, 2, 3, 4};
xt::xtensor<size_t, 1> nodesTopOpenEdge_ = {19, 20, 21, 22};
xt::xtensor<size_t, 1> nodesLeftOpenEdge_ = {6, 12};
xt::xtensor<size_t, 1> nodesRightOpenEdge_ = {11, 17};
size_t nodesBottomLeftCorner_ = 0;
size_t nodesBottomRightCorner_ = 5;
size_t nodesTopLeftCorner_ = 18;
size_t nodesTopRightCorner_ = 23;
size_t nodesLeftBottomCorner_ = 0;
size_t nodesLeftTopCorner_ = 18;
size_t nodesRightBottomCorner_ = 5;
size_t nodesRightTopCorner_ = 23;
xt::xtensor<size_t, 2> dofs_ = {{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11},
{12, 13}, {14, 15}, {16, 17}, {18, 19}, {20, 21}, {22, 23},
{24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35},
{36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}};
xt::xtensor<size_t, 2> nodesPeriodic_ = {
{0, 5}, {0, 23}, {0, 18}, {1, 19}, {2, 20}, {3, 21}, {4, 22}, {6, 11}, {12, 17}};
size_t nodesOrigin_ = 0;
xt::xtensor<size_t, 2> dofsPeriodic_ = {
{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}, {10, 11}, {12, 13},
{14, 15}, {16, 17}, {18, 19}, {10, 11}, {20, 21}, {22, 23}, {24, 25}, {26, 27},
{28, 29}, {20, 21}, {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}};
GooseFEM::Mesh::Quad4::Regular mesh(5, 3);
xt::xtensor<double, 2> coor = mesh.coor();
xt::xtensor<size_t, 2> conn = mesh.conn();
size_t nelem = mesh.nelem();
size_t nnode = mesh.nnode();
size_t nne = mesh.nne();
size_t ndim = mesh.ndim();
size_t nnodePeriodic = mesh.nnode() - mesh.nodesPeriodic().shape(0);
xt::xtensor<size_t, 1> nodesBottomEdge = mesh.nodesBottomEdge();
xt::xtensor<size_t, 1> nodesTopEdge = mesh.nodesTopEdge();
xt::xtensor<size_t, 1> nodesLeftEdge = mesh.nodesLeftEdge();
xt::xtensor<size_t, 1> nodesRightEdge = mesh.nodesRightEdge();
xt::xtensor<size_t, 1> nodesBottomOpenEdge = mesh.nodesBottomOpenEdge();
xt::xtensor<size_t, 1> nodesTopOpenEdge = mesh.nodesTopOpenEdge();
xt::xtensor<size_t, 1> nodesLeftOpenEdge = mesh.nodesLeftOpenEdge();
xt::xtensor<size_t, 1> nodesRightOpenEdge = mesh.nodesRightOpenEdge();
size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner();
size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner();
size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner();
size_t nodesTopRightCorner = mesh.nodesTopRightCorner();
size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner();
size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner();
size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner();
size_t nodesRightTopCorner = mesh.nodesRightTopCorner();
xt::xtensor<size_t, 2> dofs = mesh.dofs();
xt::xtensor<size_t, 2> nodesPeriodic = mesh.nodesPeriodic();
size_t nodesOrigin = mesh.nodesOrigin();
xt::xtensor<size_t, 2> dofsPeriodic = mesh.dofsPeriodic();
REQUIRE(nelem == nelem_);
REQUIRE(nnode == nnode_);
REQUIRE(nne == nne_);
REQUIRE(ndim == ndim_);
REQUIRE(nnodePeriodic == nnodePeriodic_);
REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_);
REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_);
REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_);
REQUIRE(nodesTopRightCorner == nodesTopRightCorner_);
REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_);
REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_);
REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_);
REQUIRE(nodesRightTopCorner == nodesRightTopCorner_);
REQUIRE(nodesOrigin == nodesOrigin_);
REQUIRE(xt::allclose(coor, coor_));
REQUIRE(xt::all(xt::equal(conn, conn_)));
REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_)));
REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_)));
REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_)));
}
SECTION("FineLayer")
{
xt::xtensor<double, 2> coor_ = {
{0., 0.}, {3., 0.}, {6., 0.}, {9., 0.}, {0., 3.}, {3., 3.}, {6., 3.}, {9., 3.},
{0., 6.}, {3., 6.}, {6., 6.}, {9., 6.}, {1., 7.}, {2., 7.}, {4., 7.}, {5., 7.},
{7., 7.}, {8., 7.}, {0., 8.}, {1., 8.}, {2., 8.}, {3., 8.}, {4., 8.}, {5., 8.},
{6., 8.}, {7., 8.}, {8., 8.}, {9., 8.}, {0., 9.}, {1., 9.}, {2., 9.}, {3., 9.},
{4., 9.}, {5., 9.}, {6., 9.}, {7., 9.}, {8., 9.}, {9., 9.}, {0., 10.}, {1., 10.},
{2., 10.}, {3., 10.}, {4., 10.}, {5., 10.}, {6., 10.}, {7., 10.}, {8., 10.}, {9., 10.},
{0., 11.}, {1., 11.}, {2., 11.}, {3., 11.}, {4., 11.}, {5., 11.}, {6., 11.}, {7., 11.},
{8., 11.}, {9., 11.}, {0., 12.}, {1., 12.}, {2., 12.}, {3., 12.}, {4., 12.}, {5., 12.},
{6., 12.}, {7., 12.}, {8., 12.}, {9., 12.}, {0., 13.}, {1., 13.}, {2., 13.}, {3., 13.},
{4., 13.}, {5., 13.}, {6., 13.}, {7., 13.}, {8., 13.}, {9., 13.}, {1., 14.}, {2., 14.},
{4., 14.}, {5., 14.}, {7., 14.}, {8., 14.}, {0., 15.}, {3., 15.}, {6., 15.}, {9., 15.},
{0., 18.}, {3., 18.}, {6., 18.}, {9., 18.}, {0., 21.}, {3., 21.}, {6., 21.}, {9., 21.}};
xt::xtensor<double, 2> conn_ = {
{0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 9, 8},
{5, 6, 10, 9}, {6, 7, 11, 10}, {8, 9, 13, 12}, {9, 21, 20, 13},
{12, 13, 20, 19}, {8, 12, 19, 18}, {9, 10, 15, 14}, {10, 24, 23, 15},
{14, 15, 23, 22}, {9, 14, 22, 21}, {10, 11, 17, 16}, {11, 27, 26, 17},
{16, 17, 26, 25}, {10, 16, 25, 24}, {18, 19, 29, 28}, {19, 20, 30, 29},
{20, 21, 31, 30}, {21, 22, 32, 31}, {22, 23, 33, 32}, {23, 24, 34, 33},
{24, 25, 35, 34}, {25, 26, 36, 35}, {26, 27, 37, 36}, {28, 29, 39, 38},
{29, 30, 40, 39}, {30, 31, 41, 40}, {31, 32, 42, 41}, {32, 33, 43, 42},
{33, 34, 44, 43}, {34, 35, 45, 44}, {35, 36, 46, 45}, {36, 37, 47, 46},
{38, 39, 49, 48}, {39, 40, 50, 49}, {40, 41, 51, 50}, {41, 42, 52, 51},
{42, 43, 53, 52}, {43, 44, 54, 53}, {44, 45, 55, 54}, {45, 46, 56, 55},
{46, 47, 57, 56}, {48, 49, 59, 58}, {49, 50, 60, 59}, {50, 51, 61, 60},
{51, 52, 62, 61}, {52, 53, 63, 62}, {53, 54, 64, 63}, {54, 55, 65, 64},
{55, 56, 66, 65}, {56, 57, 67, 66}, {58, 59, 69, 68}, {59, 60, 70, 69},
{60, 61, 71, 70}, {61, 62, 72, 71}, {62, 63, 73, 72}, {63, 64, 74, 73},
{64, 65, 75, 74}, {65, 66, 76, 75}, {66, 67, 77, 76}, {68, 69, 78, 84},
{69, 70, 79, 78}, {70, 71, 85, 79}, {78, 79, 85, 84}, {71, 72, 80, 85},
{72, 73, 81, 80}, {73, 74, 86, 81}, {80, 81, 86, 85}, {74, 75, 82, 86},
{75, 76, 83, 82}, {76, 77, 87, 83}, {82, 83, 87, 86}, {84, 85, 89, 88},
{85, 86, 90, 89}, {86, 87, 91, 90}, {88, 89, 93, 92}, {89, 90, 94, 93},
{90, 91, 95, 94}};
size_t nelem_ = 81;
size_t nnode_ = 96;
size_t nne_ = 4;
size_t ndim_ = 2;
size_t shape_x_ = 9;
size_t shape_y_ = 21;
xt::xtensor<size_t, 1> elementsMiddleLayer_ = {36, 37, 38, 39, 40, 41, 42, 43, 44};
xt::xtensor<size_t, 1> nodesBottomEdge_ = {0, 1, 2, 3};
xt::xtensor<size_t, 1> nodesTopEdge_ = {92, 93, 94, 95};
xt::xtensor<size_t, 1> nodesLeftEdge_ = {0, 4, 8, 18, 28, 38, 48, 58, 68, 84, 88, 92};
xt::xtensor<size_t, 1> nodesRightEdge_ = {3, 7, 11, 27, 37, 47, 57, 67, 77, 87, 91, 95};
xt::xtensor<size_t, 1> nodesBottomOpenEdge_ = {1, 2};
xt::xtensor<size_t, 1> nodesTopOpenEdge_ = {93, 94};
xt::xtensor<size_t, 1> nodesLeftOpenEdge_ = {4, 8, 18, 28, 38, 48, 58, 68, 84, 88};
xt::xtensor<size_t, 1> nodesRightOpenEdge_ = {7, 11, 27, 37, 47, 57, 67, 77, 87, 91};
size_t nodesBottomLeftCorner_ = 0;
size_t nodesBottomRightCorner_ = 3;
size_t nodesTopLeftCorner_ = 92;
size_t nodesTopRightCorner_ = 95;
size_t nodesLeftBottomCorner_ = 0;
size_t nodesLeftTopCorner_ = 92;
size_t nodesRightBottomCorner_ = 3;
size_t nodesRightTopCorner_ = 95;
xt::xtensor<size_t, 2> dofs_ = {
{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11}, {12, 13},
{14, 15}, {16, 17}, {18, 19}, {20, 21}, {22, 23}, {24, 25}, {26, 27},
{28, 29}, {30, 31}, {32, 33}, {34, 35}, {36, 37}, {38, 39}, {40, 41},
{42, 43}, {44, 45}, {46, 47}, {48, 49}, {50, 51}, {52, 53}, {54, 55},
{56, 57}, {58, 59}, {60, 61}, {62, 63}, {64, 65}, {66, 67}, {68, 69},
{70, 71}, {72, 73}, {74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83},
{84, 85}, {86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97},
{98, 99}, {100, 101}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111},
{112, 113}, {114, 115}, {116, 117}, {118, 119}, {120, 121}, {122, 123}, {124, 125},
{126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137}, {138, 139},
{140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149}, {150, 151}, {152, 153},
{154, 155}, {156, 157}, {158, 159}, {160, 161}, {162, 163}, {164, 165}, {166, 167},
{168, 169}, {170, 171}, {172, 173}, {174, 175}, {176, 177}, {178, 179}, {180, 181},
{182, 183}, {184, 185}, {186, 187}, {188, 189}, {190, 191}};
xt::xtensor<size_t, 2> nodesPeriodic_ = {{0, 3},
{0, 95},
{0, 92},
{1, 93},
{2, 94},
{4, 7},
{8, 11},
{18, 27},
{28, 37},
{38, 47},
{48, 57},
{58, 67},
{68, 77},
{84, 87},
{88, 91}};
size_t nodesOrigin_ = 0;
xt::xtensor<size_t, 2> dofsPeriodic_ = {
{0, 1}, {2, 3}, {4, 5}, {0, 1}, {6, 7}, {8, 9}, {10, 11},
{6, 7}, {12, 13}, {14, 15}, {16, 17}, {12, 13}, {18, 19}, {20, 21},
{22, 23}, {24, 25}, {26, 27}, {28, 29}, {30, 31}, {32, 33}, {34, 35},
{36, 37}, {38, 39}, {40, 41}, {42, 43}, {44, 45}, {46, 47}, {30, 31},
{48, 49}, {50, 51}, {52, 53}, {54, 55}, {56, 57}, {58, 59}, {60, 61},
{62, 63}, {64, 65}, {48, 49}, {66, 67}, {68, 69}, {70, 71}, {72, 73},
{74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83}, {66, 67}, {84, 85},
{86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97}, {98, 99},
{100, 101}, {84, 85}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111},
{112, 113}, {114, 115}, {116, 117}, {118, 119}, {102, 103}, {120, 121}, {122, 123},
{124, 125}, {126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137},
{120, 121}, {138, 139}, {140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149},
{150, 151}, {152, 153}, {154, 155}, {150, 151}, {156, 157}, {158, 159}, {160, 161},
{156, 157}, {0, 1}, {2, 3}, {4, 5}, {0, 1}};
GooseFEM::Mesh::Quad4::FineLayer mesh(9, 18);
xt::xtensor<double, 2> coor = mesh.coor();
xt::xtensor<size_t, 2> conn = mesh.conn();
size_t nelem = mesh.nelem();
size_t nnode = mesh.nnode();
size_t nne = mesh.nne();
size_t ndim = mesh.ndim();
size_t shape_x = mesh.nelx();
size_t shape_y = mesh.nely();
xt::xtensor<size_t, 1> elementsMiddleLayer = mesh.elementsMiddleLayer();
xt::xtensor<size_t, 1> nodesBottomEdge = mesh.nodesBottomEdge();
xt::xtensor<size_t, 1> nodesTopEdge = mesh.nodesTopEdge();
xt::xtensor<size_t, 1> nodesLeftEdge = mesh.nodesLeftEdge();
xt::xtensor<size_t, 1> nodesRightEdge = mesh.nodesRightEdge();
xt::xtensor<size_t, 1> nodesBottomOpenEdge = mesh.nodesBottomOpenEdge();
xt::xtensor<size_t, 1> nodesTopOpenEdge = mesh.nodesTopOpenEdge();
xt::xtensor<size_t, 1> nodesLeftOpenEdge = mesh.nodesLeftOpenEdge();
xt::xtensor<size_t, 1> nodesRightOpenEdge = mesh.nodesRightOpenEdge();
size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner();
size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner();
size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner();
size_t nodesTopRightCorner = mesh.nodesTopRightCorner();
size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner();
size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner();
size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner();
size_t nodesRightTopCorner = mesh.nodesRightTopCorner();
xt::xtensor<size_t, 2> dofs = mesh.dofs();
xt::xtensor<size_t, 2> nodesPeriodic = mesh.nodesPeriodic();
size_t nodesOrigin = mesh.nodesOrigin();
xt::xtensor<size_t, 2> dofsPeriodic = mesh.dofsPeriodic();
REQUIRE(nelem == nelem_);
REQUIRE(nnode == nnode_);
REQUIRE(nne == nne_);
REQUIRE(ndim == ndim_);
REQUIRE(shape_x == shape_x_);
REQUIRE(shape_y == shape_y_);
REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_);
REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_);
REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_);
REQUIRE(nodesTopRightCorner == nodesTopRightCorner_);
REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_);
REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_);
REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_);
REQUIRE(nodesRightTopCorner == nodesRightTopCorner_);
REQUIRE(nodesOrigin == nodesOrigin_);
REQUIRE(xt::allclose(coor, coor_));
REQUIRE(xt::all(xt::equal(elementsMiddleLayer, elementsMiddleLayer_)));
REQUIRE(xt::all(xt::equal(conn, conn_)));
REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_)));
REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_)));
REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_)));
REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_)));
}
SECTION("FineLayer::elementgrid_ravel - uniform")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);
xt::xtensor<size_t, 1> a = {
0, 1, 2, 3, 4,
5, 6, 7, 8, 9,
10, 11, 12, 13, 14,
15, 16, 17, 18, 19,
20, 21, 22, 23, 24};
REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({0, 5}, {0, 5}))));
REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({0, 5}, {}))));
REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({}, {0, 5}))));
REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({}, {}))));
xt::xtensor<size_t, 1> b = {
5, 6, 7, 8, 9,
10, 11, 12, 13, 14,
15, 16, 17, 18, 19};
REQUIRE(xt::all(xt::equal(b, mesh.elementgrid_ravel({1, 4}, {0, 5}))));
xt::xtensor<size_t, 1> c = {
6, 7, 8,
11, 12, 13,
16, 17, 18};
REQUIRE(xt::all(xt::equal(c, mesh.elementgrid_ravel({1, 4}, {1, 4}))));
}
SECTION("FineLayer::elementgrid_ravel - refined")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(6, 18);
xt::xtensor<size_t, 1> a = {
19, 20, 21, 22,
25, 26, 27, 28,
31, 32, 33, 34};
REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({9, 12}, {1, 5}))));
xt::xtensor<size_t, 1> b = {
0, 1,
2, 3};
REQUIRE(xt::all(xt::equal(b, mesh.elementgrid_ravel({0, 6}, {0, 6}))));
xt::xtensor<size_t, 1> c = {
50, 51,
52, 53};
REQUIRE(xt::all(xt::equal(c, mesh.elementgrid_ravel({15, 21}, {0, 6}))));
xt::xtensor<size_t, 1> d = {
4, 5, 6, 7, 8, 9, 10, 11};
REQUIRE(xt::all(xt::equal(d, mesh.elementgrid_ravel({6, 8}, {0, 6}))));
}
SECTION("FineLayer::elementgrid_ravel - uniform")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);
xt::xtensor<size_t, 1> a = {
0, 1, 2, 3, 4,
5, 6, 7, 8, 9,
10, 11, 12, 13, 14,
15, 16, 17, 18, 19,
20, 21, 22, 23, 24};
REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(10, 2)))));
REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(11, 2)))));
REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(12, 2)))));
REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(13, 2)))));
REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(14, 2)))));
xt::xtensor<size_t, 1> b10 = {
5, 6, 9,
10, 11, 14,
15, 16, 19};
xt::xtensor<size_t, 1> b11 = {
5, 6, 7,
10, 11, 12,
15, 16, 17};
xt::xtensor<size_t, 1> b12 = {
6, 7, 8,
11, 12, 13,
16, 17, 18};
xt::xtensor<size_t, 1> b13 = {
7, 8, 9,
12, 13, 14,
17, 18, 19};
xt::xtensor<size_t, 1> b14 = {
5, 8, 9,
10, 13, 14,
15, 18, 19};
REQUIRE(xt::all(xt::equal(xt::sort(b10), xt::sort(mesh.elementgrid_around_ravel(10, 1)))));
REQUIRE(xt::all(xt::equal(xt::sort(b11), xt::sort(mesh.elementgrid_around_ravel(11, 1)))));
REQUIRE(xt::all(xt::equal(xt::sort(b12), xt::sort(mesh.elementgrid_around_ravel(12, 1)))));
REQUIRE(xt::all(xt::equal(xt::sort(b13), xt::sort(mesh.elementgrid_around_ravel(13, 1)))));
REQUIRE(xt::all(xt::equal(xt::sort(b14), xt::sort(mesh.elementgrid_around_ravel(14, 1)))));
REQUIRE(xt::all(xt::equal(10, xt::sort(mesh.elementgrid_around_ravel(10, 0)))));
REQUIRE(xt::all(xt::equal(11, xt::sort(mesh.elementgrid_around_ravel(11, 0)))));
REQUIRE(xt::all(xt::equal(12, xt::sort(mesh.elementgrid_around_ravel(12, 0)))));
REQUIRE(xt::all(xt::equal(13, xt::sort(mesh.elementgrid_around_ravel(13, 0)))));
REQUIRE(xt::all(xt::equal(14, xt::sort(mesh.elementgrid_around_ravel(14, 0)))));
}
SECTION("FineLayer::elementgrid_around_ravel")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(12, 18);
xt::xtensor<size_t, 1> r0 = {
36, 37, 71,
48, 49, 59,
60, 61, 47,
};
xt::xtensor<size_t, 1> r1 = {
36, 37, 38,
48, 49, 50,
60, 61, 62,
};
xt::xtensor<size_t, 1> r12 = {
46, 47, 36,
58, 59, 48,
70, 71, 60,
};
REQUIRE(xt::all(xt::equal(xt::sort(r0), xt::sort(mesh.elementgrid_around_ravel(48, 1)))));
for (size_t n = 0; n < 10; ++n) {
REQUIRE(xt::all(xt::equal(xt::sort(r1) + n, xt::sort(mesh.elementgrid_around_ravel(49 + n, 1)))));
}
REQUIRE(xt::all(xt::equal(xt::sort(r12), xt::sort(mesh.elementgrid_around_ravel(59, 1)))));
}
SECTION("FineLayer::elementgrid_leftright")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(12, 18);
xt::xtensor<size_t, 1> r0 = {
48, 49, 59,
};
xt::xtensor<size_t, 1> r1 = {
48, 49, 50,
};
xt::xtensor<size_t, 1> r12 = {
58, 59, 48,
};
REQUIRE(xt::all(xt::equal(xt::sort(r0), xt::sort(mesh.elementgrid_leftright(48, 1, 1)))));
for (size_t n = 0; n < 10; ++n) {
REQUIRE(xt::all(xt::equal(xt::sort(r1) + n, xt::sort(mesh.elementgrid_leftright(49 + n, 1, 1)))));
}
REQUIRE(xt::all(xt::equal(xt::sort(r12), xt::sort(mesh.elementgrid_leftright(59, 1, 1)))));
}
SECTION("FineLayer - replica - trivial")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(1, 1);
GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
}
SECTION("FineLayer - replica - equi-sized")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(4, 4);
GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
}
SECTION("FineLayer - replica")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27);
GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
}
SECTION("FineLayer - replica - one")
{
double h = xt::numeric_constants<double>::PI;
GooseFEM::Mesh::Quad4::FineLayer mesh(1, 1, h);
GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
}
SECTION("FineLayer - replica - large")
{
size_t N = 2 * std::pow(3, 6);
double h = xt::numeric_constants<double>::PI;
GooseFEM::Mesh::Quad4::FineLayer mesh(N, N, h);
GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
}
SECTION("FineLayer - roll - trivial")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3);
xt::xtensor<size_t, 1> m0 = {
0, 1, 2,
3, 4, 5,
6, 7, 8
};
xt::xtensor<size_t, 1> m1 = {
2, 0, 1,
5, 3, 4,
8, 6, 7
};
xt::xtensor<size_t, 1> m2 = {
1, 2, 0,
4, 5, 3,
7, 8, 6,
};
REQUIRE(xt::all(xt::equal(mesh.roll(0), m0)));
REQUIRE(xt::all(xt::equal(mesh.roll(1), m1)));
REQUIRE(xt::all(xt::equal(mesh.roll(2), m2)));
REQUIRE(xt::all(xt::equal(mesh.roll(3), m0)));
REQUIRE(xt::all(xt::equal(mesh.roll(4), m1)));
REQUIRE(xt::all(xt::equal(mesh.roll(5), m2)));
}
SECTION("FineLayer - roll - a")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(6, 18);
xt::xtensor<size_t, 1> m0 = {
0, 1,
2, 3,
4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49,
50, 51,
52, 53
};
xt::xtensor<size_t, 1> m1 = {
1, 0,
3, 2,
8, 9, 10, 11, 4, 5, 6, 7,
15, 16, 17, 12, 13, 14,
21, 22, 23, 18, 19, 20,
27, 28, 29, 24, 25, 26,
33, 34, 35, 30, 31, 32,
39, 40, 41, 36, 37, 38,
46, 47, 48, 49, 42, 43, 44, 45,
51, 50,
53, 52
};
REQUIRE(xt::all(xt::equal(mesh.roll(0), m0)));
REQUIRE(xt::all(xt::equal(mesh.roll(1), m1)));
REQUIRE(xt::all(xt::equal(mesh.roll(2), m0)));
REQUIRE(xt::all(xt::equal(mesh.roll(3), m1)));
}
SECTION("FineLayer - roll - b")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(9, 18);
xt::xtensor<size_t, 1> m0 = {
0, 1, 2,
3, 4, 5,
6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50, 51, 52, 53,
54, 55, 56, 57, 58, 59, 60, 61, 62,
63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74,
75, 76, 77,
78, 79, 80
};
xt::xtensor<size_t, 1> m1 = {
2, 0, 1,
5, 3, 4,
14, 15, 16, 17, 6, 7, 8, 9, 10, 11, 12, 13,
24, 25, 26, 18, 19, 20, 21, 22, 23,
33, 34, 35, 27, 28, 29, 30, 31, 32,
42, 43, 44, 36, 37, 38, 39, 40, 41,
51, 52, 53, 45, 46, 47, 48, 49, 50,
60, 61, 62, 54, 55, 56, 57, 58, 59,
71, 72, 73, 74, 63, 64, 65, 66, 67, 68, 69, 70,
77, 75, 76,
80, 78, 79
};
xt::xtensor<size_t, 1> m2 = {
1, 2, 0,
4, 5, 3,
10, 11, 12, 13, 14, 15, 16, 17, 6, 7, 8, 9,
21, 22, 23, 24, 25, 26, 18, 19, 20,
30, 31, 32, 33, 34, 35, 27, 28, 29,
39, 40, 41, 42, 43, 44, 36, 37, 38,
48, 49, 50, 51, 52, 53, 45, 46, 47,
57, 58, 59, 60, 61, 62, 54, 55, 56,
67, 68, 69, 70, 71, 72, 73, 74, 63, 64, 65, 66,
76, 77, 75,
79, 80, 78
};
REQUIRE(xt::all(xt::equal(mesh.roll(0), m0)));
REQUIRE(xt::all(xt::equal(mesh.roll(1), m1)));
REQUIRE(xt::all(xt::equal(mesh.roll(2), m2)));
REQUIRE(xt::all(xt::equal(mesh.roll(3), m0)));
REQUIRE(xt::all(xt::equal(mesh.roll(4), m1)));
REQUIRE(xt::all(xt::equal(mesh.roll(5), m2)));
}
SECTION("Map::RefineRegular")
{
GooseFEM::Mesh::Quad4::Regular mesh(5, 4);
GooseFEM::Mesh::Quad4::Map::RefineRegular refine(mesh, 5, 3);
- xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nelem()});
+ std::array<size_t, 1> as = {mesh.nelem()};
+ xt::xtensor<double, 1> a = xt::random::rand<double>(as);
auto a_ = refine.mapToCoarse(refine.mapToFine(a));
REQUIRE(xt::allclose(a, xt::mean(a_, {1})));
- xt::xtensor<double, 2> b =
- xt::random::rand<double>(std::array<size_t, 2>{mesh.nelem(), 4ul});
+ std::array<size_t, 2> bs = {mesh.nelem(), 4};
+ xt::xtensor<double, 2> b = xt::random::rand<double>(bs);
auto b_ = refine.mapToCoarse(refine.mapToFine(b));
REQUIRE(xt::allclose(xt::mean(b, {1}), xt::mean(b_, {1})));
- xt::xtensor<double, 4> c =
- xt::random::rand<double>(std::array<size_t, 4>{mesh.nelem(), 4ul, 3ul, 3ul});
+ std::array<size_t, 4> cs = {mesh.nelem(), 4, 3, 3};
+ xt::xtensor<double, 4> c = xt::random::rand<double>(cs);
auto c_ = refine.mapToCoarse(refine.mapToFine(c));
REQUIRE(xt::allclose(xt::mean(c, {1}), xt::mean(c_, {1})));
}
SECTION("Map::FineLayer2Regular")
{
GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);
GooseFEM::Mesh::Quad4::Map::FineLayer2Regular map(mesh);
- xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nelem()});
+ std::array<size_t, 1> as = {mesh.nelem()};
+ xt::xtensor<double, 1> a = xt::random::rand<double>(as);
auto a_ = map.mapToRegular(a);
REQUIRE(xt::allclose(a, a_));
- xt::xtensor<double, 2> b =
- xt::random::rand<double>(std::array<size_t, 2>{mesh.nelem(), 4ul});
+ std::array<size_t, 2> bs = {mesh.nelem(), 4};
+ xt::xtensor<double, 2> b = xt::random::rand<double>(bs);
auto b_ = map.mapToRegular(b);
REQUIRE(xt::allclose(b, b_));
- xt::xtensor<double, 4> c =
- xt::random::rand<double>(std::array<size_t, 4>{mesh.nelem(), 4ul, 3ul, 3ul});
+ std::array<size_t, 4> cs = {mesh.nelem(), 4, 3, 3};
+ xt::xtensor<double, 4> c = xt::random::rand<double>(cs);
auto c_ = map.mapToRegular(c);
REQUIRE(xt::allclose(c, c_));
}
}