diff --git a/docs/details/MeshQuad4.rst b/docs/details/MeshQuad4.rst index 303b430..cb105bb 100644 --- a/docs/details/MeshQuad4.rst +++ b/docs/details/MeshQuad4.rst @@ -1,285 +1,292 @@ .. _MeshQuad4: *********** Mesh::Quad4 *********** | :download:`GooseFEM/MeshQuad4.h <../../include/GooseFEM/MeshQuad4.h>` | :download:`GooseFEM/MeshQuad4.hpp <../../include/GooseFEM/MeshQuad4.hpp>` Naming convention ================= .. image:: figures/MeshQuad4/naming_convention.svg :width: 350px :align: center Mesh::Quad4::Regular ==================== Regular mesh. .. seealso:: | :download:`Python - example ` Mesh::Quad4::Regular::nelem() ----------------------------- Return number of elements. Mesh::Quad4::Regular::nnode() ----------------------------- Return number of nodes. Mesh::Quad4::Regular::nne() --------------------------- Return number of nodes-per-element (= 3). Mesh::Quad4::Regular::ndim() ---------------------------- Return number of dimensions (= 2). Mesh::Quad4::Regular::getElementType() -------------------------------------- Return element-type. Mesh::Quad4::Regular::coor() ---------------------------- Return nodal coordinates [nnode, ndim]. Mesh::Quad4::Regular::conn() ---------------------------- Return connectivity [nelem, nne]. Mesh::Quad4::Regular::nodesXXXEdge() ------------------------------------ Node numbers along the "Bottom", "Top", "Left", or "Right" edge. Mesh::Quad4::Regular::nodesXXXOpenEdge() ---------------------------------------- Node numbers along the "Bottom", "Top", "Left", or "Right" edge, excluding the corners. Mesh::Quad4::Regular::nodesXXXCorner() -------------------------------------- Node number of one of the corners (e.g. "BottomLeft"). Mesh::Quad4::Regular::nodesPeriodic() ------------------------------------- Periodic node pairs. Each row contains on pair of (independent, dependent) node numbers. The output shape is thus [n_pairs, 2]. Mesh::Quad4::Regular::nodesOrigin() ----------------------------------- Bottom-left node, used as reference for periodicity. Mesh::Quad4::Regular::dofs() ---------------------------- DOF-numbers for each component of each node (sequential). The output shape is thus [nnode, ndim]. Mesh::Quad4::Regular::dofsPeriodic() ------------------------------------ DOF-numbers for each component of each node, for the case that the periodicity if fully eliminated. The output shape is thus [nnode, ndim]. Mesh::Quad4::Regular::elementMatrix() ------------------------------------- Return element numbers as matrix [nely, nelx]. Mesh::Quad4::FineLayer ====================== Mesh with a fine layer in the middle, and that becomes course away from this plane (see image below). Note coursening can only be done if the number of elements in horizontal direction is dividable by 3, and that it is only optimal if the number of elements in horizontal direction is a factor of 3. Note that the number of elements in the vertical direction is specified as the number of times the unit element (the number of times "h" the height should be), and that this number is only a target: the algorithm chooses in accordance with the applied coursing. .. image:: figures/MeshQuad4/FineLayer/behaviour.svg :width: 800px :align: center .. seealso:: | :download:`Python - example ` | :download:`Python - behaviour 'nx' ` | :download:`Python - element numbers ` Mesh::Quad4::FineLayer::nelem() ------------------------------- Return number of elements. Mesh::Quad4::FineLayer::nnode() ------------------------------- Return number of nodes. Mesh::Quad4::FineLayer::nne() ----------------------------- Return number of nodes-per-element (= 3). Mesh::Quad4::FineLayer::ndim() ------------------------------ Return number of dimensions (= 2). Mesh::Quad4::FineLayer::nelx() ------------------------------ Number of elements in horizontal direction (along the weak layer) (matches input). Mesh::Quad4::FineLayer::nely() ------------------------------ Actual number of elements unit elements in vertical direction (actual number of times "h" the mesh is heigh). Mesh::Quad4::FineLayer::h() --------------------------- Unit edge size (matches input). Mesh::Quad4::FineLayer::getElementType() ---------------------------------------- Return element-type. Mesh::Quad4::FineLayer::coor() ------------------------------ Return nodal coordinates [nnode, ndim]. Mesh::Quad4::FineLayer::conn() ------------------------------ Return connectivity [nelem, nne]. Mesh::Quad4::FineLayer::nodesXXXEdge() -------------------------------------- Node numbers along the "Bottom", "Top", "Left", or "Right" edge. Mesh::Quad4::FineLayer::nodesXXXOpenEdge() ------------------------------------------ Node numbers along the "Bottom", "Top", "Left", or "Right" edge, excluding the corners. Mesh::Quad4::FineLayer::nodesXXXCorner() ---------------------------------------- Node number of one of the corners (e.g. "BottomLeft"). Mesh::Quad4::FineLayer::nodesPeriodic() --------------------------------------- Periodic node pairs. Each row contains on pair of (independent, dependent) node numbers. The output shape is thus [n_pairs, 2]. Mesh::Quad4::FineLayer::nodesOrigin() ------------------------------------- Bottom-left node, used as reference for periodicity. Mesh::Quad4::FineLayer::dofs() ------------------------------ DOF-numbers for each component of each node (sequential). The output shape is thus [nnode, ndim]. Mesh::Quad4::FineLayer::dofsPeriodic() -------------------------------------- DOF-numbers for each component of each node, for the case that the periodicity if fully eliminated. The output shape is thus [nnode, ndim]. Mesh::Quad4::FineLayer::elementsMiddleLayer() --------------------------------------------- Element numbers of the middle, fine, layer +Details +------- + +.. image:: figures/MeshQuad4/FineLayer/element_numbers.svg + :width: 400px + :align: center + Mesh::Quad4::Map::RefineRegular =============================== Refine a "Regular" mesh. Mesh::Quad4::Map::RefineRegular::getCoarseMesh() ------------------------------------------------ Return course mesh as "Mesh::Quad4::Regular". Mesh::Quad4::Map::RefineRegular::getFineMesh() ---------------------------------------------- Return fine mesh as "Mesh::Quad4::Regular". Mesh::Quad4::Map::RefineRegular::getMap() ----------------------------------------- Elements of the fine mesh per element of the coarse mesh (rows). Mesh::Quad4::Map::RefineRegular::mapToCoarse(...) ------------------------------------------------- Map field to the course mesh: * Scalar per element. * Scalar per integration point. * Tensor per integration point. Mesh::Quad4::Map::RefineRegular::mapToFine(...) ----------------------------------------------- Map field to the fine mesh: * Scalar per element. * Scalar per integration point. * Tensor per integration point. Mesh::Quad4::Map::FineLayer2Regular =================================== Map "Regular" mesh to "FineLayer" mesh. .. image:: figures/MeshQuad4/Map/FineLayer2Regular/map.svg :width: 350px :align: center .. seealso:: | :download:`Python - map ` | :download:`Python - element numbers ` Mesh::Quad4::Map::FineLayer2Regular::getCoarseMesh() ---------------------------------------------------- Return course mesh as "Mesh::Quad4::Regular". Mesh::Quad4::Map::FineLayer2Regular::getFineMesh() -------------------------------------------------- Return fine mesh as "Mesh::Quad4::Regular". Mesh::Quad4::Map::FineLayer2Regular::getMap() --------------------------------------------- Elements of the fine mesh per element of the coarse mesh (rows). Mesh::Quad4::Map::FineLayer2Regular::getMapFraction() ----------------------------------------------------- Get the fraction of overlap for the output of "getMap()". Mesh::Quad4::Map::FineLayer2Regular::mapToRegular(...) ------------------------------------------------------ Map field to the course mesh: * Scalar per element. * Scalar per integration point. * Tensor per integration point. diff --git a/docs/details/figures/MeshQuad4/FineLayer/element_numbers.svg b/docs/details/figures/MeshQuad4/FineLayer/element_numbers.svg new file mode 100644 index 0000000..8f996d4 --- /dev/null +++ b/docs/details/figures/MeshQuad4/FineLayer/element_numbers.svg @@ -0,0 +1,1021 @@ + + + + + + + + + 2020-12-09T17:29:49.371989 + image/svg+xml + + + Matplotlib v3.3.2, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h index 000e9d7..1494415 100644 --- a/include/GooseFEM/MeshQuad4.h +++ b/include/GooseFEM/MeshQuad4.h @@ -1,254 +1,258 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_MESHQUAD4_H #define GOOSEFEM_MESHQUAD4_H #include "config.h" namespace GooseFEM { namespace Mesh { namespace Quad4 { namespace Map { class FineLayer2Regular; } // namespace Map 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 // type ElementType getElementType() const; // mesh xt::xtensor coor() const; // nodal positions [nnode, ndim] xt::xtensor conn() const; // connectivity [nelem, nne] // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor 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 dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; // element numbers as matrix xt::xtensor elementMatrix() 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 }; 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& coor, const xt::xtensor& 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 // type ElementType getElementType() const; // mesh xt::xtensor coor() const; // nodal positions [nnode, ndim] xt::xtensor conn() const; // connectivity [nelem, nne] // element sets xt::xtensor elementsMiddleLayer() const; // elements in the middle (fine) layer // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor 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 dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; + // mapping to 'roll' periodically, + // return element mapping, to get the new connectivity use "comm[mesh.roll(), :]" + xt::xtensor roll(size_t n, size_t axis = 0); + 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 m_nelx; // number of elements in "x" (*) xt::xtensor m_nnd; // total number of nodes in the main node layer (**) xt::xtensor m_nhx; // element size in x-direction (*) xt::xtensor m_nhy; // element size in y-direction (*) xt::xtensor m_refine; // refine direction (-1:no refine, 0:"x" (*) xt::xtensor m_startElem; // start element (*) xt::xtensor 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& coor, const xt::xtensor& conn); friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular; }; namespace Map { // Return "FineLayer"-class responsible for generating a connectivity // Throws if conversion is not possible GooseFEM::Mesh::Quad4::FineLayer FineLayer( const xt::xtensor& coor, const xt::xtensor& conn); class RefineRegular { public: // constructor RefineRegular() = default; RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny); // return the one of the two meshes 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 getMap() const; // map field xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToCoarse(const xt::xtensor& data) const; // tensor per intpnt // map field xt::xtensor mapToFine(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToFine(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToFine(const xt::xtensor& data) const; // tensor per intpnt private: // the meshes GooseFEM::Mesh::Quad4::Regular m_coarse; GooseFEM::Mesh::Quad4::Regular m_fine; // mapping xt::xtensor m_fine2coarse; xt::xtensor m_fine2coarse_index; xt::xtensor 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> getMap() const; std::vector> getMapFraction() const; // map field xt::xtensor mapToRegular(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToRegular(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToRegular(const xt::xtensor& data) const; // tensor per intpnt 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> m_elem_regular; std::vector> 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 0fc8465..b6d151b 100644 --- a/include/GooseFEM/MeshQuad4.hpp +++ b/include/GooseFEM/MeshQuad4.hpp @@ -1,1383 +1,1411 @@ /* (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 Regular::coor() const { xt::xtensor ret = xt::empty({m_nnode, m_ndim}); xt::xtensor x = xt::linspace(0.0, m_h * static_cast(m_nelx), m_nelx + 1); xt::xtensor y = xt::linspace(0.0, m_h * static_cast(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 Regular::conn() const { xt::xtensor ret = xt::empty({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 Regular::nodesBottomEdge() const { return xt::arange(m_nelx + 1); } inline xt::xtensor Regular::nodesTopEdge() const { return xt::arange(m_nelx + 1) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1) + m_nelx; } inline xt::xtensor Regular::nodesBottomOpenEdge() const { return xt::arange(1, m_nelx); } inline xt::xtensor Regular::nodesTopOpenEdge() const { return xt::arange(1, m_nelx) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftOpenEdge() const { return xt::arange(1, m_nely) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightOpenEdge() const { return xt::arange(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 Regular::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(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 Regular::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor Regular::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor 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 Regular::elementMatrix() const { return xt::arange(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& coor, const xt::xtensor& 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(nelx); // compute element size in y-direction (use symmetry, compute upper half) // temporary variables size_t nmin, ntot; xt::xtensor nhx = xt::ones({nely}); xt::xtensor nhy = xt::ones({nely}); xt::xtensor refine = -1 * xt::ones({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({nely * 2 - 1}); m_nhy = xt::empty({nely * 2 - 1}); m_refine = xt::empty({nely * 2 - 1}); m_nelx = xt::empty({nely * 2 - 1}); m_nnd = xt::empty({nely * 2}); m_startElem = xt::empty({nely * 2 - 1}); m_startNode = xt::empty({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 ElementType FineLayer::getElementType() const { return ElementType::Quad4; } inline xt::xtensor FineLayer::coor() const { // allocate output xt::xtensor ret = xt::empty({m_nnode, m_ndim}); // current node, number of element layers size_t inode = 0; size_t nely = static_cast(m_nhy.size()); // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) // - allocate xt::xtensor y = xt::empty({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 x = xt::linspace(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(m_nhx(iy) / 3); double dy = m_h * static_cast(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(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 x = xt::linspace(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(m_nhx(iy) / 3); double dy = m_h * static_cast(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(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 FineLayer::conn() const { // allocate output xt::xtensor ret = xt::empty({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(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 FineLayer::elementsMiddleLayer() const { size_t nely = m_nhy.size(); size_t iy = (nely - 1) / 2; return m_startElem(iy) + xt::arange(m_nelx(iy)); } inline xt::xtensor FineLayer::nodesBottomEdge() const { return m_startNode(0) + xt::arange(m_nelx(0) + 1); } inline xt::xtensor FineLayer::nodesTopEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(m_nelx(nely - 1) + 1); } inline xt::xtensor FineLayer::nodesLeftEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesRightEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesBottomOpenEdge() const { return m_startNode(0) + xt::arange(1, m_nelx(0)); } inline xt::xtensor FineLayer::nodesTopOpenEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(1, m_nelx(nely - 1)); } inline xt::xtensor FineLayer::nodesLeftOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesRightOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({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 FineLayer::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(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 FineLayer::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor FineLayer::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor 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 FineLayer::roll(size_t n, size_t axis) +{ + GOOSEFEM_ASSERT(axis == 0); + auto conn = this->conn(); + size_t nely = static_cast(m_nhy.size()); + xt::xtensor ret = xt::empty({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(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& coor, const xt::xtensor& 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)); auto n0 = xt::view(conn, xt::all(), 0); auto n1 = xt::view(conn, xt::all(), 1); auto n2 = xt::view(conn, xt::all(), 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(xt::where(dx > 0.0, dx, std::numeric_limits::max()))(); auto hy = xt::amin(xt::where(dy > 0.0, dx, std::numeric_limits::max()))(); GOOSEFEM_CHECK(xt::allclose(hx, hy)); if (conn.shape(0) == 1) { this->init(1, 1, hx); GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn))); GOOSEFEM_CHECK(xt::allclose(this->coor(), coor)); return; } 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)); size_t nelx = eright - eleft + 1; size_t nely = static_cast((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(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 elmat_coarse = m_coarse.elementMatrix(); xt::xtensor elmat_fine = m_fine.elementMatrix(); m_coarse2fine = xt::empty({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 RefineRegular::getMap() const { return m_coarse2fine; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& 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 ret = xt::empty({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 RefineRegular::mapToCoarse(const xt::xtensor& 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 ret = xt::empty({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 RefineRegular::mapToCoarse(const xt::xtensor& 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 ret = xt::empty({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 RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({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 RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({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 RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({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 nhx = m_finelayer.m_nhx; xt::xtensor nhy = m_finelayer.m_nhy; xt::xtensor nelx = m_finelayer.m_nelx; xt::xtensor start = m_finelayer.m_startElem; // 'matrix' of element numbers of the Regular-mesh xt::xtensor elementMatrix = m_regular.elementMatrix(); // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh xt::xtensor cum_nhy = xt::concatenate(xt::xtuple(xt::zeros({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(elementMatrix, 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 el_old = start(iy) + xt::arange(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 el_old = start(iy) + xt::arange(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 el_old = start(iy) + xt::arange(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> FineLayer2Regular::getMap() const { return m_elem_regular; } inline std::vector> FineLayer2Regular::getMapFraction() const { return m_frac_regular; } inline xt::xtensor FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); xt::xtensor ret = xt::zeros({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 FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); xt::xtensor ret = xt::zeros({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 FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); xt::xtensor ret = xt::zeros({m_regular.nelem(), data.shape(1), data.shape(2), data.shape(3)}); 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/test/basic/MeshQuad4.cpp b/test/basic/MeshQuad4.cpp index e1e89ea..c858a82 100644 --- a/test/basic/MeshQuad4.cpp +++ b/test/basic/MeshQuad4.cpp @@ -1,396 +1,507 @@ #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::MeshQuad4", "MeshQuad4.h") { SECTION("Regular") { xt::xtensor 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 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 nodesBottomEdge_ = {0, 1, 2, 3, 4, 5}; xt::xtensor nodesTopEdge_ = {18, 19, 20, 21, 22, 23}; xt::xtensor nodesLeftEdge_ = {0, 6, 12, 18}; xt::xtensor nodesRightEdge_ = {5, 11, 17, 23}; xt::xtensor nodesBottomOpenEdge_ = {1, 2, 3, 4}; xt::xtensor nodesTopOpenEdge_ = {19, 20, 21, 22}; xt::xtensor nodesLeftOpenEdge_ = {6, 12}; xt::xtensor 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 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 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 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 coor = mesh.coor(); xt::xtensor 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 nodesBottomEdge = mesh.nodesBottomEdge(); xt::xtensor nodesTopEdge = mesh.nodesTopEdge(); xt::xtensor nodesLeftEdge = mesh.nodesLeftEdge(); xt::xtensor nodesRightEdge = mesh.nodesRightEdge(); xt::xtensor nodesBottomOpenEdge = mesh.nodesBottomOpenEdge(); xt::xtensor nodesTopOpenEdge = mesh.nodesTopOpenEdge(); xt::xtensor nodesLeftOpenEdge = mesh.nodesLeftOpenEdge(); xt::xtensor 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 dofs = mesh.dofs(); xt::xtensor nodesPeriodic = mesh.nodesPeriodic(); size_t nodesOrigin = mesh.nodesOrigin(); xt::xtensor 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 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 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 elementsMiddleLayer_ = {36, 37, 38, 39, 40, 41, 42, 43, 44}; xt::xtensor nodesBottomEdge_ = {0, 1, 2, 3}; xt::xtensor nodesTopEdge_ = {92, 93, 94, 95}; xt::xtensor nodesLeftEdge_ = {0, 4, 8, 18, 28, 38, 48, 58, 68, 84, 88, 92}; xt::xtensor nodesRightEdge_ = {3, 7, 11, 27, 37, 47, 57, 67, 77, 87, 91, 95}; xt::xtensor nodesBottomOpenEdge_ = {1, 2}; xt::xtensor nodesTopOpenEdge_ = {93, 94}; xt::xtensor nodesLeftOpenEdge_ = {4, 8, 18, 28, 38, 48, 58, 68, 84, 88}; xt::xtensor 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 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 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 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 coor = mesh.coor(); xt::xtensor 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 elementsMiddleLayer = mesh.elementsMiddleLayer(); xt::xtensor nodesBottomEdge = mesh.nodesBottomEdge(); xt::xtensor nodesTopEdge = mesh.nodesTopEdge(); xt::xtensor nodesLeftEdge = mesh.nodesLeftEdge(); xt::xtensor nodesRightEdge = mesh.nodesRightEdge(); xt::xtensor nodesBottomOpenEdge = mesh.nodesBottomOpenEdge(); xt::xtensor nodesTopOpenEdge = mesh.nodesTopOpenEdge(); xt::xtensor nodesLeftOpenEdge = mesh.nodesLeftOpenEdge(); xt::xtensor 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 dofs = mesh.dofs(); xt::xtensor nodesPeriodic = mesh.nodesPeriodic(); size_t nodesOrigin = mesh.nodesOrigin(); xt::xtensor 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 - 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 - roll - trivial") + { + GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3); + xt::xtensor m0 = { + 0, 1, 2, + 3, 4, 5, + 6, 7, 8 + }; + xt::xtensor m1 = { + 2, 0, 1, + 5, 3, 4, + 8, 6, 7 + }; + xt::xtensor 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 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 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 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 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 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 a = xt::random::rand({mesh.nelem()}); auto a_ = refine.mapToCoarse(refine.mapToFine(a)); REQUIRE(xt::allclose(a, xt::mean(a_, {1}))); xt::xtensor b = xt::random::rand(std::array{mesh.nelem(), 4ul}); auto b_ = refine.mapToCoarse(refine.mapToFine(b)); REQUIRE(xt::allclose(xt::mean(b, {1}), xt::mean(b_, {1}))); xt::xtensor c = xt::random::rand(std::array{mesh.nelem(), 4ul, 3ul, 3ul}); 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 a = xt::random::rand({mesh.nelem()}); auto a_ = map.mapToRegular(a); REQUIRE(xt::allclose(a, a_)); xt::xtensor b = xt::random::rand(std::array{mesh.nelem(), 4ul}); auto b_ = map.mapToRegular(b); REQUIRE(xt::allclose(b, b_)); xt::xtensor c = xt::random::rand(std::array{mesh.nelem(), 4ul, 3ul, 3ul}); auto c_ = map.mapToRegular(c); REQUIRE(xt::allclose(c, c_)); } }