diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h index a9d8958..a5b5507 100644 --- a/include/GooseFEM/Mesh.h +++ b/include/GooseFEM/Mesh.h @@ -1,772 +1,1573 @@ /** Generic mesh operations. \file Mesh.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESH_H #define GOOSEFEM_MESH_H #include "config.h" namespace GooseFEM { namespace Mesh { /** Enumerator for element-types */ enum class ElementType { Unknown, ///< Unknown element-type Quad4, ///< Quadrilateral: 4-noded element in 2-d Hex8, ///< Hexahedron: 8-noded element in 3-d Tri3 ///< Triangle: 3-noded element in 2-d }; /** Extract the element type based on the connectivity. \param coor Nodal coordinates. \param conn Connectivity. \return ElementType(). */ inline ElementType defaultElementType( const xt::xtensor& coor, const xt::xtensor& conn); /** Base-class for regular meshes in 2-d. This class does not have a specific element-type in mind, it is used mostly internally to derive from such that common methods do not have to be reimplementation. */ class RegularBase2d { public: RegularBase2d() = default; /** \return Number of elements. */ size_t nelem() const; /** \return Number of nodes. */ size_t nnode() const; /** \return Number of nodes-per-element == 4. */ size_t nne() const; /** \return Number of dimensions == 2. */ size_t ndim() const; /** \return Number of elements in x-direction == width of the mesh in units of #h. */ virtual size_t nelx() const; /** \return Number of elements in y-direction == height of the mesh, in units of #h, */ virtual size_t nely() const; /** \return Linear edge size of one 'block'. */ double h() const; /** \return The ElementType(). */ virtual ElementType getElementType() const; /** \return Nodal coordinates [#nnode, #ndim]. */ virtual xt::xtensor coor() const; /** \return Connectivity [#nelem, #nne]. */ virtual xt::xtensor conn() const; /** \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. */ virtual xt::xtensor dofs() const; /** Nodes along the bottom edge (y = 0), in order of increasing x. \return List of node numbers. */ virtual xt::xtensor nodesBottomEdge() const; /** Nodes along the top edge (y = #nely * #h), in order of increasing x. \return List of node numbers. */ virtual xt::xtensor nodesTopEdge() const; /** Nodes along the left edge (x = 0), in order of increasing y. \return List of node numbers. */ virtual xt::xtensor nodesLeftEdge() const; /** Nodes along the right edge (x = #nelx * #h), in order of increasing y. \return List of node numbers. */ virtual xt::xtensor nodesRightEdge() const; /** Nodes along the bottom edge (y = 0), without the corners (at x = 0 and x = #nelx * #h). Same as: nodesBottomEdge()[1: -1]. \return List of node numbers. */ virtual xt::xtensor nodesBottomOpenEdge() const; /** Nodes along the top edge (y = #nely * #h), without the corners (at x = 0 and x = #nelx * #h). Same as: nodesTopEdge()[1: -1]. \return List of node numbers. */ virtual xt::xtensor nodesTopOpenEdge() const; /** Nodes along the left edge (x = 0), without the corners (at y = 0 and y = #nely * #h). Same as: nodesLeftEdge()[1: -1]. \return List of node numbers. */ virtual xt::xtensor nodesLeftOpenEdge() const; /** Nodes along the right edge (x = #nelx * #h), without the corners (at y = 0 and y = #nely * #h). Same as: nodesRightEdge()[1: -1]. \return List of node numbers. */ virtual xt::xtensor nodesRightOpenEdge() const; /** The bottom-left corner node (at x = 0, y = 0). Same as nodesBottomEdge()[0] and nodesLeftEdge()[0]. \return Node number. */ virtual size_t nodesBottomLeftCorner() const; /** The bottom-right corner node (at x = #nelx * #h, y = 0). Same as nodesBottomEdge()[-1] and nodesRightEdge()[0]. \return Node number. */ virtual size_t nodesBottomRightCorner() const; /** The top-left corner node (at x = 0, y = #nely * #h). Same as nodesTopEdge()[0] and nodesRightEdge()[-1]. \return Node number. */ virtual size_t nodesTopLeftCorner() const; /** The top-right corner node (at x = #nelx * #h, y = #nely * #h). Same as nodesTopEdge()[-1] and nodesRightEdge()[-1]. \return Node number. */ virtual size_t nodesTopRightCorner() const; /** \return Alias of nodesBottomLeftCorner(). */ size_t nodesLeftBottomCorner() const; /** \return Alias of nodesTopLeftCorner(). */ size_t nodesLeftTopCorner() const; /** \return Alias of nodesBottomRightCorner(). */ size_t nodesRightBottomCorner() const; /** \return Alias of nodesTopRightCorner(). */ size_t nodesRightTopCorner() const; /** DOF-numbers for the case that the periodicity if fully eliminated. Such that: dofs[nodesRightOpenEdge(), :] = dofs[nodesLeftOpenEdge(), :] dofs[nodesTopOpenEdge(), :] = dofs[nodesBottomOpenEdge(), :] dofs[nodesBottomRightCorner(), :] = dofs[nodesBottomLeftCorner(), :] dofs[nodesTopRightCorner(), :] = dofs[nodesBottomLeftCorner(), :] dofs[nodesTopLeftCorner(), :] = dofs[nodesBottomLeftCorner(), :] \return DOF numbers for each node [#nnode, #ndim]. */ xt::xtensor dofsPeriodic() const; /** - Periodic node pairs, in two columns: (independent, dependent) + Periodic node pairs, in two columns: (independent, dependent). + + - nodesRightOpenEdge() are tied to nodesLeftOpenEdge(). + - nodesTopOpenEdge() are tied to nodesBottomOpenEdge(). + - nodesBottomRightCorner() are tied to nodesBottomLeftCorner(). + - nodesTopRightCorner() are tied to nodesBottomLeftCorner(). + - nodesTopLeftCorner() are tied to nodesBottomLeftCorner(). \return [ntyings, #ndim]. */ xt::xtensor nodesPeriodic() const; /** Reference node to use for periodicity, because all corners are tied to it. \return Alias of nodesBottomLeftCorner(). */ size_t nodesOrigin() const; protected: double m_h; ///< See h() size_t m_nelx; ///< See nelx() size_t m_nely; ///< See nely() size_t m_nelem; ///< See nelem() size_t m_nnode; ///< See nnode() size_t m_nne; ///< See nne() size_t m_ndim; ///< See ndim() }; +/** +Base-class for regular meshes in 3-d. +This class does not have a specific element-type in mind, it is used mostly internally +to derive from such that common methods do not have to be reimplementation. +*/ +class RegularBase3d { +public: + + RegularBase3d() = default; + + /** + \return Number of elements. + */ + size_t nelem() const; + + /** + \return Number of nodes. + */ + size_t nnode() const; + + /** + \return Number of nodes-per-element == 4. + */ + size_t nne() const; + + /** + \return Number of dimensions == 2. + */ + size_t ndim() const; + + /** + \return Number of elements in x-direction == width of the mesh in units of #h. + */ + virtual size_t nelx() const; + + /** + \return Number of elements in y-direction == height of the mesh, in units of #h, + */ + virtual size_t nely() const; + + /** + \return Number of elements in y-direction == height of the mesh, in units of #h, + */ + virtual size_t nelz() const; + + /** + \return Linear edge size of one 'block'. + */ + double h() const; + + /** + \return The ElementType(). + */ + virtual ElementType getElementType() const; + + /** + \return Nodal coordinates [#nnode, #ndim]. + */ + virtual xt::xtensor coor() const; + + /** + \return Connectivity [#nelem, #nne]. + */ + virtual xt::xtensor conn() const; + + /** + \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. + */ + virtual xt::xtensor dofs() const; + + /** + Nodes along the bottom face (y = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottom() const; + + /** + Nodes along the top face (y = #nely * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesTop() const; + + /** + Nodes along the left face (x = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesLeft() const; + + /** + Nodes along the right face (x = #nelx * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesRight() const; + + /** + Nodes along the front face (z = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesFront() const; + + /** + Nodes along the back face (z = #nelz * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBack() const; + + /** + Nodes along the edge at the intersection of the front and bottom faces + (z = 0 and y = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontBottomEdge() const; + + /** + Nodes along the edge at the intersection of the front and top faces + (z = 0 and y = #nely * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontTopEdge() const; + + /** + Nodes along the edge at the intersection of the front and left faces + (z = 0 and x = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontLeftEdge() const; + + /** + Nodes along the edge at the intersection of the front and right faces + (z = 0 and x = #nelx * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontRightEdge() const; + + /** + Nodes along the edge at the intersection of the back and bottom faces + (z = #nelz * #h and y = #nely * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackBottomEdge() const; + + /** + Nodes along the edge at the intersection of the back and top faces + (z = #nelz * #h and x = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackTopEdge() const; + + /** + Nodes along the edge at the intersection of the back and left faces + (z = #nelz * #h and x = #nelx * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackLeftEdge() const; + + /** + Nodes along the edge at the intersection of the back and right faces + (? = #nelz * #h and ? = ?). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackRightEdge() const; + + /** + Nodes along the edge at the intersection of the bottom and left faces + (y = 0 and x = 0). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottomLeftEdge() const; + + /** + Nodes along the edge at the intersection of the bottom and right faces + (y = 0 and x = #nelx * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottomRightEdge() const; + + /** + Nodes along the edge at the intersection of the top and left faces + (y = 0 and x = #nelx * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesTopLeftEdge() const; + + /** + Nodes along the edge at the intersection of the top and right faces + (y = #nely * #h and x = #nelx * #h). + + \return List of node numbers. + */ + virtual xt::xtensor nodesTopRightEdge() const; + + /** + \return Alias of nodesFrontBottomEdge() + */ + xt::xtensor nodesBottomFrontEdge() const; + + /** + \return Alias of nodesBackBottomEdge() + */ + xt::xtensor nodesBottomBackEdge() const; + + /** + \return Alias of nodesFrontTopEdge() + */ + xt::xtensor nodesTopFrontEdge() const; + + /** + \return Alias of nodesBackTopEdge() + */ + xt::xtensor nodesTopBackEdge() const; + + /** + \return Alias of nodesBottomLeftEdge() + */ + xt::xtensor nodesLeftBottomEdge() const; + + /** + \return Alias of nodesFrontLeftEdge() + */ + xt::xtensor nodesLeftFrontEdge() const; + + /** + \return Alias of nodesBackLeftEdge() + */ + xt::xtensor nodesLeftBackEdge() const; + + /** + \return Alias of nodesTopLeftEdge() + */ + xt::xtensor nodesLeftTopEdge() const; + + /** + \return Alias of nodesBottomRightEdge() + */ + xt::xtensor nodesRightBottomEdge() const; + + /** + \return Alias of nodesTopRightEdge() + */ + xt::xtensor nodesRightTopEdge() const; + + /** + \return Alias of nodesFrontRightEdge() + */ + xt::xtensor nodesRightFrontEdge() const; + + /** + \return Alias of nodesBackRightEdge() + */ + xt::xtensor nodesRightBackEdge() const; + + + /** + Nodes along the front face excluding edges. + Same as different between nodesFront() and + [nodesFrontBottomEdge(), nodesFrontTopEdge(), nodesFrontLeftEdge(), nodesFrontRightEdge()] + + \return list of node numbers. + */ + virtual xt::xtensor nodesFrontFace() const; + + /** + Nodes along the back face excluding edges. + Same as different between nodesBack() and + [nodesBackBottomEdge(), nodesBackTopEdge(), nodesBackLeftEdge(), nodesBackRightEdge()] + + \return list of node numbers. + */ + virtual xt::xtensor nodesBackFace() const; + + /** + Nodes along the left face excluding edges. + Same as different between nodesLeft() and + [nodesFrontLeftEdge(), nodesBackLeftEdge(), nodesBottomLeftEdge(), nodesTopLeftEdge()] + + \return list of node numbers. + */ + virtual xt::xtensor nodesLeftFace() const; + + /** + Nodes along the right face excluding edges. + Same as different between nodesRight() and + [nodesFrontRightEdge(), nodesBackRightEdge(), nodesBottomRightEdge(), nodesTopRightEdge()] + + \return list of node numbers. + */ + virtual xt::xtensor nodesRightFace() const; + + /** + Nodes along the bottom face excluding edges. + Same as different between nodesBottom() and + [nodesBackBottomEdge(), nodesBackTopEdge(), nodesBackLeftEdge(), nodesBackRightEdge()] + + \return list of node numbers. + */ + virtual xt::xtensor nodesBottomFace() const; + + /** + Nodes along the top face excluding edges. + Same as different between nodesTop() and + [nodesFrontBottomEdge(), nodesFrontTopEdge(), nodesFrontLeftEdge(), nodesFrontRightEdge()] + + \return list of node numbers. + */ + virtual xt::xtensor nodesTopFace() const; + + /** + Same as nodesFrontBottomEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontBottomOpenEdge() const; + + /** + Same as nodesFrontTopEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontTopOpenEdge() const; + + /** + Same as nodesFrontLeftEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontLeftOpenEdge() const; + + /** + Same as nodesFrontRightEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesFrontRightOpenEdge() const; + + /** + Same as nodesBackBottomEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackBottomOpenEdge() const; + + /** + Same as nodesBackTopEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackTopOpenEdge() const; + + /** + Same as nodesBackLeftEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackLeftOpenEdge() const; + + /** + Same as nodesBackRightEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBackRightOpenEdge() const; + + /** + Same as nodesBottomLeftEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottomLeftOpenEdge() const; + + /** + Same as nodesBottomRightEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottomRightOpenEdge() const; + + /** + Same as nodesTopLeftEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesTopLeftOpenEdge() const; + + /** + Same as nodesTopRightEdge() but without corners. + + \return List of node numbers. + */ + virtual xt::xtensor nodesTopRightOpenEdge() const; + + /** + \return Alias of nodesFrontBottomOpenEdge(). + */ + xt::xtensor nodesBottomFrontOpenEdge() const; + + /** + \return Alias of nodesBackBottomOpenEdge(). + */ + xt::xtensor nodesBottomBackOpenEdge() const; + + /** + \return Alias of nodesFrontTopOpenEdge(). + */ + xt::xtensor nodesTopFrontOpenEdge() const; + + /** + \return Alias of nodesBackTopOpenEdge(). + */ + xt::xtensor nodesTopBackOpenEdge() const; + + /** + \return Alias of nodesBottomLeftOpenEdge(). + */ + xt::xtensor nodesLeftBottomOpenEdge() const; + + /** + \return Alias of nodesFrontLeftOpenEdge(). + */ + xt::xtensor nodesLeftFrontOpenEdge() const; + + /** + \return Alias of nodesBackLeftOpenEdge(). + */ + xt::xtensor nodesLeftBackOpenEdge() const; + + /** + \return Alias of nodesTopLeftOpenEdge(). + */ + xt::xtensor nodesLeftTopOpenEdge() const; + + /** + \return Alias of nodesBottomRightOpenEdge(). + */ + xt::xtensor nodesRightBottomOpenEdge() const; + + /** + \return Alias of nodesTopRightOpenEdge(). + */ + xt::xtensor nodesRightTopOpenEdge() const; + + /** + \return Alias of nodesFrontRightOpenEdge(). + */ + xt::xtensor nodesRightFrontOpenEdge() const; + + /** + \return Alias of nodesBackRightOpenEdge(). + */ + xt::xtensor nodesRightBackOpenEdge() const; + + /** + Front-Bottom-Left corner node. + + \return Node number. + */ + virtual size_t nodesFrontBottomLeftCorner() const; + + /** + Front-Bottom-Right corner node. + + \return Node number. + */ + virtual size_t nodesFrontBottomRightCorner() const; + + /** + Front-Top-Left corner node. + + \return Node number. + */ + virtual size_t nodesFrontTopLeftCorner() const; + + /** + Front-Top-Right corner node. + + \return Node number. + */ + virtual size_t nodesFrontTopRightCorner() const; + + /** + Back-Bottom-Left corner node. + + \return Node number. + */ + virtual size_t nodesBackBottomLeftCorner() const; + + /** + Back-Bottom-Right corner node. + + \return Node number. + */ + virtual size_t nodesBackBottomRightCorner() const; + + /** + Back-Top-Left corner node. + + \return Node number. + */ + virtual size_t nodesBackTopLeftCorner() const; + + /** + Back-Top-Right corner node. + + \return Node number. + */ + virtual size_t nodesBackTopRightCorner() const; + + /** + \return Alias of nodesFrontBottomLeftCorner(). + */ + size_t nodesFrontLeftBottomCorner() const; + + /** + \return Alias of nodesFrontBottomLeftCorner(). + */ + size_t nodesBottomFrontLeftCorner() const; + + /** + \return Alias of nodesFrontBottomLeftCorner(). + */ + size_t nodesBottomLeftFrontCorner() const; + + /** + \return Alias of nodesFrontBottomLeftCorner(). + */ + size_t nodesLeftFrontBottomCorner() const; + + /** + \return Alias of nodesFrontBottomLeftCorner(). + */ + size_t nodesLeftBottomFrontCorner() const; + + /** + \return Alias of nodesFrontBottomRightCorner(). + */ + size_t nodesFrontRightBottomCorner() const; + + /** + \return Alias of nodesFrontBottomRightCorner(). + */ + size_t nodesBottomFrontRightCorner() const; + + /** + \return Alias of nodesFrontBottomRightCorner(). + */ + size_t nodesBottomRightFrontCorner() const; + + /** + \return Alias of nodesFrontBottomRightCorner(). + */ + size_t nodesRightFrontBottomCorner() const; + + /** + \return Alias of nodesFrontBottomRightCorner(). + */ + size_t nodesRightBottomFrontCorner() const; + + /** + \return Alias of nodesFrontTopLeftCorner(). + */ + size_t nodesFrontLeftTopCorner() const; + + /** + \return Alias of nodesFrontTopLeftCorner(). + */ + size_t nodesTopFrontLeftCorner() const; + + /** + \return Alias of nodesFrontTopLeftCorner(). + */ + size_t nodesTopLeftFrontCorner() const; + + /** + \return Alias of nodesFrontTopLeftCorner(). + */ + size_t nodesLeftFrontTopCorner() const; + + /** + \return Alias of nodesFrontTopLeftCorner(). + */ + size_t nodesLeftTopFrontCorner() const; + + /** + \return Alias of nodesFrontTopRightCorner(). + */ + size_t nodesFrontRightTopCorner() const; + + /** + \return Alias of nodesFrontTopRightCorner(). + */ + size_t nodesTopFrontRightCorner() const; + + /** + \return Alias of nodesFrontTopRightCorner(). + */ + size_t nodesTopRightFrontCorner() const; + + /** + \return Alias of nodesFrontTopRightCorner(). + */ + size_t nodesRightFrontTopCorner() const; + + /** + \return Alias of nodesFrontTopRightCorner(). + */ + size_t nodesRightTopFrontCorner() const; + + /** + \return Alias of nodesBackBottomLeftCorner(). + */ + size_t nodesBackLeftBottomCorner() const; + + /** + \return Alias of nodesBackBottomLeftCorner(). + */ + size_t nodesBottomBackLeftCorner() const; + + /** + \return Alias of nodesBackBottomLeftCorner(). + */ + size_t nodesBottomLeftBackCorner() const; + + /** + \return Alias of nodesBackBottomLeftCorner(). + */ + size_t nodesLeftBackBottomCorner() const; + + /** + \return Alias of nodesBackBottomLeftCorner(). + */ + size_t nodesLeftBottomBackCorner() const; + + /** + \return Alias of nodesBackBottomRightCorner(). + */ + size_t nodesBackRightBottomCorner() const; + + /** + \return Alias of nodesBackBottomRightCorner(). + */ + size_t nodesBottomBackRightCorner() const; + + /** + \return Alias of nodesBackBottomRightCorner(). + */ + size_t nodesBottomRightBackCorner() const; + + /** + \return Alias of nodesBackBottomRightCorner(). + */ + size_t nodesRightBackBottomCorner() const; + + /** + \return Alias of nodesBackBottomRightCorner(). + */ + size_t nodesRightBottomBackCorner() const; + + /** + \return Alias of nodesBackTopLeftCorner(). + */ + size_t nodesBackLeftTopCorner() const; + + /** + \return Alias of nodesBackTopLeftCorner(). + */ + size_t nodesTopBackLeftCorner() const; + + /** + \return Alias of nodesBackTopLeftCorner(). + */ + size_t nodesTopLeftBackCorner() const; + + /** + \return Alias of nodesBackTopLeftCorner(). + */ + size_t nodesLeftBackTopCorner() const; + + /** + \return Alias of nodesBackTopLeftCorner(). + */ + size_t nodesLeftTopBackCorner() const; + + /** + \return Alias of nodesBackTopRightCorner(). + */ + size_t nodesBackRightTopCorner() const; + + /** + \return Alias of nodesBackTopRightCorner(). + */ + size_t nodesTopBackRightCorner() const; + + /** + \return Alias of nodesBackTopRightCorner(). + */ + size_t nodesTopRightBackCorner() const; + + /** + \return Alias of nodesBackTopRightCorner(). + */ + size_t nodesRightBackTopCorner() const; + + /** + \return Alias of nodesBackTopRightCorner(). + */ + size_t nodesRightTopBackCorner() const; + + /** + DOF-numbers for the case that the periodicity if fully eliminated. Such that: + + dofs[nodesBackFace(), :] = dofs[nodesFrontFace(), :] + dofs[nodesRightFace(), :] = dofs[nodesLeftFace(), :] + dofs[nodesTopFace(), :] = dofs[nodesBottomFace(), :] + dofs[nodesBackBottomOpenEdge(), :] = dofs[nodesFrontBottomOpenEdge(), :] + dofs[nodesBackTopOpenEdge(), :] = dofs[nodesFrontBottomOpenEdge(), :] + dofs[nodesFrontTopOpenEdge(), :] = dofs[nodesFrontBottomOpenEdge(), :] + dofs[nodesBottomRightOpenEdge(), :] = dofs[nodesBottomLeftOpenEdge(), :] + dofs[nodesTopRightOpenEdge(), :] = dofs[nodesBottomLeftOpenEdge(), :] + dofs[nodesTopLeftOpenEdge(), :] = dofs[nodesBottomLeftOpenEdge(), :] + dofs[nodesFrontRightOpenEdge(), :] = dofs[nodesFrontLeftOpenEdge(), :] + dofs[nodesBackRightOpenEdge(), :] = dofs[nodesFrontLeftOpenEdge(), :] + dofs[nodesBackLeftOpenEdge(), :] = dofs[nodesFrontLeftOpenEdge(), :] + dofs[nodesFrontBottomRightCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + dofs[nodesBackBottomRightCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + dofs[nodesBackBottomLeftCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + dofs[nodesFrontTopLeftCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + dofs[nodesFrontTopRightCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + dofs[nodesBackTopRightCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + dofs[nodesBackTopLeftCorner(), :] = dofs[nodesFrontBottomLeftCorner(), :] + + \return DOF numbers for each node [#nnode, #ndim]. + */ + xt::xtensor dofsPeriodic() const; + + /** + Periodic node pairs, in two columns: (independent, dependent). + + - nodesBackFace() are tied to nodesFrontFace(). + - nodesRightFace() are tied to nodesLeftFace(). + - nodesTopFace() are tied to nodesBottomFace(). + - nodesBackBottomOpenEdge() are tied to nodesFrontBottomOpenEdge(). + - nodesBackTopOpenEdge() are tied to nodesFrontBottomOpenEdge(). + - nodesFrontTopOpenEdge() are tied to nodesFrontBottomOpenEdge(). + - nodesBottomRightOpenEdge() are tied to nodesBottomLeftOpenEdge(). + - nodesTopRightOpenEdge() are tied to nodesBottomLeftOpenEdge(). + - nodesTopLeftOpenEdge() are tied to nodesBottomLeftOpenEdge(). + - nodesFrontRightOpenEdge() are tied to nodesFrontLeftOpenEdge(). + - nodesBackRightOpenEdge() are tied to nodesFrontLeftOpenEdge(). + - nodesBackLeftOpenEdge() are tied to nodesFrontLeftOpenEdge(). + - nodesFrontBottomRightCorner() are tied to nodesFrontBottomLeftCorner(). + - nodesBackBottomRightCorner() are tied to nodesFrontBottomLeftCorner(). + - nodesBackBottomLeftCorner() are tied to nodesFrontBottomLeftCorner(). + - nodesFrontTopLeftCorner() are tied to nodesFrontBottomLeftCorner(). + - nodesFrontTopRightCorner() are tied to nodesFrontBottomLeftCorner(). + - nodesBackTopRightCorner() are tied to nodesFrontBottomLeftCorner(). + - nodesBackTopLeftCorner() are tied to nodesFrontBottomLeftCorner(). + + \return [ntyings, #ndim]. + */ + xt::xtensor nodesPeriodic() const; + + /** + Reference node to use for periodicity, because all corners are tied to it. + + \return Alias of nodesFrontBottomLeftCorner(). + */ + size_t nodesOrigin() const; + +protected: + double m_h; ///< See h() + size_t m_nelx; ///< See nelx() + size_t m_nely; ///< See nely() + size_t m_nelz; ///< See nely() + size_t m_nelem; ///< See nelem() + size_t m_nnode; ///< See nnode() + size_t m_nne; ///< See nne() + size_t m_ndim; ///< See ndim() +}; + /** Find overlapping nodes. The output has the following structure: [[nodes_from_mesh_a], [nodes_from_mesh_b]] \param coor_a Nodal coordinates of mesh "a". \param coor_b Nodal coordinates of mesh "b". \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. \return Overlapping nodes. */ inline xt::xtensor overlapping( const xt::xtensor& coor_a, const xt::xtensor& coor_b, double rtol = 1e-5, double atol = 1e-8); /** Stitch two mesh objects, specifying overlapping nodes by hand. */ class ManualStitch { public: ManualStitch() = default; /** \param coor_a Nodal coordinates of mesh "a". \param conn_a Connectivity of mesh "a". \param overlapping_nodes_a Node-numbers of mesh "a" that overlap with mesh "b". \param coor_b Nodal coordinates of mesh "b". \param conn_b Connectivity of mesh "b". \param overlapping_nodes_b Node-numbers of mesh "b" that overlap with mesh "a". \param check_position If ``true`` the nodes are checked for position overlap. \param rtol Relative tolerance for check on position overlap. \param atol Absolute tolerance for check on position overlap. */ ManualStitch( const xt::xtensor& coor_a, const xt::xtensor& conn_a, const xt::xtensor& overlapping_nodes_a, const xt::xtensor& coor_b, const xt::xtensor& conn_b, const xt::xtensor& overlapping_nodes_b, bool check_position = true, double rtol = 1e-5, double atol = 1e-8); /** \return Number of sub meshes == 2. */ size_t nmesh() const; /** \return Number of elements. */ size_t nelem() const; /** \return Number of nodes. */ size_t nnode() const; /** \return Number of nodes-per-element. */ size_t nne() const; /** \return Number of dimensions. */ size_t ndim() const; /** \return Nodal coordinates [#nnode, #ndim]. */ xt::xtensor coor() const; /** \return Connectivity [#nelem, #nne]. */ xt::xtensor conn() const; /** \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. */ xt::xtensor dofs() const; /** \return Node-map per sub-mesh. */ std::vector> nodemap() const; /** \return Element-map per sub-mesh. */ std::vector> elemmap() const; /** \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return Node-map for a given mesh. */ xt::xtensor nodemap(size_t mesh_index) const; /** \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return Element-map for a given mesh. */ xt::xtensor elemmap(size_t mesh_index) const; /** Convert set of node numbers for an original mesh to the stitched mesh. \param set List of node numbers. \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return List of node numbers for the stitched mesh. */ xt::xtensor nodeset(const xt::xtensor& set, size_t mesh_index) const; /** Convert set of element numbers for an original mesh to the stitched mesh. \param set List of element numbers. \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return List of element numbers for the stitched mesh. */ xt::xtensor elemset(const xt::xtensor& set, size_t mesh_index) const; private: xt::xtensor m_coor; xt::xtensor m_conn; xt::xtensor m_map_b; size_t m_nnd_a; size_t m_nel_a; size_t m_nel_b; }; /** Stitch mesh objects, automatically searching for overlapping nodes. */ class Stitch { public: /** \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. */ Stitch(double rtol = 1e-5, double atol = 1e-8); /** Add mesh to be stitched. \param coor Nodal coordinates. \param conn Connectivity. */ void push_back(const xt::xtensor& coor, const xt::xtensor& conn); /** \return Number of sub meshes. */ size_t nmesh() const; /** \return Number of elements. */ size_t nelem() const; /** \return Number of nodes. */ size_t nnode() const; /** \return Number of nodes-per-element. */ size_t nne() const; /** \return Number of dimensions. */ size_t ndim() const; /** \return Nodal coordinates [#nnode, #ndim]. */ xt::xtensor coor() const; /** \return Connectivity [#nelem, #nne]. */ xt::xtensor conn() const; /** \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. */ xt::xtensor dofs() const; /** \return Node-map per sub-mesh. */ std::vector> nodemap() const; /** \return Element-map per sub-mesh. */ std::vector> elemmap() const; /** The node numbers in the stitched mesh that are coming from a specific sub-mesh. \param mesh_index Index of the sub-mesh. \return List of node numbers. */ xt::xtensor nodemap(size_t mesh_index) const; /** The element numbers in the stitched mesh that are coming from a specific sub-mesh. \param mesh_index Index of the sub-mesh. \return List of element numbers. */ xt::xtensor elemmap(size_t mesh_index) const; /** Convert set of node-numbers for a sub-mesh to the stitched mesh. \param set List of node numbers. \param mesh_index Index of the sub-mesh. \return List of node numbers for the stitched mesh. */ xt::xtensor nodeset(const xt::xtensor& set, size_t mesh_index) const; /** Convert set of element-numbers for a sub-mesh to the stitched mesh. \param set List of element numbers. \param mesh_index Index of the sub-mesh. \return List of element numbers for the stitched mesh. */ xt::xtensor elemset(const xt::xtensor& set, size_t mesh_index) const; /** Combine set of node numbers for an original to the final mesh (removes duplicates). \param set List of node numbers per mesh. \return List of node numbers for the stitched mesh. */ xt::xtensor nodeset(const std::vector>& set) const; /** Combine set of element numbers for an original to the final mesh. \param set List of element numbers per mesh. \return List of element numbers for the stitched mesh. */ xt::xtensor elemset(const std::vector>& set) const; private: xt::xtensor m_coor; xt::xtensor m_conn; std::vector> m_map; std::vector m_nel; ///< Number of elements per sub-mesh. std::vector m_el_offset; double m_rtol; double m_atol; }; /** \rst Renumber indices to lowest possible index. For example: .. math:: \begin{bmatrix} 0 & 1 \\ 5 & 4 \end{bmatrix} is renumbered to .. math:: \begin{bmatrix} 0 & 1 \\ 3 & 2 \end{bmatrix} Or, in pseudo-code, the result of this function is that: .. code-block:: python dofs = renumber(dofs) sort(unique(dofs[:])) == range(max(dofs+1)) .. tip:: One can use the wrapper function :cpp:func:`GooseFEM::Mesh::renumber`. This class gives more advanced features. \endrst */ class Renumber { public: Renumber() = default; /** \param dofs DOF-numbers. */ template Renumber(const T& dofs); /** Get renumbered DOFs (same as ``Renumber::apply(dofs)``). \param dofs List of (DOF-)numbers. \return Renumbered list of (DOF-)numbers. */ [[deprecated]] xt::xtensor get(const xt::xtensor& dofs) const; /** Apply renumbering to other set. \param list List of (DOF-)numbers. \return Renumbered list of (DOF-)numbers. */ template T apply(const T& list) const; /** Get the list needed to renumber, e.g.: dofs_renumbered(i, j) = index(dofs(i, j)) \return Renumber-index. */ xt::xtensor index() const; private: xt::xtensor m_renum; }; /** Renumber to lowest possible index (see GooseFEM::Mesh::Renumber). \param dofs DOF-numbers. \return Renumbered DOF-numbers. */ inline xt::xtensor renumber(const xt::xtensor& dofs); /** Reorder to lowest possible index, in specific order. For example for ``Reorder({iiu, iip})`` after reordering: iiu = xt::range(nnu); iip = xt::range(nnp) + nnu; */ class Reorder { public: Reorder() = default; /** \param args List of (DOF-)numbers. */ Reorder(const std::initializer_list> args); /** Get reordered DOFs (same as ``Reorder::apply(dofs)``). \param dofs List of (DOF-)numbers. \return Reordered list of (DOF-)numbers. */ [[deprecated]] xt::xtensor get(const xt::xtensor& dofs) const; /** Apply reordering to other set. \param list List of (DOF-)numbers. \return Reordered list of (DOF-)numbers. */ template T apply(const T& list) const; /** Get the list needed to reorder, e.g.: dofs_reordered(i, j) = index(dofs(i, j)) \return Reorder-index. */ xt::xtensor index() const; private: xt::xtensor m_renum; }; /** List with DOF-numbers in sequential order. The output is a sequential list of DOF-numbers for each vector-component of each node. For example for 3 nodes in 2 dimensions the output is \rst .. math:: \begin{bmatrix} 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \endrst \param nnode Number of nodes. \param ndim Number of dimensions. \return DOF-numbers. */ inline xt::xtensor dofs(size_t nnode, size_t ndim); /** Number of elements connected to each node. \param conn Connectivity. \return Coordination per node. */ inline xt::xtensor coordination(const xt::xtensor& conn); /** Elements connected to each node. \param conn Connectivity. \param sorted If ``true`` the output is sorted. \return Elements per node. */ inline std::vector> elem2node( const xt::xtensor& conn, bool sorted = true); /** Return size of each element edge. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Edge-sizes per element. */ inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type); /** Return size of each element edge. The element-type is automatically determined, see defaultElementType(). \param coor Nodal coordinates. \param conn Connectivity. \return Edge-sizes per element. */ inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn); /** Coordinates of the center of each element. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Center of each element. */ inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type); /** Coordinates of the center of each element. The element-type is automatically determined, see defaultElementType(). \param coor Nodal coordinates. \param conn Connectivity. \return Center of each element. */ inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn); /** Convert an element-map to a node-map. \param elem_map Element-map such that ``new_elvar = elvar[elem_map]``. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Node-map such that ``new_nodevar = nodevar[node_map]`` */ inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn, ElementType type); /** Convert an element-map to a node-map. The element-type is automatically determined, see defaultElementType(). \param elem_map Element-map such that ``new_elvar = elvar[elem_map]``. \param coor Nodal coordinates. \param conn Connectivity. \return Node-map such that ``new_nodevar = nodevar[node_map]`` */ inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn); } // namespace Mesh } // namespace GooseFEM #include "Mesh.hpp" #endif diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp index 9259e51..f211a7f 100644 --- a/include/GooseFEM/Mesh.hpp +++ b/include/GooseFEM/Mesh.hpp @@ -1,841 +1,1554 @@ /** Implementation of Mesh.h \file Mesh.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESH_HPP #define GOOSEFEM_MESH_HPP #include "Mesh.h" namespace GooseFEM { namespace Mesh { inline ElementType defaultElementType( const xt::xtensor& coor, const xt::xtensor& conn) { if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) { return ElementType::Tri3; } if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) { return ElementType::Quad4; } if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) { return ElementType::Hex8; } throw std::runtime_error("Element-type not implemented"); } inline size_t RegularBase2d::nelem() const { return m_nelem; } inline size_t RegularBase2d::nnode() const { return m_nnode; } inline size_t RegularBase2d::nne() const { return m_nne; } inline size_t RegularBase2d::ndim() const { return m_ndim; } inline size_t RegularBase2d::nelx() const { return m_nelx; } inline size_t RegularBase2d::nely() const { return m_nely; } inline double RegularBase2d::h() const { return m_h; } inline ElementType RegularBase2d::getElementType() const { return ElementType::Unknown; } inline xt::xtensor RegularBase2d::coor() const { return xt::empty({0, 0}); } inline xt::xtensor RegularBase2d::conn() const { return xt::empty({0, 0}); } inline xt::xtensor RegularBase2d::nodesBottomEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesTopEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesLeftEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesRightEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesBottomOpenEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesTopOpenEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesLeftOpenEdge() const { return xt::empty({0}); } inline xt::xtensor RegularBase2d::nodesRightOpenEdge() const { return xt::empty({0}); } inline size_t RegularBase2d::nodesBottomLeftCorner() const { return 0; } inline size_t RegularBase2d::nodesBottomRightCorner() const { return 0; } inline size_t RegularBase2d::nodesTopLeftCorner() const { return 0; } inline size_t RegularBase2d::nodesTopRightCorner() const { return 0; } inline size_t RegularBase2d::nodesLeftBottomCorner() const { - return nodesBottomLeftCorner(); + return this->nodesBottomLeftCorner(); } inline size_t RegularBase2d::nodesLeftTopCorner() const { - return nodesTopLeftCorner(); + return this->nodesTopLeftCorner(); } inline size_t RegularBase2d::nodesRightBottomCorner() const { - return nodesBottomRightCorner(); + return this->nodesBottomRightCorner(); } inline size_t RegularBase2d::nodesRightTopCorner() const { - return nodesTopRightCorner(); + return this->nodesTopRightCorner(); } inline xt::xtensor RegularBase2d::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() + size_t(3), size_t(2)}; 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 RegularBase2d::nodesOrigin() const { return nodesBottomLeftCorner(); } inline xt::xtensor RegularBase2d::dofs() const { return GooseFEM::Mesh::dofs(this->nnode(), this->ndim()); } inline xt::xtensor RegularBase2d::dofsPeriodic() const { xt::xtensor ret = this->dofs(); xt::xtensor nodePer = this->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 < this->ndim(); ++j) { xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j); } return GooseFEM::Mesh::renumber(ret); } +inline size_t RegularBase3d::nelem() const +{ + return m_nelem; +} + +inline size_t RegularBase3d::nnode() const +{ + return m_nnode; +} + +inline size_t RegularBase3d::nne() const +{ + return m_nne; +} + +inline size_t RegularBase3d::ndim() const +{ + return m_ndim; +} + +inline size_t RegularBase3d::nelx() const +{ + return m_nelx; +} + +inline size_t RegularBase3d::nely() const +{ + return m_nely; +} + +inline size_t RegularBase3d::nelz() const +{ + return m_nelz; +} + +inline double RegularBase3d::h() const +{ + return m_h; +} + +inline ElementType RegularBase3d::getElementType() const +{ + return ElementType::Unknown; +} + +inline xt::xtensor RegularBase3d::coor() const +{ + return xt::empty({0, 0}); +} + +inline xt::xtensor RegularBase3d::conn() const +{ + return xt::empty({0, 0}); +} + +inline xt::xtensor RegularBase3d::nodesBottom() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesTop() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesLeft() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesRight() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFront() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBack() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontBottomEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontTopEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontLeftEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontRightEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackBottomEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackTopEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackLeftEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackRightEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomLeftEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomRightEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesTopLeftEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesTopRightEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomFrontEdge() const +{ + return this->nodesFrontBottomEdge(); +} +inline xt::xtensor RegularBase3d::nodesBottomBackEdge() const +{ + return this->nodesBackBottomEdge(); +} +inline xt::xtensor RegularBase3d::nodesTopFrontEdge() const +{ + return this->nodesFrontTopEdge(); +} +inline xt::xtensor RegularBase3d::nodesTopBackEdge() const +{ + return this->nodesBackTopEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftBottomEdge() const +{ + return this->nodesBottomLeftEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftFrontEdge() const +{ + return this->nodesFrontLeftEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftBackEdge() const +{ + return this->nodesBackLeftEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftTopEdge() const +{ + return this->nodesTopLeftEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightBottomEdge() const +{ + return this->nodesBottomRightEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightTopEdge() const +{ + return this->nodesTopRightEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightFrontEdge() const +{ + return this->nodesFrontRightEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightBackEdge() const +{ + return this->nodesBackRightEdge(); +} + +inline xt::xtensor RegularBase3d::nodesFrontFace() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackFace() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesLeftFace() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesRightFace() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomFace() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesTopFace() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontBottomOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontTopOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontLeftOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesFrontRightOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackBottomOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackTopOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackLeftOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBackRightOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomLeftOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomRightOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesTopLeftOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesTopRightOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase3d::nodesBottomFrontOpenEdge() const +{ + return this->nodesFrontBottomOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesBottomBackOpenEdge() const +{ + return this->nodesBackBottomOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesTopFrontOpenEdge() const +{ + return this->nodesFrontTopOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesTopBackOpenEdge() const +{ + return this->nodesBackTopOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftBottomOpenEdge() const +{ + return this->nodesBottomLeftOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftFrontOpenEdge() const +{ + return this->nodesFrontLeftOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftBackOpenEdge() const +{ + return this->nodesBackLeftOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesLeftTopOpenEdge() const +{ + return this->nodesTopLeftOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightBottomOpenEdge() const +{ + return this->nodesBottomRightOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightTopOpenEdge() const +{ + return this->nodesTopRightOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightFrontOpenEdge() const +{ + return this->nodesFrontRightOpenEdge(); +} +inline xt::xtensor RegularBase3d::nodesRightBackOpenEdge() const +{ + return this->nodesBackRightOpenEdge(); +} + +inline size_t RegularBase3d::nodesFrontBottomLeftCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesFrontBottomRightCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesFrontTopLeftCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesFrontTopRightCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesBackBottomLeftCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesBackBottomRightCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesBackTopLeftCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesBackTopRightCorner() const +{ + return 0; +} + +inline size_t RegularBase3d::nodesFrontLeftBottomCorner() const +{ + return this->nodesFrontBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesBottomFrontLeftCorner() const +{ + return this->nodesFrontBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesBottomLeftFrontCorner() const +{ + return this->nodesFrontBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftFrontBottomCorner() const +{ + return this->nodesFrontBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftBottomFrontCorner() const +{ + return this->nodesFrontBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesFrontRightBottomCorner() const +{ + return this->nodesFrontBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesBottomFrontRightCorner() const +{ + return this->nodesFrontBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesBottomRightFrontCorner() const +{ + return this->nodesFrontBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesRightFrontBottomCorner() const +{ + return this->nodesFrontBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesRightBottomFrontCorner() const +{ + return this->nodesFrontBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesFrontLeftTopCorner() const +{ + return this->nodesFrontTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesTopFrontLeftCorner() const +{ + return this->nodesFrontTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesTopLeftFrontCorner() const +{ + return this->nodesFrontTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftFrontTopCorner() const +{ + return this->nodesFrontTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftTopFrontCorner() const +{ + return this->nodesFrontTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesFrontRightTopCorner() const +{ + return this->nodesFrontTopRightCorner(); +} + +inline size_t RegularBase3d::nodesTopFrontRightCorner() const +{ + return this->nodesFrontTopRightCorner(); +} + +inline size_t RegularBase3d::nodesTopRightFrontCorner() const +{ + return this->nodesFrontTopRightCorner(); +} + +inline size_t RegularBase3d::nodesRightFrontTopCorner() const +{ + return this->nodesFrontTopRightCorner(); +} + +inline size_t RegularBase3d::nodesRightTopFrontCorner() const +{ + return this->nodesFrontTopRightCorner(); +} + +inline size_t RegularBase3d::nodesBackLeftBottomCorner() const +{ + return this->nodesBackBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesBottomBackLeftCorner() const +{ + return this->nodesBackBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesBottomLeftBackCorner() const +{ + return this->nodesBackBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftBackBottomCorner() const +{ + return this->nodesBackBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftBottomBackCorner() const +{ + return this->nodesBackBottomLeftCorner(); +} + +inline size_t RegularBase3d::nodesBackRightBottomCorner() const +{ + return this->nodesBackBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesBottomBackRightCorner() const +{ + return this->nodesBackBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesBottomRightBackCorner() const +{ + return this->nodesBackBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesRightBackBottomCorner() const +{ + return this->nodesBackBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesRightBottomBackCorner() const +{ + return this->nodesBackBottomRightCorner(); +} + +inline size_t RegularBase3d::nodesBackLeftTopCorner() const +{ + return this->nodesBackTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesTopBackLeftCorner() const +{ + return this->nodesBackTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesTopLeftBackCorner() const +{ + return this->nodesBackTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftBackTopCorner() const +{ + return this->nodesBackTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesLeftTopBackCorner() const +{ + return this->nodesBackTopLeftCorner(); +} + +inline size_t RegularBase3d::nodesBackRightTopCorner() const +{ + return this->nodesBackTopRightCorner(); +} + +inline size_t RegularBase3d::nodesTopBackRightCorner() const +{ + return this->nodesBackTopRightCorner(); +} + +inline size_t RegularBase3d::nodesTopRightBackCorner() const +{ + return this->nodesBackTopRightCorner(); +} + +inline size_t RegularBase3d::nodesRightBackTopCorner() const +{ + return this->nodesBackTopRightCorner(); +} + +inline size_t RegularBase3d::nodesRightTopBackCorner() const +{ + return this->nodesBackTopRightCorner(); +} + +inline xt::xtensor RegularBase3d::nodesPeriodic() const +{ + xt::xtensor fro = nodesFrontFace(); + xt::xtensor bck = nodesBackFace(); + xt::xtensor lft = nodesLeftFace(); + xt::xtensor rgt = nodesRightFace(); + xt::xtensor bot = nodesBottomFace(); + xt::xtensor top = nodesTopFace(); + + xt::xtensor froBot = nodesFrontBottomOpenEdge(); + xt::xtensor froTop = nodesFrontTopOpenEdge(); + xt::xtensor froLft = nodesFrontLeftOpenEdge(); + xt::xtensor froRgt = nodesFrontRightOpenEdge(); + xt::xtensor bckBot = nodesBackBottomOpenEdge(); + xt::xtensor bckTop = nodesBackTopOpenEdge(); + xt::xtensor bckLft = nodesBackLeftOpenEdge(); + xt::xtensor bckRgt = nodesBackRightOpenEdge(); + xt::xtensor botLft = nodesBottomLeftOpenEdge(); + xt::xtensor botRgt = nodesBottomRightOpenEdge(); + xt::xtensor topLft = nodesTopLeftOpenEdge(); + xt::xtensor topRgt = nodesTopRightOpenEdge(); + + size_t tface = fro.size() + lft.size() + bot.size(); + size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size(); + size_t tnode = 7; + xt::xtensor ret = xt::empty({tface + tedge + tnode, std::size_t(2)}); + + size_t i = 0; + + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesFrontBottomRightCorner(); + ++i; + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesBackBottomRightCorner(); + ++i; + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesBackBottomLeftCorner(); + ++i; + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesFrontTopLeftCorner(); + ++i; + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesFrontTopRightCorner(); + ++i; + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesBackTopRightCorner(); + ++i; + ret(i, 0) = nodesFrontBottomLeftCorner(); + ret(i, 1) = nodesBackTopLeftCorner(); + ++i; + + for (size_t j = 0; j < froBot.size(); ++j) { + ret(i, 0) = froBot(j); + ret(i, 1) = bckBot(j); + ++i; + } + for (size_t j = 0; j < froBot.size(); ++j) { + ret(i, 0) = froBot(j); + ret(i, 1) = bckTop(j); + ++i; + } + for (size_t j = 0; j < froBot.size(); ++j) { + ret(i, 0) = froBot(j); + ret(i, 1) = froTop(j); + ++i; + } + for (size_t j = 0; j < botLft.size(); ++j) { + ret(i, 0) = botLft(j); + ret(i, 1) = botRgt(j); + ++i; + } + for (size_t j = 0; j < botLft.size(); ++j) { + ret(i, 0) = botLft(j); + ret(i, 1) = topRgt(j); + ++i; + } + for (size_t j = 0; j < botLft.size(); ++j) { + ret(i, 0) = botLft(j); + ret(i, 1) = topLft(j); + ++i; + } + for (size_t j = 0; j < froLft.size(); ++j) { + ret(i, 0) = froLft(j); + ret(i, 1) = froRgt(j); + ++i; + } + for (size_t j = 0; j < froLft.size(); ++j) { + ret(i, 0) = froLft(j); + ret(i, 1) = bckRgt(j); + ++i; + } + for (size_t j = 0; j < froLft.size(); ++j) { + ret(i, 0) = froLft(j); + ret(i, 1) = bckLft(j); + ++i; + } + + for (size_t j = 0; j < fro.size(); ++j) { + ret(i, 0) = fro(j); + ret(i, 1) = bck(j); + ++i; + } + for (size_t j = 0; j < lft.size(); ++j) { + ret(i, 0) = lft(j); + ret(i, 1) = rgt(j); + ++i; + } + for (size_t j = 0; j < bot.size(); ++j) { + ret(i, 0) = bot(j); + ret(i, 1) = top(j); + ++i; + } + + return ret; +} + +inline size_t RegularBase3d::nodesOrigin() const +{ + return nodesFrontBottomLeftCorner(); +} + +inline xt::xtensor RegularBase3d::dofs() const +{ + return GooseFEM::Mesh::dofs(this->nnode(), this->ndim()); +} + +inline xt::xtensor RegularBase3d::dofsPeriodic() const +{ + xt::xtensor ret = this->dofs(); + xt::xtensor nodePer = this->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 < this->ndim(); ++j) { + xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j); + } + + return GooseFEM::Mesh::renumber(ret); +} + namespace detail { template inline T renum(const T& arg, const R& mapping) { T ret = T::from_shape(arg.shape()); auto jt = ret.begin(); for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) { *jt = mapping(*it); } return ret; } } // namespace detail inline xt::xtensor overlapping( const xt::xtensor& coor_a, const xt::xtensor& coor_b, double rtol, double atol) { GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); std::vector ret_a; std::vector ret_b; for (size_t i = 0; i < coor_a.shape(0); ++i) { auto idx = xt::flatten_indices(xt::argwhere(xt::prod(xt::isclose( coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1))); for (auto& j : idx) { ret_a.push_back(i); ret_b.push_back(j); } } xt::xtensor ret = xt::empty({size_t(2), ret_a.size()}); for (size_t i = 0; i < ret_a.size(); ++i) { ret(0, i) = ret_a[i]; ret(1, i) = ret_b[i]; } return ret; } inline ManualStitch::ManualStitch( const xt::xtensor& coor_a, const xt::xtensor& conn_a, const xt::xtensor& overlapping_nodes_a, const xt::xtensor& coor_b, const xt::xtensor& conn_b, const xt::xtensor& overlapping_nodes_b, bool check_position, double rtol, double atol) { UNUSED(rtol); UNUSED(atol); GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape())); GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1)); if (check_position) { GOOSEFEM_ASSERT(xt::allclose( xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()), xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()), rtol, atol)); } size_t nnda = coor_a.shape(0); size_t nndb = coor_b.shape(0); size_t ndim = coor_a.shape(1); size_t nelim = overlapping_nodes_a.size(); size_t nela = conn_a.shape(0); size_t nelb = conn_b.shape(0); size_t nne = conn_a.shape(1); m_nel_a = nela; m_nel_b = nelb; m_nnd_a = nnda; xt::xtensor keep_b = xt::setdiff1d(xt::arange(nndb), overlapping_nodes_b); m_map_b = xt::empty({nndb}); xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a; xt::view(m_map_b, xt::keep(keep_b)) = xt::arange(keep_b.size()) + nnda; m_conn = xt::empty({nela + nelb, nne}); xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a; xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b); m_coor = xt::empty({nnda + nndb - nelim, ndim}); xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a; xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) = xt::view(coor_b, xt::keep(keep_b), xt::all()); } inline xt::xtensor ManualStitch::coor() const { return m_coor; } inline xt::xtensor ManualStitch::conn() const { return m_conn; } inline size_t ManualStitch::nmesh() const { return 2; } inline size_t ManualStitch::nelem() const { return m_conn.shape(0); } inline size_t ManualStitch::nnode() const { return m_coor.shape(0); } inline size_t ManualStitch::nne() const { return m_conn.shape(1); } inline size_t ManualStitch::ndim() const { return m_coor.shape(1); } inline xt::xtensor ManualStitch::dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline std::vector> ManualStitch::nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } inline std::vector> ManualStitch::elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } inline xt::xtensor ManualStitch::nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nnd_a); } return m_map_b; } inline xt::xtensor ManualStitch::elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nel_a); } return xt::arange(m_nel_b) + m_nel_a; } inline xt::xtensor ManualStitch::nodeset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a); return set; } GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size()); return detail::renum(set, m_map_b); } inline xt::xtensor ManualStitch::elemset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a); return set; } GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b); return set + m_nel_a; } inline Stitch::Stitch(double rtol, double atol) { m_rtol = rtol; m_atol = atol; } inline void Stitch::push_back( const xt::xtensor& coor, const xt::xtensor& conn) { if (m_map.size() == 0) { m_coor = coor; m_conn = conn; m_map.push_back(xt::eval(xt::arange(coor.shape(0)))); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(0); return; } auto overlap = overlapping(m_coor, coor, m_rtol, m_atol); size_t index = m_map.size(); ManualStitch stich( m_coor, m_conn, xt::view(overlap, 0, xt::all()), coor, conn, xt::view(overlap, 1, xt::all()), false); m_coor = stich.coor(); m_conn = stich.conn(); m_map.push_back(stich.nodemap(1)); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]); } inline size_t Stitch::nmesh() const { return m_map.size(); } inline xt::xtensor Stitch::coor() const { return m_coor; } inline xt::xtensor Stitch::conn() const { return m_conn; } inline size_t Stitch::nelem() const { return m_conn.shape(0); } inline size_t Stitch::nnode() const { return m_coor.shape(0); } inline size_t Stitch::nne() const { return m_conn.shape(1); } inline size_t Stitch::ndim() const { return m_coor.shape(1); } inline xt::xtensor Stitch::dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline std::vector> Stitch::nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } inline std::vector> Stitch::elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } inline xt::xtensor Stitch::nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return m_map[mesh_index]; } inline xt::xtensor Stitch::elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return xt::arange(m_nel[mesh_index]) + m_el_offset[mesh_index]; } inline xt::xtensor Stitch::nodeset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size()); return detail::renum(set, m_map[mesh_index]); } inline xt::xtensor Stitch::elemset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]); return set + m_el_offset[mesh_index]; } inline xt::xtensor Stitch::nodeset(const std::vector>& set) const { GOOSEFEM_ASSERT(set.size() == m_map.size()); size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->nodeset(set[i], i); n += set[i].size(); } return xt::unique(ret); } inline xt::xtensor Stitch::elemset(const std::vector>& set) const { GOOSEFEM_ASSERT(set.size() == m_map.size()); size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i); n += set[i].size(); } return ret; } template inline Renumber::Renumber(const T& dofs) { size_t n = xt::amax(dofs)() + 1; size_t i = 0; xt::xtensor unique = xt::unique(dofs); m_renum = xt::empty({n}); for (auto& j : unique) { m_renum(j) = i; ++i; } } inline xt::xtensor Renumber::get(const xt::xtensor& dofs) const { GOOSEFEM_WARNING("Renumber::get is deprecated, use Renumber::apply"); return this->apply(dofs); } template inline T Renumber::apply(const T& list) const { return detail::renum(list, m_renum); } inline xt::xtensor Renumber::index() const { return m_renum; } inline xt::xtensor renumber(const xt::xtensor& dofs) { return Renumber(dofs).apply(dofs); } inline Reorder::Reorder(const std::initializer_list> args) { size_t n = 0; size_t i = 0; for (auto& arg : args) { if (arg.size() == 0) { continue; } n = std::max(n, xt::amax(arg)() + 1); } #ifdef GOOSEFEM_ENABLE_ASSERT for (auto& arg : args) { GOOSEFEM_ASSERT(xt::unique(arg) == xt::sort(arg)); } #endif m_renum = xt::empty({n}); for (auto& arg : args) { for (auto& j : arg) { m_renum(j) = i; ++i; } } } inline xt::xtensor Reorder::get(const xt::xtensor& dofs) const { GOOSEFEM_WARNING("Reorder::get is deprecated, use Reorder::apply"); return this->apply(dofs); } template inline T Reorder::apply(const T& list) const { T ret = T::from_shape(list.shape()); auto jt = ret.begin(); for (auto it = list.begin(); it != list.end(); ++it, ++jt) { *jt = m_renum(*it); } return ret; } inline xt::xtensor Reorder::index() const { return m_renum; } inline xt::xtensor dofs(size_t nnode, size_t ndim) { return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline xt::xtensor coordination(const xt::xtensor& conn) { size_t nnode = xt::amax(conn)() + 1; xt::xtensor N = xt::zeros({nnode}); for (auto it = conn.begin(); it != conn.end(); ++it) { N(*it) += 1; } return N; } inline std::vector> elem2node(const xt::xtensor& conn, bool sorted) { auto N = coordination(conn); auto nnode = N.size(); std::vector> ret(nnode); for (size_t i = 0; i < nnode; ++i) { ret[i].reserve(N(i)); } for (size_t e = 0; e < conn.shape(0); ++e) { for (size_t m = 0; m < conn.shape(1); ++m) { ret[conn(e, m)].push_back(e); } } if (sorted) { for (auto& row : ret) { std::sort(row.begin(), row.end()); } } return ret; } inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); if (type == ElementType::Quad4) { GOOSEFEM_ASSERT(coor.shape(1) == 2ul); GOOSEFEM_ASSERT(conn.shape(1) == 4ul); xt::xtensor n0 = xt::view(conn, xt::all(), 0); xt::xtensor n1 = xt::view(conn, xt::all(), 1); xt::xtensor n2 = xt::view(conn, xt::all(), 2); xt::xtensor n3 = xt::view(conn, xt::all(), 3); xt::xtensor x0 = xt::view(coor, xt::keep(n0), 0); xt::xtensor x1 = xt::view(coor, xt::keep(n1), 0); xt::xtensor x2 = xt::view(coor, xt::keep(n2), 0); xt::xtensor x3 = xt::view(coor, xt::keep(n3), 0); xt::xtensor y0 = xt::view(coor, xt::keep(n0), 1); xt::xtensor y1 = xt::view(coor, xt::keep(n1), 1); xt::xtensor y2 = xt::view(coor, xt::keep(n2), 1); xt::xtensor y3 = xt::view(coor, xt::keep(n3), 1); xt::xtensor ret = xt::empty(conn.shape()); xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0)); xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0)); xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0)); xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0)); return ret; } throw std::runtime_error("Element-type not implemented"); } inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn) { return edgesize(coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); xt::xtensor ret = xt::zeros({conn.shape(0), coor.shape(1)}); if (type == ElementType::Quad4) { GOOSEFEM_ASSERT(coor.shape(1) == 2); GOOSEFEM_ASSERT(conn.shape(1) == 4); for (size_t i = 0; i < 4; ++i) { auto n = xt::view(conn, xt::all(), i); ret += xt::view(coor, xt::keep(n), xt::all()); } ret /= 4.0; return ret; } throw std::runtime_error("Element-type not implemented"); } inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn) { return centers(coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0)); size_t N = coor.shape(0); xt::xtensor ret = N * xt::ones({N}); if (type == ElementType::Quad4) { GOOSEFEM_ASSERT(coor.shape(1) == 2); GOOSEFEM_ASSERT(conn.shape(1) == 4); for (size_t i = 0; i < 4; ++i) { xt::xtensor t = N * xt::ones({N}); auto old_nd = xt::view(conn, xt::all(), i); auto new_nd = xt::view(conn, xt::keep(elem_map), i); xt::view(t, xt::keep(old_nd)) = new_nd; ret = xt::where(xt::equal(ret, N), t, ret); } return ret; } throw std::runtime_error("Element-type not implemented"); } inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn) { return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn)); } } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/MeshHex8.h b/include/GooseFEM/MeshHex8.h index 5a9d925..c185303 100644 --- a/include/GooseFEM/MeshHex8.h +++ b/include/GooseFEM/MeshHex8.h @@ -1,368 +1,183 @@ /** Generate simple meshes of 8-noded hexahedral elements in 3d (GooseFEM::Mesh::ElementType::Hex8). \file MeshHex8.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHHEX8_H #define GOOSEFEM_MESHHEX8_H #include "config.h" +#include "Mesh.h" namespace GooseFEM { namespace Mesh { namespace Hex8 { -class Regular { +/** +Regular mesh: equi-sized elements. +*/ +class Regular : public RegularBase3d { public: - Regular(size_t nelx, size_t nely, size_t nelz, double h = 1.); - - // 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 - - // type - ElementType getElementType() const; - - // mesh - xt::xtensor coor() const; // nodal positions [nnode, ndim] - xt::xtensor conn() const; // connectivity [nelem, nne] - - // boundary nodes: planes - xt::xtensor nodesFront() const; - xt::xtensor nodesBack() const; - xt::xtensor nodesLeft() const; - xt::xtensor nodesRight() const; - xt::xtensor nodesBottom() const; - xt::xtensor nodesTop() const; - - // boundary nodes: faces - xt::xtensor nodesFrontFace() const; - xt::xtensor nodesBackFace() const; - xt::xtensor nodesLeftFace() const; - xt::xtensor nodesRightFace() const; - xt::xtensor nodesBottomFace() const; - xt::xtensor nodesTopFace() const; - - // boundary nodes: edges - xt::xtensor nodesFrontBottomEdge() const; - xt::xtensor nodesFrontTopEdge() const; - xt::xtensor nodesFrontLeftEdge() const; - xt::xtensor nodesFrontRightEdge() const; - xt::xtensor nodesBackBottomEdge() const; - xt::xtensor nodesBackTopEdge() const; - xt::xtensor nodesBackLeftEdge() const; - xt::xtensor nodesBackRightEdge() const; - xt::xtensor nodesBottomLeftEdge() const; - xt::xtensor nodesBottomRightEdge() const; - xt::xtensor nodesTopLeftEdge() const; - xt::xtensor nodesTopRightEdge() const; - - // boundary nodes: edges (aliases) - xt::xtensor nodesBottomFrontEdge() const; - xt::xtensor nodesBottomBackEdge() const; - xt::xtensor nodesTopFrontEdge() const; - xt::xtensor nodesTopBackEdge() const; - xt::xtensor nodesLeftBottomEdge() const; - xt::xtensor nodesLeftFrontEdge() const; - xt::xtensor nodesLeftBackEdge() const; - xt::xtensor nodesLeftTopEdge() const; - xt::xtensor nodesRightBottomEdge() const; - xt::xtensor nodesRightTopEdge() const; - xt::xtensor nodesRightFrontEdge() const; - xt::xtensor nodesRightBackEdge() const; - - // boundary nodes: edges, without corners - xt::xtensor nodesFrontBottomOpenEdge() const; - xt::xtensor nodesFrontTopOpenEdge() const; - xt::xtensor nodesFrontLeftOpenEdge() const; - xt::xtensor nodesFrontRightOpenEdge() const; - xt::xtensor nodesBackBottomOpenEdge() const; - xt::xtensor nodesBackTopOpenEdge() const; - xt::xtensor nodesBackLeftOpenEdge() const; - xt::xtensor nodesBackRightOpenEdge() const; - xt::xtensor nodesBottomLeftOpenEdge() const; - xt::xtensor nodesBottomRightOpenEdge() const; - xt::xtensor nodesTopLeftOpenEdge() const; - xt::xtensor nodesTopRightOpenEdge() const; - - // boundary nodes: edges, without corners (aliases) - xt::xtensor nodesBottomFrontOpenEdge() const; - xt::xtensor nodesBottomBackOpenEdge() const; - xt::xtensor nodesTopFrontOpenEdge() const; - xt::xtensor nodesTopBackOpenEdge() const; - xt::xtensor nodesLeftBottomOpenEdge() const; - xt::xtensor nodesLeftFrontOpenEdge() const; - xt::xtensor nodesLeftBackOpenEdge() const; - xt::xtensor nodesLeftTopOpenEdge() const; - xt::xtensor nodesRightBottomOpenEdge() const; - xt::xtensor nodesRightTopOpenEdge() const; - xt::xtensor nodesRightFrontOpenEdge() const; - xt::xtensor nodesRightBackOpenEdge() const; - - // boundary nodes: corners - size_t nodesFrontBottomLeftCorner() const; - size_t nodesFrontBottomRightCorner() const; - size_t nodesFrontTopLeftCorner() const; - size_t nodesFrontTopRightCorner() const; - size_t nodesBackBottomLeftCorner() const; - size_t nodesBackBottomRightCorner() const; - size_t nodesBackTopLeftCorner() const; - size_t nodesBackTopRightCorner() const; - - // boundary nodes: corners (aliases) - size_t nodesFrontLeftBottomCorner() const; - size_t nodesBottomFrontLeftCorner() const; - size_t nodesBottomLeftFrontCorner() const; - size_t nodesLeftFrontBottomCorner() const; - size_t nodesLeftBottomFrontCorner() const; - size_t nodesFrontRightBottomCorner() const; - size_t nodesBottomFrontRightCorner() const; - size_t nodesBottomRightFrontCorner() const; - size_t nodesRightFrontBottomCorner() const; - size_t nodesRightBottomFrontCorner() const; - size_t nodesFrontLeftTopCorner() const; - size_t nodesTopFrontLeftCorner() const; - size_t nodesTopLeftFrontCorner() const; - size_t nodesLeftFrontTopCorner() const; - size_t nodesLeftTopFrontCorner() const; - size_t nodesFrontRightTopCorner() const; - size_t nodesTopFrontRightCorner() const; - size_t nodesTopRightFrontCorner() const; - size_t nodesRightFrontTopCorner() const; - size_t nodesRightTopFrontCorner() const; - size_t nodesBackLeftBottomCorner() const; - size_t nodesBottomBackLeftCorner() const; - size_t nodesBottomLeftBackCorner() const; - size_t nodesLeftBackBottomCorner() const; - size_t nodesLeftBottomBackCorner() const; - size_t nodesBackRightBottomCorner() const; - size_t nodesBottomBackRightCorner() const; - size_t nodesBottomRightBackCorner() const; - size_t nodesRightBackBottomCorner() const; - size_t nodesRightBottomBackCorner() const; - size_t nodesBackLeftTopCorner() const; - size_t nodesTopBackLeftCorner() const; - size_t nodesTopLeftBackCorner() const; - size_t nodesLeftBackTopCorner() const; - size_t nodesLeftTopBackCorner() const; - size_t nodesBackRightTopCorner() const; - size_t nodesTopBackRightCorner() const; - size_t nodesTopRightBackCorner() const; - size_t nodesRightBackTopCorner() const; - size_t nodesRightTopBackCorner() 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; -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_nelz; // number of elements in z-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 = 8; // number of nodes-per-element - static const size_t m_ndim = 3; // number of dimensions + /** + Constructor. + + \param nelx Number of elements in horizontal (x) direction. + \param nely Number of elements in vertical (y) direction. + \param nelz Number of elements in vertical (z) direction. + \param h Edge size (width == height == depth). + */ + Regular(size_t nelx, size_t nely, size_t nelz, double h = 1.0); + + ElementType getElementType() const override; + xt::xtensor coor() const override; + xt::xtensor conn() const override; + xt::xtensor nodesFront() const override; + xt::xtensor nodesBack() const override; + xt::xtensor nodesLeft() const override; + xt::xtensor nodesRight() const override; + xt::xtensor nodesBottom() const override; + xt::xtensor nodesTop() const override; + xt::xtensor nodesFrontFace() const override; + xt::xtensor nodesBackFace() const override; + xt::xtensor nodesLeftFace() const override; + xt::xtensor nodesRightFace() const override; + xt::xtensor nodesBottomFace() const override; + xt::xtensor nodesTopFace() const override; + xt::xtensor nodesFrontBottomEdge() const override; + xt::xtensor nodesFrontTopEdge() const override; + xt::xtensor nodesFrontLeftEdge() const override; + xt::xtensor nodesFrontRightEdge() const override; + xt::xtensor nodesBackBottomEdge() const override; + xt::xtensor nodesBackTopEdge() const override; + xt::xtensor nodesBackLeftEdge() const override; + xt::xtensor nodesBackRightEdge() const override; + xt::xtensor nodesBottomLeftEdge() const override; + xt::xtensor nodesBottomRightEdge() const override; + xt::xtensor nodesTopLeftEdge() const override; + xt::xtensor nodesTopRightEdge() const override; + xt::xtensor nodesFrontBottomOpenEdge() const override; + xt::xtensor nodesFrontTopOpenEdge() const override; + xt::xtensor nodesFrontLeftOpenEdge() const override; + xt::xtensor nodesFrontRightOpenEdge() const override; + xt::xtensor nodesBackBottomOpenEdge() const override; + xt::xtensor nodesBackTopOpenEdge() const override; + xt::xtensor nodesBackLeftOpenEdge() const override; + xt::xtensor nodesBackRightOpenEdge() const override; + xt::xtensor nodesBottomLeftOpenEdge() const override; + xt::xtensor nodesBottomRightOpenEdge() const override; + xt::xtensor nodesTopLeftOpenEdge() const override; + xt::xtensor nodesTopRightOpenEdge() const override; + size_t nodesFrontBottomLeftCorner() const override; + size_t nodesFrontBottomRightCorner() const override; + size_t nodesFrontTopLeftCorner() const override; + size_t nodesFrontTopRightCorner() const override; + size_t nodesBackBottomLeftCorner() const override; + size_t nodesBackBottomRightCorner() const override; + size_t nodesBackTopLeftCorner() const override; + size_t nodesBackTopRightCorner() const override; }; -class FineLayer { +/** +Mesh with fine middle layer, and coarser elements towards the top and bottom. +*/ +class FineLayer : public RegularBase3d { public: - FineLayer(size_t nelx, size_t nely, size_t nelz, double h = 1.0, size_t nfine = 1); - - // 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 - size_t nelz() const; // number of elements in y-direction - - // 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: planes - xt::xtensor nodesFront() const; - xt::xtensor nodesBack() const; - xt::xtensor nodesLeft() const; - xt::xtensor nodesRight() const; - xt::xtensor nodesBottom() const; - xt::xtensor nodesTop() const; - // boundary nodes: faces - xt::xtensor nodesFrontFace() const; - xt::xtensor nodesBackFace() const; - xt::xtensor nodesLeftFace() const; - xt::xtensor nodesRightFace() const; - xt::xtensor nodesBottomFace() const; - xt::xtensor nodesTopFace() const; + /** + Constructor. - // boundary nodes: edges - xt::xtensor nodesFrontBottomEdge() const; - xt::xtensor nodesFrontTopEdge() const; - xt::xtensor nodesFrontLeftEdge() const; - xt::xtensor nodesFrontRightEdge() const; - xt::xtensor nodesBackBottomEdge() const; - xt::xtensor nodesBackTopEdge() const; - xt::xtensor nodesBackLeftEdge() const; - xt::xtensor nodesBackRightEdge() const; - xt::xtensor nodesBottomLeftEdge() const; - xt::xtensor nodesBottomRightEdge() const; - xt::xtensor nodesTopLeftEdge() const; - xt::xtensor nodesTopRightEdge() const; + \param nelx Number of elements (along the middle layer) in horizontal (x) direction. + \param nely Approximate equivalent number of elements in vertical (y) direction. + \param nelz Number of elements (along the middle layer) in depth (z) direction. + \param h Edge size (width == height == depth) of elements along the weak layer. - // boundary nodes: edges (aliases) - xt::xtensor nodesBottomFrontEdge() const; - xt::xtensor nodesBottomBackEdge() const; - xt::xtensor nodesTopFrontEdge() const; - xt::xtensor nodesTopBackEdge() const; - xt::xtensor nodesLeftBottomEdge() const; - xt::xtensor nodesLeftFrontEdge() const; - xt::xtensor nodesLeftBackEdge() const; - xt::xtensor nodesLeftTopEdge() const; - xt::xtensor nodesRightBottomEdge() const; - xt::xtensor nodesRightTopEdge() const; - xt::xtensor nodesRightFrontEdge() const; - xt::xtensor nodesRightBackEdge() const; - - // boundary nodes: edges, without corners - xt::xtensor nodesFrontBottomOpenEdge() const; - xt::xtensor nodesFrontTopOpenEdge() const; - xt::xtensor nodesFrontLeftOpenEdge() const; - xt::xtensor nodesFrontRightOpenEdge() const; - xt::xtensor nodesBackBottomOpenEdge() const; - xt::xtensor nodesBackTopOpenEdge() const; - xt::xtensor nodesBackLeftOpenEdge() const; - xt::xtensor nodesBackRightOpenEdge() const; - xt::xtensor nodesBottomLeftOpenEdge() const; - xt::xtensor nodesBottomRightOpenEdge() const; - xt::xtensor nodesTopLeftOpenEdge() const; - xt::xtensor nodesTopRightOpenEdge() const; - - // boundary nodes: edges, without corners (aliases) - xt::xtensor nodesBottomFrontOpenEdge() const; - xt::xtensor nodesBottomBackOpenEdge() const; - xt::xtensor nodesTopFrontOpenEdge() const; - xt::xtensor nodesTopBackOpenEdge() const; - xt::xtensor nodesLeftBottomOpenEdge() const; - xt::xtensor nodesLeftFrontOpenEdge() const; - xt::xtensor nodesLeftBackOpenEdge() const; - xt::xtensor nodesLeftTopOpenEdge() const; - xt::xtensor nodesRightBottomOpenEdge() const; - xt::xtensor nodesRightTopOpenEdge() const; - xt::xtensor nodesRightFrontOpenEdge() const; - xt::xtensor nodesRightBackOpenEdge() const; - - // boundary nodes: corners - size_t nodesFrontBottomLeftCorner() const; - size_t nodesFrontBottomRightCorner() const; - size_t nodesFrontTopLeftCorner() const; - size_t nodesFrontTopRightCorner() const; - size_t nodesBackBottomLeftCorner() const; - size_t nodesBackBottomRightCorner() const; - size_t nodesBackTopLeftCorner() const; - size_t nodesBackTopRightCorner() const; - - // boundary nodes: corners (aliases) - size_t nodesFrontLeftBottomCorner() const; - size_t nodesBottomFrontLeftCorner() const; - size_t nodesBottomLeftFrontCorner() const; - size_t nodesLeftFrontBottomCorner() const; - size_t nodesLeftBottomFrontCorner() const; - size_t nodesFrontRightBottomCorner() const; - size_t nodesBottomFrontRightCorner() const; - size_t nodesBottomRightFrontCorner() const; - size_t nodesRightFrontBottomCorner() const; - size_t nodesRightBottomFrontCorner() const; - size_t nodesFrontLeftTopCorner() const; - size_t nodesTopFrontLeftCorner() const; - size_t nodesTopLeftFrontCorner() const; - size_t nodesLeftFrontTopCorner() const; - size_t nodesLeftTopFrontCorner() const; - size_t nodesFrontRightTopCorner() const; - size_t nodesTopFrontRightCorner() const; - size_t nodesTopRightFrontCorner() const; - size_t nodesRightFrontTopCorner() const; - size_t nodesRightTopFrontCorner() const; - size_t nodesBackLeftBottomCorner() const; - size_t nodesBottomBackLeftCorner() const; - size_t nodesBottomLeftBackCorner() const; - size_t nodesLeftBackBottomCorner() const; - size_t nodesLeftBottomBackCorner() const; - size_t nodesBackRightBottomCorner() const; - size_t nodesBottomBackRightCorner() const; - size_t nodesBottomRightBackCorner() const; - size_t nodesRightBackBottomCorner() const; - size_t nodesRightBottomBackCorner() const; - size_t nodesBackLeftTopCorner() const; - size_t nodesTopBackLeftCorner() const; - size_t nodesTopLeftBackCorner() const; - size_t nodesLeftBackTopCorner() const; - size_t nodesLeftTopBackCorner() const; - size_t nodesBackRightTopCorner() const; - size_t nodesTopBackRightCorner() const; - size_t nodesTopRightBackCorner() const; - size_t nodesRightBackTopCorner() const; - size_t nodesRightTopBackCorner() 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; + \param nfine + Extra number of fine layers around the middle layer. + By default the element size is kept smaller than the distance to the middle layer. + */ + FineLayer(size_t nelx, size_t nely, size_t nelz, double h = 1.0, size_t nfine = 1); - // front-bottom-left node, used as reference for periodicity - size_t nodesOrigin() const; + ElementType getElementType() const override; + xt::xtensor coor() const override; + xt::xtensor conn() const override; + xt::xtensor nodesFront() const override; + xt::xtensor nodesBack() const override; + xt::xtensor nodesLeft() const override; + xt::xtensor nodesRight() const override; + xt::xtensor nodesBottom() const override; + xt::xtensor nodesTop() const override; + xt::xtensor nodesFrontFace() const override; + xt::xtensor nodesBackFace() const override; + xt::xtensor nodesLeftFace() const override; + xt::xtensor nodesRightFace() const override; + xt::xtensor nodesBottomFace() const override; + xt::xtensor nodesTopFace() const override; + xt::xtensor nodesFrontBottomEdge() const override; + xt::xtensor nodesFrontTopEdge() const override; + xt::xtensor nodesFrontLeftEdge() const override; + xt::xtensor nodesFrontRightEdge() const override; + xt::xtensor nodesBackBottomEdge() const override; + xt::xtensor nodesBackTopEdge() const override; + xt::xtensor nodesBackLeftEdge() const override; + xt::xtensor nodesBackRightEdge() const override; + xt::xtensor nodesBottomLeftEdge() const override; + xt::xtensor nodesBottomRightEdge() const override; + xt::xtensor nodesTopLeftEdge() const override; + xt::xtensor nodesTopRightEdge() const override; + xt::xtensor nodesFrontBottomOpenEdge() const override; + xt::xtensor nodesFrontTopOpenEdge() const override; + xt::xtensor nodesFrontLeftOpenEdge() const override; + xt::xtensor nodesFrontRightOpenEdge() const override; + xt::xtensor nodesBackBottomOpenEdge() const override; + xt::xtensor nodesBackTopOpenEdge() const override; + xt::xtensor nodesBackLeftOpenEdge() const override; + xt::xtensor nodesBackRightOpenEdge() const override; + xt::xtensor nodesBottomLeftOpenEdge() const override; + xt::xtensor nodesBottomRightOpenEdge() const override; + xt::xtensor nodesTopLeftOpenEdge() const override; + xt::xtensor nodesTopRightOpenEdge() const override; + size_t nodesFrontBottomLeftCorner() const override; + size_t nodesFrontBottomRightCorner() const override; + size_t nodesFrontTopLeftCorner() const override; + size_t nodesFrontTopRightCorner() const override; + size_t nodesBackBottomLeftCorner() const override; + size_t nodesBackBottomRightCorner() const override; + size_t nodesBackTopLeftCorner() const override; + size_t nodesBackTopRightCorner() const override; + + size_t nelx() const override; + size_t nely() const override; + size_t nelz() const override; + + /** + Elements in the middle (fine) layer. + + \return List of element numbers. + */ + xt::xtensor elementsMiddleLayer() const; private: - double m_h; // elementary element edge-size (in all directions) - double m_Lx; // mesh size in "x" - double m_Lz; // mesh size in "z" - size_t m_nelem; // number of elements - size_t m_nnode; // number of nodes - static const size_t m_nne = 8; // number of nodes-per-element - static const size_t m_ndim = 3; // number of dimensions - xt::xtensor m_nelx; // number of elements in "x" (*) - xt::xtensor m_nelz; // number of elements in "z" (*) - xt::xtensor m_nnd; // 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_nhz; // element size in z-direction (*) - xt::xtensor m_refine; // refine direction (-1:no refine, 0:"x", 2:"z") (*) - xt::xtensor m_startElem; // start element (*) - xt::xtensor m_startNode; // start node (**) - // (*) per element layer in "y" - // (**) per node layer in "y" + double m_Lx; ///< mesh size in "x" + double m_Lz; ///< mesh size in "z" + xt::xtensor m_layer_nelx; ///< number of elements in "x" per element layer in "y" + xt::xtensor m_layer_nelz; ///< number of elements in "z" per element layer in "y" + xt::xtensor m_nnd; ///< number of nodes in the main node layer per node layer in "y" + xt::xtensor m_nhx; ///< element size in x-direction per element layer in "y" + xt::xtensor m_nhy; ///< element size in y-direction per element layer in "y" + xt::xtensor m_nhz; ///< element size in z-direction per element layer in "y" + xt::xtensor m_refine; ///< refine direction (-1:no refine, 0:"x", 2:"z") per element layer in "y" + xt::xtensor m_startElem; ///< start element per element layer in "y" + xt::xtensor m_startNode; ///< start node per node layer in "y" }; } // namespace Hex8 } // namespace Mesh } // namespace GooseFEM #include "MeshHex8.hpp" #endif diff --git a/include/GooseFEM/MeshHex8.hpp b/include/GooseFEM/MeshHex8.hpp index 0e340ea..3d638a6 100644 --- a/include/GooseFEM/MeshHex8.hpp +++ b/include/GooseFEM/MeshHex8.hpp @@ -1,3017 +1,2189 @@ /** Implementation of MeshHex8.h \file MeshHex8.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHHEX8_HPP #define GOOSEFEM_MESHHEX8_HPP #include "MeshHex8.h" namespace GooseFEM { namespace Mesh { namespace Hex8 { inline Regular::Regular(size_t nelx, size_t nely, size_t nelz, double h) - : m_h(h), m_nelx(nelx), m_nely(nely), m_nelz(nelz) { + m_h = h; + m_nelx = nelx; + m_nely = nely; + m_nelz = nelz; + m_nne = 8; + m_ndim = 3; + GOOSEFEM_ASSERT(m_nelx >= 1ul); GOOSEFEM_ASSERT(m_nely >= 1ul); GOOSEFEM_ASSERT(m_nelz >= 1ul); m_nnode = (m_nelx + 1) * (m_nely + 1) * (m_nelz + 1); m_nelem = m_nelx * m_nely * m_nelz; } -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 ElementType Regular::getElementType() const { return ElementType::Hex8; } 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); xt::xtensor z = xt::linspace(0.0, m_h * static_cast(m_nelz), m_nelz + 1); size_t inode = 0; for (size_t iz = 0; iz < m_nelz + 1; ++iz) { 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); ret(inode, 2) = z(iz); ++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 iz = 0; iz < m_nelz; ++iz) { 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 + iz * (m_nely + 1) * (m_nelx + 1); ret(ielem, 1) = iy * (m_nelx + 1) + (ix + 1) + iz * (m_nely + 1) * (m_nelx + 1); ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + ix + iz * (m_nely + 1) * (m_nelx + 1); ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1) + iz * (m_nely + 1) * (m_nelx + 1); ret(ielem, 4) = iy * (m_nelx + 1) + ix + (iz + 1) * (m_nely + 1) * (m_nelx + 1); ret(ielem, 5) = (iy) * (m_nelx + 1) + (ix + 1) + (iz + 1) * (m_nely + 1) * (m_nelx + 1); ret(ielem, 7) = (iy + 1) * (m_nelx + 1) + ix + (iz + 1) * (m_nely + 1) * (m_nelx + 1); ret(ielem, 6) = (iy + 1) * (m_nelx + 1) + (ix + 1) + (iz + 1) * (m_nely + 1) * (m_nelx + 1); ++ielem; } } } return ret; } inline xt::xtensor Regular::nodesFront() const { xt::xtensor ret = xt::empty({(m_nelx + 1) * (m_nely + 1)}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(iy * (m_nelx + 1) + ix) = iy * (m_nelx + 1) + ix; } } return ret; } inline xt::xtensor Regular::nodesBack() const { xt::xtensor ret = xt::empty({(m_nelx + 1) * (m_nely + 1)}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(iy * (m_nelx + 1) + ix) = iy * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1); } } return ret; } inline xt::xtensor Regular::nodesLeft() const { xt::xtensor ret = xt::empty({(m_nely + 1) * (m_nelz + 1)}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iz * (m_nely + 1) + iy) = iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1); } } return ret; } inline xt::xtensor Regular::nodesRight() const { xt::xtensor ret = xt::empty({(m_nely + 1) * (m_nelz + 1)}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iz * (m_nely + 1) + iy) = iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } } return ret; } inline xt::xtensor Regular::nodesBottom() const { xt::xtensor ret = xt::empty({(m_nelx + 1) * (m_nelz + 1)}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(iz * (m_nelx + 1) + ix) = ix + iz * (m_nelx + 1) * (m_nely + 1); } } return ret; } inline xt::xtensor Regular::nodesTop() const { xt::xtensor ret = xt::empty({(m_nelx + 1) * (m_nelz + 1)}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(iz * (m_nelx + 1) + ix) = ix + m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1); } } return ret; } inline xt::xtensor Regular::nodesFrontFace() const { xt::xtensor ret = xt::empty({(m_nelx - 1) * (m_nely - 1)}); for (size_t iy = 1; iy < m_nely; ++iy) { for (size_t ix = 1; ix < m_nelx; ++ix) { ret((iy - 1) * (m_nelx - 1) + (ix - 1)) = iy * (m_nelx + 1) + ix; } } return ret; } inline xt::xtensor Regular::nodesBackFace() const { xt::xtensor ret = xt::empty({(m_nelx - 1) * (m_nely - 1)}); for (size_t iy = 1; iy < m_nely; ++iy) { for (size_t ix = 1; ix < m_nelx; ++ix) { ret((iy - 1) * (m_nelx - 1) + (ix - 1)) = iy * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1); } } return ret; } inline xt::xtensor Regular::nodesLeftFace() const { xt::xtensor ret = xt::empty({(m_nely - 1) * (m_nelz - 1)}); for (size_t iz = 1; iz < m_nelz; ++iz) { for (size_t iy = 1; iy < m_nely; ++iy) { ret((iz - 1) * (m_nely - 1) + (iy - 1)) = iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1); } } return ret; } inline xt::xtensor Regular::nodesRightFace() const { xt::xtensor ret = xt::empty({(m_nely - 1) * (m_nelz - 1)}); for (size_t iz = 1; iz < m_nelz; ++iz) { for (size_t iy = 1; iy < m_nely; ++iy) { ret((iz - 1) * (m_nely - 1) + (iy - 1)) = iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } } return ret; } inline xt::xtensor Regular::nodesBottomFace() const { xt::xtensor ret = xt::empty({(m_nelx - 1) * (m_nelz - 1)}); for (size_t iz = 1; iz < m_nelz; ++iz) { for (size_t ix = 1; ix < m_nelx; ++ix) { ret((iz - 1) * (m_nelx - 1) + (ix - 1)) = ix + iz * (m_nelx + 1) * (m_nely + 1); } } return ret; } inline xt::xtensor Regular::nodesTopFace() const { xt::xtensor ret = xt::empty({(m_nelx - 1) * (m_nelz - 1)}); for (size_t iz = 1; iz < m_nelz; ++iz) { for (size_t ix = 1; ix < m_nelx; ++ix) { ret((iz - 1) * (m_nelx - 1) + (ix - 1)) = ix + m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1); } } return ret; } inline xt::xtensor Regular::nodesFrontBottomEdge() const { xt::xtensor ret = xt::empty({m_nelx + 1}); for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(ix) = ix; } return ret; } inline xt::xtensor Regular::nodesFrontTopEdge() const { xt::xtensor ret = xt::empty({m_nelx + 1}); for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(ix) = ix + m_nely * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesFrontLeftEdge() const { xt::xtensor ret = xt::empty({m_nely + 1}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iy) = iy * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesFrontRightEdge() const { xt::xtensor ret = xt::empty({m_nely + 1}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iy) = iy * (m_nelx + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesBackBottomEdge() const { xt::xtensor ret = xt::empty({m_nelx + 1}); for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(ix) = ix + m_nelz * (m_nely + 1) * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesBackTopEdge() const { xt::xtensor ret = xt::empty({m_nelx + 1}); for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(ix) = m_nely * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesBackLeftEdge() const { xt::xtensor ret = xt::empty({m_nely + 1}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iy) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1); } return ret; } inline xt::xtensor Regular::nodesBackRightEdge() const { xt::xtensor ret = xt::empty({m_nely + 1}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iy) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesBottomLeftEdge() const { xt::xtensor ret = xt::empty({m_nelz + 1}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { ret(iz) = iz * (m_nelx + 1) * (m_nely + 1); } return ret; } inline xt::xtensor Regular::nodesBottomRightEdge() const { xt::xtensor ret = xt::empty({m_nelz + 1}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { ret(iz) = iz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesTopLeftEdge() const { xt::xtensor ret = xt::empty({m_nelz + 1}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { ret(iz) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1); } return ret; } inline xt::xtensor Regular::nodesTopRightEdge() const { xt::xtensor ret = xt::empty({m_nelz + 1}); for (size_t iz = 0; iz < m_nelz + 1; ++iz) { ret(iz) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } return ret; } -inline xt::xtensor Regular::nodesBottomFrontEdge() const -{ - return nodesFrontBottomEdge(); -} -inline xt::xtensor Regular::nodesBottomBackEdge() const -{ - return nodesBackBottomEdge(); -} -inline xt::xtensor Regular::nodesTopFrontEdge() const -{ - return nodesFrontTopEdge(); -} -inline xt::xtensor Regular::nodesTopBackEdge() const -{ - return nodesBackTopEdge(); -} -inline xt::xtensor Regular::nodesLeftBottomEdge() const -{ - return nodesBottomLeftEdge(); -} -inline xt::xtensor Regular::nodesLeftFrontEdge() const -{ - return nodesFrontLeftEdge(); -} -inline xt::xtensor Regular::nodesLeftBackEdge() const -{ - return nodesBackLeftEdge(); -} -inline xt::xtensor Regular::nodesLeftTopEdge() const -{ - return nodesTopLeftEdge(); -} -inline xt::xtensor Regular::nodesRightBottomEdge() const -{ - return nodesBottomRightEdge(); -} -inline xt::xtensor Regular::nodesRightTopEdge() const -{ - return nodesTopRightEdge(); -} -inline xt::xtensor Regular::nodesRightFrontEdge() const -{ - return nodesFrontRightEdge(); -} -inline xt::xtensor Regular::nodesRightBackEdge() const -{ - return nodesBackRightEdge(); -} - inline xt::xtensor Regular::nodesFrontBottomOpenEdge() const { xt::xtensor ret = xt::empty({m_nelx - 1}); for (size_t ix = 1; ix < m_nelx; ++ix) { ret(ix - 1) = ix; } return ret; } inline xt::xtensor Regular::nodesFrontTopOpenEdge() const { xt::xtensor ret = xt::empty({m_nelx - 1}); for (size_t ix = 1; ix < m_nelx; ++ix) { ret(ix - 1) = ix + m_nely * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesFrontLeftOpenEdge() const { xt::xtensor ret = xt::empty({m_nely - 1}); for (size_t iy = 1; iy < m_nely; ++iy) { ret(iy - 1) = iy * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesFrontRightOpenEdge() const { xt::xtensor ret = xt::empty({m_nely - 1}); for (size_t iy = 1; iy < m_nely; ++iy) { ret(iy - 1) = iy * (m_nelx + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesBackBottomOpenEdge() const { xt::xtensor ret = xt::empty({m_nelx - 1}); for (size_t ix = 1; ix < m_nelx; ++ix) { ret(ix - 1) = ix + m_nelz * (m_nely + 1) * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesBackTopOpenEdge() const { xt::xtensor ret = xt::empty({m_nelx - 1}); for (size_t ix = 1; ix < m_nelx; ++ix) { ret(ix - 1) = m_nely * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesBackLeftOpenEdge() const { xt::xtensor ret = xt::empty({m_nely - 1}); for (size_t iy = 1; iy < m_nely; ++iy) { ret(iy - 1) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1); } return ret; } inline xt::xtensor Regular::nodesBackRightOpenEdge() const { xt::xtensor ret = xt::empty({m_nely - 1}); for (size_t iy = 1; iy < m_nely; ++iy) { ret(iy - 1) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesBottomLeftOpenEdge() const { xt::xtensor ret = xt::empty({m_nelz - 1}); for (size_t iz = 1; iz < m_nelz; ++iz) { ret(iz - 1) = iz * (m_nelx + 1) * (m_nely + 1); } return ret; } inline xt::xtensor Regular::nodesBottomRightOpenEdge() const { xt::xtensor ret = xt::empty({m_nelz - 1}); for (size_t iz = 1; iz < m_nelz; ++iz) { ret(iz - 1) = iz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesTopLeftOpenEdge() const { xt::xtensor ret = xt::empty({m_nelz - 1}); for (size_t iz = 1; iz < m_nelz; ++iz) { ret(iz - 1) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1); } return ret; } inline xt::xtensor Regular::nodesTopRightOpenEdge() const { xt::xtensor ret = xt::empty({m_nelz - 1}); for (size_t iz = 1; iz < m_nelz; ++iz) { ret(iz - 1) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx; } return ret; } -inline xt::xtensor Regular::nodesBottomFrontOpenEdge() const -{ - return nodesFrontBottomOpenEdge(); -} -inline xt::xtensor Regular::nodesBottomBackOpenEdge() const -{ - return nodesBackBottomOpenEdge(); -} -inline xt::xtensor Regular::nodesTopFrontOpenEdge() const -{ - return nodesFrontTopOpenEdge(); -} -inline xt::xtensor Regular::nodesTopBackOpenEdge() const -{ - return nodesBackTopOpenEdge(); -} -inline xt::xtensor Regular::nodesLeftBottomOpenEdge() const -{ - return nodesBottomLeftOpenEdge(); -} -inline xt::xtensor Regular::nodesLeftFrontOpenEdge() const -{ - return nodesFrontLeftOpenEdge(); -} -inline xt::xtensor Regular::nodesLeftBackOpenEdge() const -{ - return nodesBackLeftOpenEdge(); -} -inline xt::xtensor Regular::nodesLeftTopOpenEdge() const -{ - return nodesTopLeftOpenEdge(); -} -inline xt::xtensor Regular::nodesRightBottomOpenEdge() const -{ - return nodesBottomRightOpenEdge(); -} -inline xt::xtensor Regular::nodesRightTopOpenEdge() const -{ - return nodesTopRightOpenEdge(); -} -inline xt::xtensor Regular::nodesRightFrontOpenEdge() const -{ - return nodesFrontRightOpenEdge(); -} -inline xt::xtensor Regular::nodesRightBackOpenEdge() const -{ - return nodesBackRightOpenEdge(); -} - inline size_t Regular::nodesFrontBottomLeftCorner() const { return 0; } inline size_t Regular::nodesFrontBottomRightCorner() const { return m_nelx; } inline size_t Regular::nodesFrontTopLeftCorner() const { return m_nely * (m_nelx + 1); } inline size_t Regular::nodesFrontTopRightCorner() const { return m_nely * (m_nelx + 1) + m_nelx; } inline size_t Regular::nodesBackBottomLeftCorner() const { return m_nelz * (m_nely + 1) * (m_nelx + 1); } inline size_t Regular::nodesBackBottomRightCorner() const { return m_nelx + m_nelz * (m_nely + 1) * (m_nelx + 1); } inline size_t Regular::nodesBackTopLeftCorner() const { return m_nely * (m_nelx + 1) + m_nelz * (m_nely + 1) * (m_nelx + 1); } inline size_t Regular::nodesBackTopRightCorner() const { return m_nely * (m_nelx + 1) + m_nelx + m_nelz * (m_nely + 1) * (m_nelx + 1); } -inline size_t Regular::nodesFrontLeftBottomCorner() const +inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine) { - return nodesFrontBottomLeftCorner(); -} -inline size_t Regular::nodesBottomFrontLeftCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t Regular::nodesBottomLeftFrontCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t Regular::nodesLeftFrontBottomCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t Regular::nodesLeftBottomFrontCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t Regular::nodesFrontRightBottomCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t Regular::nodesBottomFrontRightCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t Regular::nodesBottomRightFrontCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t Regular::nodesRightFrontBottomCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t Regular::nodesRightBottomFrontCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t Regular::nodesFrontLeftTopCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t Regular::nodesTopFrontLeftCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t Regular::nodesTopLeftFrontCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t Regular::nodesLeftFrontTopCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t Regular::nodesLeftTopFrontCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t Regular::nodesFrontRightTopCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t Regular::nodesTopFrontRightCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t Regular::nodesTopRightFrontCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t Regular::nodesRightFrontTopCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t Regular::nodesRightTopFrontCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t Regular::nodesBackLeftBottomCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t Regular::nodesBottomBackLeftCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t Regular::nodesBottomLeftBackCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t Regular::nodesLeftBackBottomCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t Regular::nodesLeftBottomBackCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t Regular::nodesBackRightBottomCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t Regular::nodesBottomBackRightCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t Regular::nodesBottomRightBackCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t Regular::nodesRightBackBottomCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t Regular::nodesRightBottomBackCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t Regular::nodesBackLeftTopCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t Regular::nodesTopBackLeftCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t Regular::nodesTopLeftBackCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t Regular::nodesLeftBackTopCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t Regular::nodesLeftTopBackCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t Regular::nodesBackRightTopCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t Regular::nodesTopBackRightCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t Regular::nodesTopRightBackCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t Regular::nodesRightBackTopCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t Regular::nodesRightTopBackCorner() const -{ - return nodesBackTopRightCorner(); -} - -inline xt::xtensor Regular::nodesPeriodic() const -{ - xt::xtensor fro = nodesFrontFace(); - xt::xtensor bck = nodesBackFace(); - xt::xtensor lft = nodesLeftFace(); - xt::xtensor rgt = nodesRightFace(); - xt::xtensor bot = nodesBottomFace(); - xt::xtensor top = nodesTopFace(); - - xt::xtensor froBot = nodesFrontBottomOpenEdge(); - xt::xtensor froTop = nodesFrontTopOpenEdge(); - xt::xtensor froLft = nodesFrontLeftOpenEdge(); - xt::xtensor froRgt = nodesFrontRightOpenEdge(); - xt::xtensor bckBot = nodesBackBottomOpenEdge(); - xt::xtensor bckTop = nodesBackTopOpenEdge(); - xt::xtensor bckLft = nodesBackLeftOpenEdge(); - xt::xtensor bckRgt = nodesBackRightOpenEdge(); - xt::xtensor botLft = nodesBottomLeftOpenEdge(); - xt::xtensor botRgt = nodesBottomRightOpenEdge(); - xt::xtensor topLft = nodesTopLeftOpenEdge(); - xt::xtensor topRgt = nodesTopRightOpenEdge(); - - size_t tface = fro.size() + lft.size() + bot.size(); - size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size(); - size_t tnode = 7; - xt::xtensor ret = xt::empty({tface + tedge + tnode, std::size_t(2)}); - - size_t i = 0; - - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesFrontBottomRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackBottomRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackBottomLeftCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesFrontTopLeftCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesFrontTopRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackTopRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackTopLeftCorner(); - ++i; - - for (size_t j = 0; j < froBot.size(); ++j) { - ret(i, 0) = froBot(j); - ret(i, 1) = bckBot(j); - ++i; - } - for (size_t j = 0; j < froBot.size(); ++j) { - ret(i, 0) = froBot(j); - ret(i, 1) = bckTop(j); - ++i; - } - for (size_t j = 0; j < froBot.size(); ++j) { - ret(i, 0) = froBot(j); - ret(i, 1) = froTop(j); - ++i; - } - for (size_t j = 0; j < botLft.size(); ++j) { - ret(i, 0) = botLft(j); - ret(i, 1) = botRgt(j); - ++i; - } - for (size_t j = 0; j < botLft.size(); ++j) { - ret(i, 0) = botLft(j); - ret(i, 1) = topRgt(j); - ++i; - } - for (size_t j = 0; j < botLft.size(); ++j) { - ret(i, 0) = botLft(j); - ret(i, 1) = topLft(j); - ++i; - } - for (size_t j = 0; j < froLft.size(); ++j) { - ret(i, 0) = froLft(j); - ret(i, 1) = froRgt(j); - ++i; - } - for (size_t j = 0; j < froLft.size(); ++j) { - ret(i, 0) = froLft(j); - ret(i, 1) = bckRgt(j); - ++i; - } - for (size_t j = 0; j < froLft.size(); ++j) { - ret(i, 0) = froLft(j); - ret(i, 1) = bckLft(j); - ++i; - } - - for (size_t j = 0; j < fro.size(); ++j) { - ret(i, 0) = fro(j); - ret(i, 1) = bck(j); - ++i; - } - for (size_t j = 0; j < lft.size(); ++j) { - ret(i, 0) = lft(j); - ret(i, 1) = rgt(j); - ++i; - } - for (size_t j = 0; j < bot.size(); ++j) { - ret(i, 0) = bot(j); - ret(i, 1) = top(j); - ++i; - } - - return ret; -} - -inline size_t Regular::nodesOrigin() const -{ - return nodesFrontBottomLeftCorner(); -} + m_h = h; + m_nne = 8; + m_ndim = 3; -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(); - - for (size_t i = 0; i < nodePer.shape(0); ++i) { - for (size_t j = 0; j < m_ndim; ++j) { - ret(nodePer(i, 1), j) = ret(nodePer(i, 0), j); - } - } - - return GooseFEM::Mesh::renumber(ret); -} - -inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine) : m_h(h) -{ // basic assumptions GOOSEFEM_ASSERT(nelx >= 1ul); GOOSEFEM_ASSERT(nely >= 1ul); GOOSEFEM_ASSERT(nelz >= 1ul); // store basic info m_Lx = m_h * static_cast(nelx); m_Lz = m_h * static_cast(nelz); // 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 nhz = 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-(z-)direction should fit the total number of elements in // x-(z-)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 (and z-direction) if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) { // - process refinement in x-direction 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; // - rule (2) satisfied: coarsen next element layer in z-direction if (iy + 1 < nely && ntot + 2 * nhy(iy) < nmin) { if (nelz % (3 * nhz(iy + 1)) == 0) { // - update the number of elements in y-direction ntot += nhy(iy); // - proceed to next element layer in y-direction ++iy; // - process refinement in z-direction refine(iy) = 2; nhy(iy) = nhy(iy - 1); auto vnhz = xt::view(nhz, xt::range(iy, _)); vnhz *= 3; } } } // rules (1,2) satisfied: coarse in z-direction else if (3 * nhy(iy) <= ntot && nelz % (3 * nhz(iy)) == 0 && ntot + nhy(iy) < nmin) { // - process refinement in z-direction refine(iy) = 2; nhy(iy) *= 2; auto vnhy = xt::view(nhy, xt::range(iy + 1, _)); auto vnhz = xt::view(nhz, xt::range(iy, _)); vnhy *= 3; vnhz *= 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_nhz = xt::empty({nely * 2 - 1}); m_refine = xt::empty({nely * 2 - 1}); m_layer_nelx = xt::empty({nely * 2 - 1}); m_layer_nelz = 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_nhz(iy) = nhz(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_nhz(iy + nely) = nhz(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_layer_nelx(iy) = nelx / m_nhx(iy); m_layer_nelz(iy) = nelz / m_nhz(iy); } // compute the number of nodes per node layer in y-direction // - bottom half for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) m_nnd(iy) = (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy) + 1); // - top half for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) m_nnd(iy + 1) = (m_layer_nelx(iy) + 1) * (m_layer_nelz(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_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1); } else if (m_refine(i) == 2) { m_nnode += (m_layer_nelx(i) + 1) * (3 * m_layer_nelz(i) + 1); } else { m_nnode += (m_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_layer_nelx(i)) * (m_layer_nelz(i)); } else if (m_refine(i) == 2) { m_nelem += (m_layer_nelx(i)) * (4 * m_layer_nelz(i)); } else { m_nelem += (m_layer_nelx(i)) * (m_layer_nelz(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_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1); } else if (m_refine(i) == 2) { m_nnode += (m_layer_nelx(i) + 1) * (5 * m_layer_nelz(i) + 1); } else { m_nnode += (m_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_layer_nelx(i)) * (m_layer_nelz(i)); } else if (m_refine(i) == 2) { m_nelem += (m_layer_nelx(i)) * (4 * m_layer_nelz(i)); } else { m_nelem += (m_layer_nelx(i)) * (m_layer_nelz(i)); } // - store the starting node of the next layer m_startNode(i + 1) = m_nnode; } // - add the top row of nodes m_nnode += (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(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_layer_nelx)(); } inline size_t FineLayer::nely() const { return xt::sum(m_nhy)(); } inline size_t FineLayer::nelz() const { return xt::amax(m_layer_nelz)(); } inline ElementType FineLayer::getElementType() const { return ElementType::Hex8; } 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_layer_nelx(iy) + 1); xt::xtensor z = xt::linspace(0.0, m_Lz, m_layer_nelz(iy) + 1); // add nodes of the bottom layer of this element for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ret(inode, 2) = z(iz); ++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 iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { for (size_t ix = 0; ix < m_layer_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; ret(inode, 2) = z(iz); ++inode; } } } } // add extra nodes of the intermediate layer, for refinement in z-direction else if (m_refine(iy) == 2) { // - get position offset in y- and z-direction double dz = m_h * static_cast(m_nhz(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t j = 0; j < 2; ++j) { for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy) + dy; ret(inode, 2) = z(iz) + dz * static_cast(j + 1); ++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_layer_nelx(iy) + 1); xt::xtensor z = xt::linspace(0.0, m_Lz, m_layer_nelz(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 iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { for (size_t ix = 0; ix < m_layer_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; ret(inode, 2) = z(iz); ++inode; } } } } // add extra nodes of the intermediate layer, for refinement in z-direction else if (m_refine(iy) == 2) { // - get position offset in y- and z-direction double dz = m_h * static_cast(m_nhz(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t j = 0; j < 2; ++j) { for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy) + dy; ret(inode, 2) = z(iz) + dz * static_cast(j + 1); ++inode; } } } } // add nodes of the top layer of this element for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy + 1); ret(inode, 2) = z(iz); ++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 iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) { ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1); 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 iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) { // -- bottom element ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 2) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 3) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 7) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy)); ielem++; // -- top-right element ret(ielem, 0) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 1) = top + (3 * ix + 3) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 3) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 4) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = top + (3 * ix + 3) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 7) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy)); ielem++; // -- top-center element ret(ielem, 0) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 1) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 2) = top + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 3) = top + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 4) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 5) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 6) = top + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 7) = top + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ielem++; // -- top-left element ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 1) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 2) = top + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 3) = top + (3 * ix) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 6) = top + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 7) = top + (3 * ix) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); 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 iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) { // -- lower-left element ret(ielem, 0) = bot + (3 * ix) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 2) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + (3 * ix) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 6) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ielem++; // -- lower-center element ret(ielem, 0) = bot + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 2) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 3) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 4) = bot + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 6) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 7) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy)); ielem++; // -- lower-right element ret(ielem, 0) = bot + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (3 * ix + 3) + iz * (3 * m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 3) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 4) = bot + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (3 * ix + 3) + (iz + 1) * (3 * m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy)); ielem++; // -- upper element ret(ielem, 0) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 1) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy)); ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 4) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 5) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy)); ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ielem++; } } } // - define connectivity: refinement along the z-direction (below the middle layer) else if (m_refine(iy) == 2 && iy <= (nely - 1) / 2) { for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) { // -- bottom element ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 3) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 4) = mid + ix + 2 * iz * (m_layer_nelx(iy) + 1); ret(ielem, 5) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = mid + (ix + 1) + 2 * iz * (m_layer_nelx(iy) + 1); ielem++; // -- top-back element ret(ielem, 0) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 1) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (3 * iz + 3) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (3 * iz + 3) * (m_layer_nelx(iy) + 1); ielem++; // -- top-center element ret(ielem, 0) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 1) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 4) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ielem++; // -- top-front element ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + (3 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + (3 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 4) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ielem++; } } } // - define connectivity: coarsening along the z-direction (above the middle layer) else if (m_refine(iy) == 2 && iy > (nely - 1) / 2) { for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) { for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) { // -- bottom-front element ret(ielem, 0) = bot + ix + (3 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (ix + 1) + (3 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1); ielem++; // -- bottom-center element ret(ielem, 0) = bot + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 3) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ielem++; // -- bottom-back element ret(ielem, 0) = bot + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 1) = bot + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 3) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 4) = bot + ix + (3 * iz + 3) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 3) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ielem++; // -- top element ret(ielem, 0) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 1) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1); ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1); ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1); ret(ielem, 4) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 5) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1); ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1); ielem++; } } } } return ret; } inline xt::xtensor FineLayer::elementsMiddleLayer() const { size_t nely = static_cast(m_nhy.size()); size_t iy = (nely - 1) / 2; xt::xtensor ret = xt::empty({m_layer_nelx(iy) * m_layer_nelz(iy)}); for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) { for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) { ret(ix + iz * m_layer_nelx(iy)) = m_startElem(iy) + ix + iz * m_layer_nelx(iy); } } return ret; } inline xt::xtensor FineLayer::nodesFront() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 + 1; } else { n += m_layer_nelx(iy) + 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 + 1; } else { n += m_layer_nelx(iy) + 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(j) = m_startNode(iy) + ix; ++j; } // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } // -- top node layer for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(j) = m_startNode(iy + 1) + ix; ++j; } } return ret; } inline xt::xtensor FineLayer::nodesBack() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 + 1; } else { n += m_layer_nelx(iy) + 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 + 1; } else { n += m_layer_nelx(iy) + 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(j) = m_startNode(iy) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy); ++j; } // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_nelz(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_nelz(iy); ++j; } } // -- top node layer for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) { ret(j) = m_startNode(iy + 1) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesLeft() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 2) { n += m_layer_nelz(iy) * 3 + 1; } else { n += m_layer_nelz(iy) + 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { if (m_refine(iy) == 2) { n += m_layer_nelz(iy) * 3 + 1; } else { n += m_layer_nelz(iy) + 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1); ++j; } // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_nnd(iy); ++j; } } // -- top node layer for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesRight() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 2) n += m_layer_nelz(iy) * 3 + 1; else n += m_layer_nelz(iy) + 1; } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { if (m_refine(iy) == 2) n += m_layer_nelz(iy) * 3 + 1; else n += m_layer_nelz(iy) + 1; } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } } // -- top node layer for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) { ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesBottom() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // allocate node list xt::xtensor ret = xt::empty({m_nnd(nely)}); // counter size_t j = 0; // fill node list for (size_t ix = 0; ix < m_layer_nelx(0) + 1; ++ix) { for (size_t iz = 0; iz < m_layer_nelz(0) + 1; ++iz) { ret(j) = m_startNode(0) + ix + iz * (m_layer_nelx(0) + 1); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesTop() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // allocate node list xt::xtensor ret = xt::empty({m_nnd(nely)}); // counter size_t j = 0; // fill node list for (size_t ix = 0; ix < m_layer_nelx(nely - 1) + 1; ++ix) { for (size_t iz = 0; iz < m_layer_nelz(nely - 1) + 1; ++iz) { ret(j) = m_startNode(nely) + ix + iz * (m_layer_nelx(nely - 1) + 1); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesFrontFace() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 - 1; } else { n += m_layer_nelx(iy) - 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 - 1; } else { n += m_layer_nelx(iy) - 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix; ++j; } // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } // -- top node layer for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy + 1) + ix; ++j; } } return ret; } inline xt::xtensor FineLayer::nodesBackFace() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 - 1; } else { n += m_layer_nelx(iy) - 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { if (m_refine(iy) == 0) { n += m_layer_nelx(iy) * 3 - 1; } else { n += m_layer_nelx(iy) - 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy); ++j; } // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_nelz(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { // -- refinement layer if (m_refine(iy) == 0) { for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_nelz(iy); ++j; } } // -- top node layer for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) { ret(j) = m_startNode(iy + 1) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesLeftFace() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 2) { n += m_layer_nelz(iy) * 3 - 1; } else { n += m_layer_nelz(iy) - 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { if (m_refine(iy) == 2) { n += m_layer_nelz(iy) * 3 - 1; } else { n += m_layer_nelz(iy) - 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1); ++j; } // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_nnd(iy); ++j; } } // -- top node layer for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesRightFace() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { if (m_refine(iy) == 2) { n += m_layer_nelz(iy) * 3 - 1; } else { n += m_layer_nelz(iy) - 1; } } // - top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { if (m_refine(iy) == 2) { n += m_layer_nelz(iy) * 3 - 1; } else { n += m_layer_nelz(iy) - 1; } } // allocate node-list xt::xtensor ret = xt::empty({n}); // initialize counter: current index in the node-list "ret" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { // -- bottom node layer for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } } } // top half: (middle node layer +) top node layer for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { // -- refinement layer if (m_refine(iy) == 2) { for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } } // -- top node layer for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) { ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesBottomFace() const { // allocate node list xt::xtensor ret = xt::empty({(m_layer_nelx(0) - 1) * (m_layer_nelz(0) - 1)}); // counter size_t j = 0; // fill node list for (size_t ix = 1; ix < m_layer_nelx(0); ++ix) { for (size_t iz = 1; iz < m_layer_nelz(0); ++iz) { ret(j) = m_startNode(0) + ix + iz * (m_layer_nelx(0) + 1); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesTopFace() const { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // allocate node list xt::xtensor ret = xt::empty({(m_layer_nelx(nely - 1) - 1) * (m_layer_nelz(nely - 1) - 1)}); // counter size_t j = 0; // fill node list for (size_t ix = 1; ix < m_layer_nelx(nely - 1); ++ix) { for (size_t iz = 1; iz < m_layer_nelz(nely - 1); ++iz) { ret(j) = m_startNode(nely) + ix + iz * (m_layer_nelx(nely - 1) + 1); ++j; } } return ret; } inline xt::xtensor FineLayer::nodesFrontBottomEdge() const { xt::xtensor ret = xt::empty({m_layer_nelx(0) + 1}); for (size_t ix = 0; ix < m_layer_nelx(0) + 1; ++ix) { ret(ix) = m_startNode(0) + ix; } return ret; } inline xt::xtensor FineLayer::nodesFrontTopEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelx(nely - 1) + 1}); for (size_t ix = 0; ix < m_layer_nelx(nely - 1) + 1; ++ix) { ret(ix) = m_startNode(nely) + ix; } return ret; } inline xt::xtensor FineLayer::nodesFrontLeftEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely + 1}); for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { ret(iy) = m_startNode(iy); } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { ret(iy + 1) = m_startNode(iy + 1); } return ret; } inline xt::xtensor FineLayer::nodesFrontRightEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely + 1}); for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { ret(iy) = m_startNode(iy) + m_layer_nelx(iy); } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { ret(iy + 1) = m_startNode(iy + 1) + m_layer_nelx(iy); } return ret; } inline xt::xtensor FineLayer::nodesBackBottomEdge() const { xt::xtensor ret = xt::empty({m_layer_nelx(0) + 1}); for (size_t ix = 0; ix < m_layer_nelx(0) + 1; ++ix) { ret(ix) = m_startNode(0) + ix + (m_layer_nelx(0) + 1) * (m_layer_nelz(0)); } return ret; } inline xt::xtensor FineLayer::nodesBackTopEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelx(nely - 1) + 1}); for (size_t ix = 0; ix < m_layer_nelx(nely - 1) + 1; ++ix) { ret(ix) = m_startNode(nely) + ix + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1)); } return ret; } inline xt::xtensor FineLayer::nodesBackLeftEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely + 1}); for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { ret(iy) = m_startNode(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { ret(iy + 1) = m_startNode(iy + 1) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } return ret; } inline xt::xtensor FineLayer::nodesBackRightEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely + 1}); for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { ret(iy) = m_startNode(iy) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { ret(iy + 1) = m_startNode(iy + 1) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } return ret; } inline xt::xtensor FineLayer::nodesBottomLeftEdge() const { xt::xtensor ret = xt::empty({m_layer_nelz(0) + 1}); for (size_t iz = 0; iz < m_layer_nelz(0) + 1; ++iz) { ret(iz) = m_startNode(0) + iz * (m_layer_nelx(0) + 1); } return ret; } inline xt::xtensor FineLayer::nodesBottomRightEdge() const { xt::xtensor ret = xt::empty({m_layer_nelz(0) + 1}); for (size_t iz = 0; iz < m_layer_nelz(0) + 1; ++iz) { ret(iz) = m_startNode(0) + m_layer_nelx(0) + iz * (m_layer_nelx(0) + 1); } return ret; } inline xt::xtensor FineLayer::nodesTopLeftEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelz(nely - 1) + 1}); for (size_t iz = 0; iz < m_layer_nelz(nely - 1) + 1; ++iz) { ret(iz) = m_startNode(nely) + iz * (m_layer_nelx(nely - 1) + 1); } return ret; } inline xt::xtensor FineLayer::nodesTopRightEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelz(nely - 1) + 1}); for (size_t iz = 0; iz < m_layer_nelz(nely - 1) + 1; ++iz) { ret(iz) = m_startNode(nely) + m_layer_nelx(nely - 1) + iz * (m_layer_nelx(nely - 1) + 1); } return ret; } -inline xt::xtensor FineLayer::nodesBottomFrontEdge() const -{ - return nodesFrontBottomEdge(); -} -inline xt::xtensor FineLayer::nodesBottomBackEdge() const -{ - return nodesBackBottomEdge(); -} -inline xt::xtensor FineLayer::nodesTopFrontEdge() const -{ - return nodesFrontTopEdge(); -} -inline xt::xtensor FineLayer::nodesTopBackEdge() const -{ - return nodesBackTopEdge(); -} -inline xt::xtensor FineLayer::nodesLeftBottomEdge() const -{ - return nodesBottomLeftEdge(); -} -inline xt::xtensor FineLayer::nodesLeftFrontEdge() const -{ - return nodesFrontLeftEdge(); -} -inline xt::xtensor FineLayer::nodesLeftBackEdge() const -{ - return nodesBackLeftEdge(); -} -inline xt::xtensor FineLayer::nodesLeftTopEdge() const -{ - return nodesTopLeftEdge(); -} -inline xt::xtensor FineLayer::nodesRightBottomEdge() const -{ - return nodesBottomRightEdge(); -} -inline xt::xtensor FineLayer::nodesRightTopEdge() const -{ - return nodesTopRightEdge(); -} -inline xt::xtensor FineLayer::nodesRightFrontEdge() const -{ - return nodesFrontRightEdge(); -} -inline xt::xtensor FineLayer::nodesRightBackEdge() const -{ - return nodesBackRightEdge(); -} - inline xt::xtensor FineLayer::nodesFrontBottomOpenEdge() const { xt::xtensor ret = xt::empty({m_layer_nelx(0) - 1}); for (size_t ix = 1; ix < m_layer_nelx(0); ++ix) { ret(ix - 1) = m_startNode(0) + ix; } return ret; } inline xt::xtensor FineLayer::nodesFrontTopOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelx(nely - 1) - 1}); for (size_t ix = 1; ix < m_layer_nelx(nely - 1); ++ix) { ret(ix - 1) = m_startNode(nely) + ix; } return ret; } inline xt::xtensor FineLayer::nodesFrontLeftOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely - 1}); for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { ret(iy - 1) = m_startNode(iy); } for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { ret(iy) = m_startNode(iy + 1); } return ret; } inline xt::xtensor FineLayer::nodesFrontRightOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely - 1}); for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { ret(iy - 1) = m_startNode(iy) + m_layer_nelx(iy); } for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { ret(iy) = m_startNode(iy + 1) + m_layer_nelx(iy); } return ret; } inline xt::xtensor FineLayer::nodesBackBottomOpenEdge() const { xt::xtensor ret = xt::empty({m_layer_nelx(0) - 1}); for (size_t ix = 1; ix < m_layer_nelx(0); ++ix) { ret(ix - 1) = m_startNode(0) + ix + (m_layer_nelx(0) + 1) * (m_layer_nelz(0)); } return ret; } inline xt::xtensor FineLayer::nodesBackTopOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelx(nely - 1) - 1}); for (size_t ix = 1; ix < m_layer_nelx(nely - 1); ++ix) { ret(ix - 1) = m_startNode(nely) + ix + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1)); } return ret; } inline xt::xtensor FineLayer::nodesBackLeftOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely - 1}); for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { ret(iy - 1) = m_startNode(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { ret(iy) = m_startNode(iy + 1) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } return ret; } inline xt::xtensor FineLayer::nodesBackRightOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({nely - 1}); for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) { ret(iy - 1) = m_startNode(iy) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) { ret(iy) = m_startNode(iy + 1) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy)); } return ret; } inline xt::xtensor FineLayer::nodesBottomLeftOpenEdge() const { xt::xtensor ret = xt::empty({m_layer_nelz(0) - 1}); for (size_t iz = 1; iz < m_layer_nelz(0); ++iz) { ret(iz - 1) = m_startNode(0) + iz * (m_layer_nelx(0) + 1); } return ret; } inline xt::xtensor FineLayer::nodesBottomRightOpenEdge() const { xt::xtensor ret = xt::empty({m_layer_nelz(0) - 1}); for (size_t iz = 1; iz < m_layer_nelz(0); ++iz) { ret(iz - 1) = m_startNode(0) + m_layer_nelx(0) + iz * (m_layer_nelx(0) + 1); } return ret; } inline xt::xtensor FineLayer::nodesTopLeftOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelz(nely - 1) - 1}); for (size_t iz = 1; iz < m_layer_nelz(nely - 1); ++iz) { ret(iz - 1) = m_startNode(nely) + iz * (m_layer_nelx(nely - 1) + 1); } return ret; } inline xt::xtensor FineLayer::nodesTopRightOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_layer_nelz(nely - 1) - 1}); for (size_t iz = 1; iz < m_layer_nelz(nely - 1); ++iz) { ret(iz - 1) = m_startNode(nely) + m_layer_nelx(nely - 1) + iz * (m_layer_nelx(nely - 1) + 1); } return ret; } -inline xt::xtensor FineLayer::nodesBottomFrontOpenEdge() const -{ - return nodesFrontBottomOpenEdge(); -} -inline xt::xtensor FineLayer::nodesBottomBackOpenEdge() const -{ - return nodesBackBottomOpenEdge(); -} -inline xt::xtensor FineLayer::nodesTopFrontOpenEdge() const -{ - return nodesFrontTopOpenEdge(); -} -inline xt::xtensor FineLayer::nodesTopBackOpenEdge() const -{ - return nodesBackTopOpenEdge(); -} -inline xt::xtensor FineLayer::nodesLeftBottomOpenEdge() const -{ - return nodesBottomLeftOpenEdge(); -} -inline xt::xtensor FineLayer::nodesLeftFrontOpenEdge() const -{ - return nodesFrontLeftOpenEdge(); -} -inline xt::xtensor FineLayer::nodesLeftBackOpenEdge() const -{ - return nodesBackLeftOpenEdge(); -} -inline xt::xtensor FineLayer::nodesLeftTopOpenEdge() const -{ - return nodesTopLeftOpenEdge(); -} -inline xt::xtensor FineLayer::nodesRightBottomOpenEdge() const -{ - return nodesBottomRightOpenEdge(); -} -inline xt::xtensor FineLayer::nodesRightTopOpenEdge() const -{ - return nodesTopRightOpenEdge(); -} -inline xt::xtensor FineLayer::nodesRightFrontOpenEdge() const -{ - return nodesFrontRightOpenEdge(); -} -inline xt::xtensor FineLayer::nodesRightBackOpenEdge() const -{ - return nodesBackRightOpenEdge(); -} - inline size_t FineLayer::nodesFrontBottomLeftCorner() const { return m_startNode(0); } inline size_t FineLayer::nodesFrontBottomRightCorner() const { return m_startNode(0) + m_layer_nelx(0); } inline size_t FineLayer::nodesFrontTopLeftCorner() const { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely); } inline size_t FineLayer::nodesFrontTopRightCorner() const { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + m_layer_nelx(nely - 1); } inline size_t FineLayer::nodesBackBottomLeftCorner() const { return m_startNode(0) + (m_layer_nelx(0) + 1) * (m_layer_nelz(0)); } inline size_t FineLayer::nodesBackBottomRightCorner() const { return m_startNode(0) + m_layer_nelx(0) + (m_layer_nelx(0) + 1) * (m_layer_nelz(0)); } inline size_t FineLayer::nodesBackTopLeftCorner() const { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1)); } inline size_t FineLayer::nodesBackTopRightCorner() const { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + m_layer_nelx(nely - 1) + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1)); } -inline size_t FineLayer::nodesFrontLeftBottomCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t FineLayer::nodesBottomFrontLeftCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t FineLayer::nodesBottomLeftFrontCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t FineLayer::nodesLeftFrontBottomCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t FineLayer::nodesLeftBottomFrontCorner() const -{ - return nodesFrontBottomLeftCorner(); -} -inline size_t FineLayer::nodesFrontRightBottomCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t FineLayer::nodesBottomFrontRightCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t FineLayer::nodesBottomRightFrontCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t FineLayer::nodesRightFrontBottomCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t FineLayer::nodesRightBottomFrontCorner() const -{ - return nodesFrontBottomRightCorner(); -} -inline size_t FineLayer::nodesFrontLeftTopCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t FineLayer::nodesTopFrontLeftCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t FineLayer::nodesTopLeftFrontCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t FineLayer::nodesLeftFrontTopCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t FineLayer::nodesLeftTopFrontCorner() const -{ - return nodesFrontTopLeftCorner(); -} -inline size_t FineLayer::nodesFrontRightTopCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t FineLayer::nodesTopFrontRightCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t FineLayer::nodesTopRightFrontCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t FineLayer::nodesRightFrontTopCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t FineLayer::nodesRightTopFrontCorner() const -{ - return nodesFrontTopRightCorner(); -} -inline size_t FineLayer::nodesBackLeftBottomCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t FineLayer::nodesBottomBackLeftCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t FineLayer::nodesBottomLeftBackCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t FineLayer::nodesLeftBackBottomCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t FineLayer::nodesLeftBottomBackCorner() const -{ - return nodesBackBottomLeftCorner(); -} -inline size_t FineLayer::nodesBackRightBottomCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t FineLayer::nodesBottomBackRightCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t FineLayer::nodesBottomRightBackCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t FineLayer::nodesRightBackBottomCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t FineLayer::nodesRightBottomBackCorner() const -{ - return nodesBackBottomRightCorner(); -} -inline size_t FineLayer::nodesBackLeftTopCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t FineLayer::nodesTopBackLeftCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t FineLayer::nodesTopLeftBackCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t FineLayer::nodesLeftBackTopCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t FineLayer::nodesLeftTopBackCorner() const -{ - return nodesBackTopLeftCorner(); -} -inline size_t FineLayer::nodesBackRightTopCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t FineLayer::nodesTopBackRightCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t FineLayer::nodesTopRightBackCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t FineLayer::nodesRightBackTopCorner() const -{ - return nodesBackTopRightCorner(); -} -inline size_t FineLayer::nodesRightTopBackCorner() const -{ - return nodesBackTopRightCorner(); -} - -inline xt::xtensor FineLayer::nodesPeriodic() const -{ - xt::xtensor fro = nodesFrontFace(); - xt::xtensor bck = nodesBackFace(); - xt::xtensor lft = nodesLeftFace(); - xt::xtensor rgt = nodesRightFace(); - xt::xtensor bot = nodesBottomFace(); - xt::xtensor top = nodesTopFace(); - - xt::xtensor froBot = nodesFrontBottomOpenEdge(); - xt::xtensor froTop = nodesFrontTopOpenEdge(); - xt::xtensor froLft = nodesFrontLeftOpenEdge(); - xt::xtensor froRgt = nodesFrontRightOpenEdge(); - xt::xtensor bckBot = nodesBackBottomOpenEdge(); - xt::xtensor bckTop = nodesBackTopOpenEdge(); - xt::xtensor bckLft = nodesBackLeftOpenEdge(); - xt::xtensor bckRgt = nodesBackRightOpenEdge(); - xt::xtensor botLft = nodesBottomLeftOpenEdge(); - xt::xtensor botRgt = nodesBottomRightOpenEdge(); - xt::xtensor topLft = nodesTopLeftOpenEdge(); - xt::xtensor topRgt = nodesTopRightOpenEdge(); - - size_t tface = fro.size() + lft.size() + bot.size(); - size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size(); - size_t tnode = 7; - xt::xtensor ret = xt::empty({tface + tedge + tnode, std::size_t(2)}); - - size_t i = 0; - - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesFrontBottomRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackBottomRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackBottomLeftCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesFrontTopLeftCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesFrontTopRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackTopRightCorner(); - ++i; - ret(i, 0) = nodesFrontBottomLeftCorner(); - ret(i, 1) = nodesBackTopLeftCorner(); - ++i; - - for (size_t j = 0; j < froBot.size(); ++j) { - ret(i, 0) = froBot(j); - ret(i, 1) = bckBot(j); - ++i; - } - for (size_t j = 0; j < froBot.size(); ++j) { - ret(i, 0) = froBot(j); - ret(i, 1) = bckTop(j); - ++i; - } - for (size_t j = 0; j < froBot.size(); ++j) { - ret(i, 0) = froBot(j); - ret(i, 1) = froTop(j); - ++i; - } - for (size_t j = 0; j < botLft.size(); ++j) { - ret(i, 0) = botLft(j); - ret(i, 1) = botRgt(j); - ++i; - } - for (size_t j = 0; j < botLft.size(); ++j) { - ret(i, 0) = botLft(j); - ret(i, 1) = topRgt(j); - ++i; - } - for (size_t j = 0; j < botLft.size(); ++j) { - ret(i, 0) = botLft(j); - ret(i, 1) = topLft(j); - ++i; - } - for (size_t j = 0; j < froLft.size(); ++j) { - ret(i, 0) = froLft(j); - ret(i, 1) = froRgt(j); - ++i; - } - for (size_t j = 0; j < froLft.size(); ++j) { - ret(i, 0) = froLft(j); - ret(i, 1) = bckRgt(j); - ++i; - } - for (size_t j = 0; j < froLft.size(); ++j) { - ret(i, 0) = froLft(j); - ret(i, 1) = bckLft(j); - ++i; - } - - for (size_t j = 0; j < fro.size(); ++j) { - ret(i, 0) = fro(j); - ret(i, 1) = bck(j); - ++i; - } - for (size_t j = 0; j < lft.size(); ++j) { - ret(i, 0) = lft(j); - ret(i, 1) = rgt(j); - ++i; - } - for (size_t j = 0; j < bot.size(); ++j) { - ret(i, 0) = bot(j); - ret(i, 1) = top(j); - ++i; - } - - return ret; -} - -inline size_t FineLayer::nodesOrigin() const -{ - return nodesFrontBottomLeftCorner(); -} - -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(); - - for (size_t i = 0; i < nodePer.shape(0); ++i) { - for (size_t j = 0; j < m_ndim; ++j) { - ret(nodePer(i, 1), j) = ret(nodePer(i, 0), j); - } - } - - return GooseFEM::Mesh::renumber(ret); -} - } // namespace Hex8 } // namespace Mesh } // namespace GooseFEM #endif