diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h index f50aabe..ffb5d84 100644 --- a/include/GooseFEM/MeshQuad4.h +++ b/include/GooseFEM/MeshQuad4.h @@ -1,258 +1,269 @@ /* (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 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 }; 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 + // - elements in the middle (fine) layer + xt::xtensor elementsMiddleLayer() const; + // - select region of elements from 'matrix' of element numbers + xt::xtensor elementgrid_ravel( + std::array start_stop_rows, + std::array start_stop_cols) const; + // - select region of elements from 'matrix' of element numbers around an element + xt::xtensor elementgrid_around_ravel( + size_t element, + std::array start_stop_rows, + std::array start_stop_cols) const; + // xt::xtensor elementsAround(size_t e, int left, int right) const; // region-of-interests // 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 in the x-direction, // returns element mapping, such that: new_elemvar = elemvar[elem_map] xt::xtensor 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 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 f3f3b7e..620dd35 100644 --- a/include/GooseFEM/MeshQuad4.hpp +++ b/include/GooseFEM/MeshQuad4.hpp @@ -1,1418 +1,1558 @@ /* (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::elementgrid() 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::elementgrid_ravel( + std::array start_stop_rows, + std::array start_stop_cols) const +{ + GOOSEFEM_ASSERT(start_stop_rows[0] <= this->nely()); + GOOSEFEM_ASSERT(start_stop_rows[1] <= this->nely()); + GOOSEFEM_ASSERT(start_stop_cols[0] <= this->nelx()); + GOOSEFEM_ASSERT(start_stop_cols[1] <= this->nelx()); + + if (start_stop_rows[0] == start_stop_rows[1] || start_stop_cols[0] == start_stop_cols[1]) { + xt::xtensor ret = xt::empty({0}); + return ret; + } + + // Compute dimensions + + auto H = xt::cumsum(m_nhy); + size_t yl = 0; + if (start_stop_rows[0] > 0) { + yl = xt::argmax(H > start_stop_rows[0])(); + } + size_t yu = xt::argmax(H >= start_stop_rows[1])(); + size_t hx = std::max(m_nhx(yl), m_nhx(yu)); + size_t xl = (start_stop_cols[0] - start_stop_cols[0] % hx) / hx; + size_t xu = (start_stop_cols[1] - start_stop_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 ret = xt::empty({N}); + + // Write output + + 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; + } + xt::xtensor e = m_startElem(i) + xl * hx / m_nhx(i) + xt::arange(n); + xt::view(ret, xt::range(N, N + n)) = e; + N += n; + } + + return ret; +} + + +// inline xt::xtensor FineLayer::elementsAround(size_t e, int left, int right) const +// { +// GOOSEFEM_ASSERT(right - left <= static_cast(this->nelx())); +// GOOSEFEM_ASSERT(right - left <= static_cast(this->nely())); +// GOOSEFEM_ASSERT(e < this->nelem()); +// GOOSEFEM_WIP_ASSERT(left <= 0); +// GOOSEFEM_WIP_ASSERT(right >= 0); + +// size_t l = static_cast(std::abs(left)); +// size_t r = static_cast(std::abs(right)); +// size_t w = l + r; + +// size_t nely = m_nhy.size(); +// size_t iy = xt::argmin(xt::greater(e, m_startElem))() - 1; +// size_t nel = m_nelx(iy); +// size_t mid = static_cast(static_cast(m_startElem(iy)) - left); +// size_t nroll = nel - (nel - nel % 2) / 2 + e - m_startElem(iy); +// GOOSEFEM_WIP_ASSERT(iy == (nely - 1) / 2); +// GOOSEFEM_WIP_ASSERT(iy >= l); +// GOOSEFEM_WIP_ASSERT(iy < r); + +// // Select elements layers around the middle layer + +// size_t nyb = 0; +// size_t nyt = 0; +// size_t nhb = 0; +// size_t nht = m_nhy(iy); +// if (left < 0) { +// while (nhb < l) { +// nyb += 1; +// nhb += m_nhy(iy - nyb); +// } +// } +// if (right > 0) { +// while (nht < r) { +// nyt += 1; +// nht += m_nhy(iy + nyt); +// } +// } +// nyt++; + +// // Allocate output + +// size_t N = 0; + +// for (size_t i = iy - nyb; i < iy + nyt; ++i) { +// if (m_refine(iy) == -1) { +// N += w / m_nhx(i); // no refinement +// } +// else { +// N += 4 * w / m_nhx(i); // refinement +// } +// } + +// xt::xtensor ret = xt::empty({N}); + +// // Fill output + +// N = 0; +// size_t n; + +// std::cout << w << std::endl; +// std::cout << m_nhx << std::endl; + +// for (size_t i = iy - nyb; i < iy + nyt; ++i) { +// if (m_refine(iy) == -1) { +// n = w / m_nhx(i); // no refinement +// } +// else { +// n = 4 * w / m_nhx(i); // refinement +// } +// xt::view(ret, xt::range(N, N + n)) = m_startElem(i) + xt::arange(n); +// N += n; +// } + +// return ret; +// } + 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) { 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)); 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((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.elementgrid(); xt::xtensor elmat_fine = m_fine.elementgrid(); 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 elementgrid = m_regular.elementgrid(); // 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(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 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 e5a43ad..d674316 100644 --- a/test/basic/MeshQuad4.cpp +++ b/test/basic/MeshQuad4.cpp @@ -1,526 +1,550 @@ #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::elementgrid_ravel") + { + GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5); + xt::xtensor 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})))); + + xt::xtensor 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 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 - 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::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::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 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_)); } }