diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h index 5400c8e..6d0be18 100644 --- a/include/GooseFEM/Mesh.h +++ b/include/GooseFEM/Mesh.h @@ -1,1634 +1,1632 @@ /** 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 { /** Generic mesh operations, and simple mesh definitions. */ 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. +\param coor Nodal coordinates [nnode, ndim]. +\param conn Connectivity [nelem, nne]. \return ElementType(). */ -inline ElementType defaultElementType( - const xt::xtensor& coor, - const xt::xtensor& conn); +template +inline ElementType defaultElementType(const S& coor, const T& conn); /** Abstract 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; virtual ~RegularBase2d() = default; /** Number of elements. \return unsigned int */ size_t nelem() const; /** Number of nodes. \return unsigned int */ size_t nnode() const; /** Number of nodes-per-element == 4. \return unsigned int */ size_t nne() const; /** Number of dimensions == 2. \return unsigned int */ size_t ndim() const; /** Number of elements in x-direction == width of the mesh in units of #h. \return unsigned int */ virtual size_t nelx() const; /** Number of elements in y-direction == height of the mesh, in units of #h, \return unsigned int */ virtual size_t nely() const; /** Linear edge size of one 'block'. \return double */ double h() const; /** The ElementType(). \return element type */ virtual ElementType getElementType() const = 0; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ virtual xt::xtensor coor() const = 0; /** Connectivity [#nelem, #nne]. \return nodes per element */ virtual xt::xtensor conn() const = 0; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ 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 = 0; /** Nodes along the top edge (y = #nely * #h), in order of increasing x. \return List of node numbers. */ virtual xt::xtensor nodesTopEdge() const = 0; /** Nodes along the left edge (x = 0), in order of increasing y. \return List of node numbers. */ virtual xt::xtensor nodesLeftEdge() const = 0; /** Nodes along the right edge (x = #nelx * #h), in order of increasing y. \return List of node numbers. */ virtual xt::xtensor nodesRightEdge() const = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** Alias of nodesBottomLeftCorner(). */ size_t nodesLeftBottomCorner() const; /** Alias of nodesTopLeftCorner(). */ size_t nodesLeftTopCorner() const; /** Alias of nodesBottomRightCorner(). */ size_t nodesRightBottomCorner() const; /** 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). - 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() }; /** Abstract 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; virtual ~RegularBase3d() = default; /** Number of elements. \return unsigned int */ size_t nelem() const; /** Number of nodes. \return unsigned int */ size_t nnode() const; /** Number of nodes-per-element == 4. \return unsigned int */ size_t nne() const; /** Number of dimensions == 2. \return unsigned int */ size_t ndim() const; /** Number of elements in x-direction == width of the mesh in units of #h. \return unsigned int */ virtual size_t nelx() const; /** Number of elements in y-direction == height of the mesh, in units of #h, \return unsigned int */ virtual size_t nely() const; /** Number of elements in y-direction == height of the mesh, in units of #h, \return unsigned int */ virtual size_t nelz() const; /** Linear edge size of one 'block'. \return double */ double h() const; /** The ElementType(). \return element type */ virtual ElementType getElementType() const = 0; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ virtual xt::xtensor coor() const = 0; /** Connectivity [#nelem, #nne]. \return nodes per element */ virtual xt::xtensor conn() const = 0; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ xt::xtensor dofs() const; /** Nodes along the bottom face (y = 0). \return List of node numbers. */ virtual xt::xtensor nodesBottom() const = 0; /** Nodes along the top face (y = #nely * #h). \return List of node numbers. */ virtual xt::xtensor nodesTop() const = 0; /** Nodes along the left face (x = 0). \return List of node numbers. */ virtual xt::xtensor nodesLeft() const = 0; /** Nodes along the right face (x = #nelx * #h). \return List of node numbers. */ virtual xt::xtensor nodesRight() const = 0; /** Nodes along the front face (z = 0). \return List of node numbers. */ virtual xt::xtensor nodesFront() const = 0; /** Nodes along the back face (z = #nelz * #h). \return List of node numbers. */ virtual xt::xtensor nodesBack() const = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** Alias of nodesFrontBottomEdge() */ xt::xtensor nodesBottomFrontEdge() const; /** Alias of nodesBackBottomEdge() */ xt::xtensor nodesBottomBackEdge() const; /** Alias of nodesFrontTopEdge() */ xt::xtensor nodesTopFrontEdge() const; /** Alias of nodesBackTopEdge() */ xt::xtensor nodesTopBackEdge() const; /** Alias of nodesBottomLeftEdge() */ xt::xtensor nodesLeftBottomEdge() const; /** Alias of nodesFrontLeftEdge() */ xt::xtensor nodesLeftFrontEdge() const; /** Alias of nodesBackLeftEdge() */ xt::xtensor nodesLeftBackEdge() const; /** Alias of nodesTopLeftEdge() */ xt::xtensor nodesLeftTopEdge() const; /** Alias of nodesBottomRightEdge() */ xt::xtensor nodesRightBottomEdge() const; /** Alias of nodesTopRightEdge() */ xt::xtensor nodesRightTopEdge() const; /** Alias of nodesFrontRightEdge() */ xt::xtensor nodesRightFrontEdge() const; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** 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 = 0; /** Same as nodesFrontBottomEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesFrontBottomOpenEdge() const = 0; /** Same as nodesFrontTopEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesFrontTopOpenEdge() const = 0; /** Same as nodesFrontLeftEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesFrontLeftOpenEdge() const = 0; /** Same as nodesFrontRightEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesFrontRightOpenEdge() const = 0; /** Same as nodesBackBottomEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesBackBottomOpenEdge() const = 0; /** Same as nodesBackTopEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesBackTopOpenEdge() const = 0; /** Same as nodesBackLeftEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesBackLeftOpenEdge() const = 0; /** Same as nodesBackRightEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesBackRightOpenEdge() const = 0; /** Same as nodesBottomLeftEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesBottomLeftOpenEdge() const = 0; /** Same as nodesBottomRightEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesBottomRightOpenEdge() const = 0; /** Same as nodesTopLeftEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesTopLeftOpenEdge() const = 0; /** Same as nodesTopRightEdge() but without corners. \return List of node numbers. */ virtual xt::xtensor nodesTopRightOpenEdge() const = 0; /** Alias of nodesFrontBottomOpenEdge(). */ xt::xtensor nodesBottomFrontOpenEdge() const; /** Alias of nodesBackBottomOpenEdge(). */ xt::xtensor nodesBottomBackOpenEdge() const; /** Alias of nodesFrontTopOpenEdge(). */ xt::xtensor nodesTopFrontOpenEdge() const; /** Alias of nodesBackTopOpenEdge(). */ xt::xtensor nodesTopBackOpenEdge() const; /** Alias of nodesBottomLeftOpenEdge(). */ xt::xtensor nodesLeftBottomOpenEdge() const; /** Alias of nodesFrontLeftOpenEdge(). */ xt::xtensor nodesLeftFrontOpenEdge() const; /** Alias of nodesBackLeftOpenEdge(). */ xt::xtensor nodesLeftBackOpenEdge() const; /** Alias of nodesTopLeftOpenEdge(). */ xt::xtensor nodesLeftTopOpenEdge() const; /** Alias of nodesBottomRightOpenEdge(). */ xt::xtensor nodesRightBottomOpenEdge() const; /** Alias of nodesTopRightOpenEdge(). */ xt::xtensor nodesRightTopOpenEdge() const; /** Alias of nodesFrontRightOpenEdge(). */ xt::xtensor nodesRightFrontOpenEdge() const; /** Alias of nodesBackRightOpenEdge(). */ xt::xtensor nodesRightBackOpenEdge() const; /** Front-Bottom-Left corner node. \return Node number. */ virtual size_t nodesFrontBottomLeftCorner() const = 0; /** Front-Bottom-Right corner node. \return Node number. */ virtual size_t nodesFrontBottomRightCorner() const = 0; /** Front-Top-Left corner node. \return Node number. */ virtual size_t nodesFrontTopLeftCorner() const = 0; /** Front-Top-Right corner node. \return Node number. */ virtual size_t nodesFrontTopRightCorner() const = 0; /** Back-Bottom-Left corner node. \return Node number. */ virtual size_t nodesBackBottomLeftCorner() const = 0; /** Back-Bottom-Right corner node. \return Node number. */ virtual size_t nodesBackBottomRightCorner() const = 0; /** Back-Top-Left corner node. \return Node number. */ virtual size_t nodesBackTopLeftCorner() const = 0; /** Back-Top-Right corner node. \return Node number. */ virtual size_t nodesBackTopRightCorner() const = 0; /** Alias of nodesFrontBottomLeftCorner(). */ size_t nodesFrontLeftBottomCorner() const; /** Alias of nodesFrontBottomLeftCorner(). */ size_t nodesBottomFrontLeftCorner() const; /** Alias of nodesFrontBottomLeftCorner(). */ size_t nodesBottomLeftFrontCorner() const; /** Alias of nodesFrontBottomLeftCorner(). */ size_t nodesLeftFrontBottomCorner() const; /** Alias of nodesFrontBottomLeftCorner(). */ size_t nodesLeftBottomFrontCorner() const; /** Alias of nodesFrontBottomRightCorner(). */ size_t nodesFrontRightBottomCorner() const; /** Alias of nodesFrontBottomRightCorner(). */ size_t nodesBottomFrontRightCorner() const; /** Alias of nodesFrontBottomRightCorner(). */ size_t nodesBottomRightFrontCorner() const; /** Alias of nodesFrontBottomRightCorner(). */ size_t nodesRightFrontBottomCorner() const; /** Alias of nodesFrontBottomRightCorner(). */ size_t nodesRightBottomFrontCorner() const; /** Alias of nodesFrontTopLeftCorner(). */ size_t nodesFrontLeftTopCorner() const; /** Alias of nodesFrontTopLeftCorner(). */ size_t nodesTopFrontLeftCorner() const; /** Alias of nodesFrontTopLeftCorner(). */ size_t nodesTopLeftFrontCorner() const; /** Alias of nodesFrontTopLeftCorner(). */ size_t nodesLeftFrontTopCorner() const; /** Alias of nodesFrontTopLeftCorner(). */ size_t nodesLeftTopFrontCorner() const; /** Alias of nodesFrontTopRightCorner(). */ size_t nodesFrontRightTopCorner() const; /** Alias of nodesFrontTopRightCorner(). */ size_t nodesTopFrontRightCorner() const; /** Alias of nodesFrontTopRightCorner(). */ size_t nodesTopRightFrontCorner() const; /** Alias of nodesFrontTopRightCorner(). */ size_t nodesRightFrontTopCorner() const; /** Alias of nodesFrontTopRightCorner(). */ size_t nodesRightTopFrontCorner() const; /** Alias of nodesBackBottomLeftCorner(). */ size_t nodesBackLeftBottomCorner() const; /** Alias of nodesBackBottomLeftCorner(). */ size_t nodesBottomBackLeftCorner() const; /** Alias of nodesBackBottomLeftCorner(). */ size_t nodesBottomLeftBackCorner() const; /** Alias of nodesBackBottomLeftCorner(). */ size_t nodesLeftBackBottomCorner() const; /** Alias of nodesBackBottomLeftCorner(). */ size_t nodesLeftBottomBackCorner() const; /** Alias of nodesBackBottomRightCorner(). */ size_t nodesBackRightBottomCorner() const; /** Alias of nodesBackBottomRightCorner(). */ size_t nodesBottomBackRightCorner() const; /** Alias of nodesBackBottomRightCorner(). */ size_t nodesBottomRightBackCorner() const; /** Alias of nodesBackBottomRightCorner(). */ size_t nodesRightBackBottomCorner() const; /** Alias of nodesBackBottomRightCorner(). */ size_t nodesRightBottomBackCorner() const; /** Alias of nodesBackTopLeftCorner(). */ size_t nodesBackLeftTopCorner() const; /** Alias of nodesBackTopLeftCorner(). */ size_t nodesTopBackLeftCorner() const; /** Alias of nodesBackTopLeftCorner(). */ size_t nodesTopLeftBackCorner() const; /** Alias of nodesBackTopLeftCorner(). */ size_t nodesLeftBackTopCorner() const; /** Alias of nodesBackTopLeftCorner(). */ size_t nodesLeftTopBackCorner() const; /** Alias of nodesBackTopRightCorner(). */ size_t nodesBackRightTopCorner() const; /** Alias of nodesBackTopRightCorner(). */ size_t nodesTopBackRightCorner() const; /** Alias of nodesBackTopRightCorner(). */ size_t nodesTopRightBackCorner() const; /** Alias of nodesBackTopRightCorner(). */ size_t nodesRightBackTopCorner() const; /** 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 coor_a Nodal coordinates of mesh "a" [nnode, ndim]. +\param coor_b Nodal coordinates of mesh "b" [nnode, ndim]. \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. \return Overlapping nodes. */ +template inline xt::xtensor overlapping( - const xt::xtensor& coor_a, - const xt::xtensor& coor_b, + const S& coor_a, + const T& 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 coor_a Nodal coordinates of mesh "a" [nnode, ndim]. + \param conn_a Connectivity of mesh "a" [nelem, nne]. + \param overlapping_nodes_a Node-numbers of mesh "a" that overlap with mesh "b" [n]. + \param coor_b Nodal coordinates of mesh "b" [nnode, ndim]. + \param conn_b Connectivity of mesh "b" [nelem, nne]. + \param overlapping_nodes_b Node-numbers of mesh "b" that overlap with mesh "a" [n]. \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. */ + template 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, + const C& coor_a, + const E& conn_a, + const N& overlapping_nodes_a, + const C& coor_b, + const E& conn_b, + const N& overlapping_nodes_b, bool check_position = true, double rtol = 1e-5, double atol = 1e-8); /** Number of sub meshes == 2. \return unsigned int */ size_t nmesh() const; /** Number of elements. \return unsigned int */ size_t nelem() const; /** Number of nodes. \return unsigned int */ size_t nnode() const; /** Number of nodes-per-element. \return unsigned int */ size_t nne() const; /** Number of dimensions. \return unsigned int */ size_t ndim() const; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ xt::xtensor coor() const; /** Connectivity [#nelem, #nne]. \return nodes per element */ xt::xtensor conn() const; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ xt::xtensor dofs() const; /** Node-map per sub-mesh. \return nodes per mesh */ std::vector> nodemap() const; /** Element-map per sub-mesh. \return elements per 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; + template + T nodeset(const T& 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; + template + T elemset(const T& 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. + \param coor Nodal coordinates [nnode, ndim]. + \param conn Connectivity [nelem, nne]. */ - void push_back(const xt::xtensor& coor, const xt::xtensor& conn); + template + void push_back(const C& coor, const E& conn); /** Number of sub meshes. \return unsigned int */ size_t nmesh() const; /** Number of elements. \return unsigned int */ size_t nelem() const; /** Number of nodes. \return unsigned int */ size_t nnode() const; /** Number of nodes-per-element. \return unsigned int */ size_t nne() const; /** Number of dimensions. \return unsigned int */ size_t ndim() const; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ xt::xtensor coor() const; /** Connectivity [#nelem, #nne]. \return nodes per element */ xt::xtensor conn() const; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ xt::xtensor dofs() const; /** Node-map per sub-mesh. \return nodes per mesh */ std::vector> nodemap() const; /** Element-map per sub-mesh. \return elements per 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; + template + T nodeset(const T& 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; + template + T elemset(const T& 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; + template + T nodeset(const std::vector& set) const; + + /** \copydoc nodeset(const std::vector&) const */ + template + T nodeset(std::initializer_list 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; + template + T elemset(const std::vector& set) const; + + /** \copydoc elemset(const std::vector&) const */ + template + T elemset(std::initializer_list set) const; protected: xt::xtensor m_coor; ///< Nodal coordinates [#nnode, #ndim] xt::xtensor m_conn; ///< Connectivity [#nelem, #nne] std::vector> m_map; ///< See nodemap(size_t) std::vector m_nel; ///< Number of elements per sub-mesh. std::vector m_el_offset; ///< First element of every sub-mesh. double m_rtol; ///< Relative tolerance to find overlapping nodes. double m_atol; ///< Absolute tolerance to find overlapping nodes. }; /** Vertically stack meshes. */ class Vstack : public Stitch { public: /** \param check_overlap Check if nodes are overlapping when adding a mesh. \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. */ Vstack(bool check_overlap = true, double rtol = 1e-5, double atol = 1e-8); /** Add a mesh to the top of the current stack. Each time the current `nodes_bottom` are stitched with the then highest `nodes_top`. - \param coor Nodal coordinates. - \param conn Connectivity. - \param nodes_bottom Nodes along the bottom edge. - \param nodes_top Nodes along the top edge. + \param coor Nodal coordinates [nnode, ndim]. + \param conn Connectivity [nelem, nne]. + \param nodes_bottom Nodes along the bottom edge [n]. + \param nodes_top Nodes along the top edge [n]. */ - void push_back( - const xt::xtensor& coor, - const xt::xtensor& conn, - const xt::xtensor& nodes_bottom, - const xt::xtensor& nodes_top); + template + void push_back(const C& coor, const E& conn, const N& nodes_bottom, const N& nodes_top); private: std::vector> m_nodes_bot; ///< Bottom nodes of each mesh (renumbered). std::vector> m_nodes_top; ///< Top nodes of each mesh (renumbered). bool m_check_overlap; ///< Check if nodes are overlapping when adding a mesh. }; /** Renumber indices to lowest possible index. For example: \f$ \begin{bmatrix} 0 & 1 \\ 5 & 4 \end{bmatrix} \f$ is renumbered to \f$ \begin{bmatrix} 0 & 1 \\ 3 & 2 \end{bmatrix} \f$ Or, in pseudo-code, the result of this function is that: dofs = renumber(dofs) sort(unique(dofs[:])) == range(max(dofs+1)) \note One can use the wrapper function renumber(). This class gives more advanced features. */ 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. +\param dofs DOF-numbers [nnode, ndim]. \return Renumbered DOF-numbers. */ -inline xt::xtensor renumber(const xt::xtensor& dofs); +template +inline T renumber(const T& 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; + template + Reorder(const std::initializer_list args); /** 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 \f$ \begin{bmatrix} 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \f$ \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. +\param conn Connectivity [nelem, nne]. \return Coordination per node. */ -inline xt::xtensor coordination(const xt::xtensor& conn); +template +inline xt::xtensor coordination(const E& 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 41ca8ea..e1f26e9 100644 --- a/include/GooseFEM/Mesh.hpp +++ b/include/GooseFEM/Mesh.hpp @@ -1,1291 +1,1311 @@ /** 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) +template +inline ElementType defaultElementType(const S& coor, const T& conn) { + GOOSEFEM_ASSERT(coor.dimension() == 2); + GOOSEFEM_ASSERT(conn.dimension() == 2); + 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 size_t RegularBase2d::nodesLeftBottomCorner() const { return this->nodesBottomLeftCorner(); } inline size_t RegularBase2d::nodesLeftTopCorner() const { return this->nodesTopLeftCorner(); } inline size_t RegularBase2d::nodesRightBottomCorner() const { return this->nodesBottomRightCorner(); } inline size_t RegularBase2d::nodesRightTopCorner() const { 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 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::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::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 +template inline xt::xtensor overlapping( - const xt::xtensor& coor_a, - const xt::xtensor& coor_b, + const S& coor_a, + const T& coor_b, double rtol, double atol) { + GOOSEFEM_ASSERT(coor_a.dimension() == 2); + GOOSEFEM_ASSERT(coor_b.dimension() == 2); 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; } +template 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, + const C& coor_a, + const E& conn_a, + const N& overlapping_nodes_a, + const C& coor_b, + const E& conn_b, + const N& overlapping_nodes_b, bool check_position, double rtol, double atol) { UNUSED(rtol); UNUSED(atol); + GOOSEFEM_ASSERT(coor_a.dimension() == 2); + GOOSEFEM_ASSERT(conn_a.dimension() == 2); + GOOSEFEM_ASSERT(overlapping_nodes_a.dimension() == 1); + GOOSEFEM_ASSERT(coor_b.dimension() == 2); + GOOSEFEM_ASSERT(conn_b.dimension() == 2); + GOOSEFEM_ASSERT(overlapping_nodes_b.dimension() == 1); 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_CHECK(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 +template +inline T ManualStitch::nodeset(const T& 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 +template +inline T ManualStitch::elemset(const T& 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) +template +inline void Stitch::push_back(const C& coor, const E& conn) { + GOOSEFEM_ASSERT(coor.dimension() == 2); + GOOSEFEM_ASSERT(conn.dimension() == 2); + 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 stitch( m_coor, m_conn, xt::view(overlap, 0, xt::all()), coor, conn, xt::view(overlap, 1, xt::all()), false); m_coor = stitch.coor(); m_conn = stitch.conn(); m_map.push_back(stitch.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 +template +inline T Stitch::nodeset(const T& 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 +template +inline T Stitch::elemset(const T& 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 +template +inline T 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 +template +inline T 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 T Stitch::nodeset(std::initializer_list set) const +{ + return this->nodeset(std::vector(set)); +} + +template +inline T Stitch::elemset(std::initializer_list set) const +{ + return this->elemset(std::vector(set)); +} + inline Vstack::Vstack(bool check_overlap, double rtol, double atol) { m_check_overlap = check_overlap; m_rtol = rtol; m_atol = atol; } -inline void Vstack::push_back( - const xt::xtensor& coor, - const xt::xtensor& conn, - const xt::xtensor& nodes_bot, - const xt::xtensor& nodes_top) +template +inline void Vstack::push_back(const C& coor, const E& conn, const N& nodes_bot, const N& nodes_top) { 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); m_nodes_bot.push_back(nodes_bot); m_nodes_top.push_back(nodes_top); return; } GOOSEFEM_ASSERT(nodes_bot.size() == m_nodes_top.back().size()); size_t index = m_map.size(); double shift = xt::amax(xt::view(m_coor, xt::all(), 1))(); auto x = coor; xt::view(x, xt::all(), 1) += shift; ManualStitch stitch( m_coor, m_conn, m_nodes_top.back(), x, conn, nodes_bot, m_check_overlap, m_rtol, m_atol); m_nodes_bot.push_back(stitch.nodeset(nodes_bot, 1)); m_nodes_top.push_back(stitch.nodeset(nodes_top, 1)); m_coor = stitch.coor(); m_conn = stitch.conn(); m_map.push_back(stitch.nodemap(1)); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]); } 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) +template +inline T renumber(const T& dofs) { return Renumber(dofs).apply(dofs); } -inline Reorder::Reorder(const std::initializer_list> args) +template +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) +template +inline xt::xtensor coordination(const E& conn) { + GOOSEFEM_ASSERT(conn.dimension() == 2); + 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/python/Allocate.hpp b/python/Allocate.hpp index 0785ddc..3e6a6dd 100644 --- a/python/Allocate.hpp +++ b/python/Allocate.hpp @@ -1,35 +1,35 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include namespace py = pybind11; void init_Allocate(py::module& m) { m.def("AsTensor", static_cast (*)( - const xt::xarray&, const std::vector&)>(&GooseFEM::AsTensor), + const xt::xarray&, const std::vector&)>(&GooseFEM::AsTensor), "See :cpp:func:`GooseFEM::AsTensor`.", py::arg("arg"), py::arg("shape")); m.def("AsTensor", static_cast (*)( - size_t, const xt::xarray&, size_t)>(&GooseFEM::AsTensor), + size_t, const xt::xarray&, size_t)>(&GooseFEM::AsTensor), "See :cpp:func:`GooseFEM::AsTensor`.", py::arg("rank"), py::arg("arg"), py::arg("n")); m.def("as3d", &GooseFEM::as3d>, "See :cpp:func:`GooseFEM::as3d`.", py::arg("arg")); } diff --git a/python/Mesh.hpp b/python/Mesh.hpp index 8a609d9..edafd26 100644 --- a/python/Mesh.hpp +++ b/python/Mesh.hpp @@ -1,433 +1,433 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_Mesh(py::module& m) { py::enum_(m, "ElementType", "ElementType." "See :cpp:enum:`GooseFEM::Mesh::ElementType`.") .value("Unknown", GooseFEM::Mesh::ElementType::Unknown) .value("Tri3", GooseFEM::Mesh::ElementType::Tri3) .value("Quad4", GooseFEM::Mesh::ElementType::Quad4) .value("Hex8", GooseFEM::Mesh::ElementType::Hex8) .export_values(); m.def("overlapping", - &GooseFEM::Mesh::overlapping, + &GooseFEM::Mesh::overlapping, xt::xtensor>, "Find overlapping nodes." "See :cpp:func:`GooseFEM::Mesh::overlapping`.", py::arg("coor_a"), py::arg("coor_b"), py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8); py::class_(m, "ManualStitch") .def(py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&, const xt::xtensor&, const xt::xtensor&, const xt::xtensor&, bool, double, double>(), "Manually stitch meshes." "See :cpp:class:`GooseFEM::Mesh::ManualStitch`.", py::arg("coor_a"), py::arg("conn_a"), py::arg("overlapping_nodes_a"), py::arg("coor_b"), py::arg("conn_b"), py::arg("overlapping_nodes_b"), py::arg("check_position") = true, py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8) .def("nmesh", &GooseFEM::Mesh::ManualStitch::nmesh, "Number of sub meshes." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nmesh`.") .def("nnode", &GooseFEM::Mesh::ManualStitch::nnode, "Number of nodes of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nnode`.") .def("nelem", &GooseFEM::Mesh::ManualStitch::nelem, "Number of elements of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nelem`.") .def("ndim", &GooseFEM::Mesh::ManualStitch::ndim, "Number of dimensions (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::ndim`.") .def("nne", &GooseFEM::Mesh::ManualStitch::nne, "Number of nodes per element (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nne`.") .def("coor", &GooseFEM::Mesh::ManualStitch::coor, "Coordinates of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::coor`.") .def("conn", &GooseFEM::Mesh::ManualStitch::conn, "Connectivity of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::conn`.") .def("dofs", &GooseFEM::Mesh::ManualStitch::dofs, "DOF numbers per node." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::dofs`.") .def("nodemap", py::overload_cast<>( &GooseFEM::Mesh::ManualStitch::nodemap, py::const_), "Node-map for givall sub-meshes." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.") .def("elemmap", py::overload_cast<>( &GooseFEM::Mesh::ManualStitch::elemmap, py::const_), "Element-map for all sub-meshes." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemmap`.") .def("nodemap", py::overload_cast( &GooseFEM::Mesh::ManualStitch::nodemap, py::const_), "Node-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.", py::arg("mesh_index")) .def("elemmap", py::overload_cast( &GooseFEM::Mesh::ManualStitch::elemmap, py::const_), "Element-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemmap`.", py::arg("mesh_index")) .def("nodeset", - &GooseFEM::Mesh::ManualStitch::nodeset, + &GooseFEM::Mesh::ManualStitch::nodeset>, "Convert node-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodeset`.", py::arg("set"), py::arg("mesh_index")) .def("elemset", - &GooseFEM::Mesh::ManualStitch::elemset, + &GooseFEM::Mesh::ManualStitch::elemset>, "Convert element-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemset`.", py::arg("set"), py::arg("mesh_index")) .def("__repr__", [](const GooseFEM::Mesh::ManualStitch&) { return ""; }); py::class_(m, "Stitch") .def(py::init(), "Stitch meshes." "See :cpp:class:`GooseFEM::Mesh::Stitch`.", py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8) .def("push_back", - &GooseFEM::Mesh::Stitch::push_back, + &GooseFEM::Mesh::Stitch::push_back, xt::xtensor>, "Add mesh to be stitched." "See :cpp:func:`GooseFEM::Mesh::Stitch::push_back`.", py::arg("coor"), py::arg("conn")) .def("nmesh", &GooseFEM::Mesh::Stitch::nmesh, "Number of sub meshes." "See :cpp:func:`GooseFEM::Mesh::Stitch::nmesh`.") .def("nnode", &GooseFEM::Mesh::Stitch::nnode, "Number of nodes of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nnode`.") .def("nelem", &GooseFEM::Mesh::Stitch::nelem, "Number of elements of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nelem`.") .def("ndim", &GooseFEM::Mesh::Stitch::ndim, "Number of dimensions (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::Stitch::ndim`.") .def("nne", &GooseFEM::Mesh::Stitch::nne, "Number of nodes per element (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::Stitch::nne`.") .def("coor", &GooseFEM::Mesh::Stitch::coor, "Coordinates of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::coor`.") .def("conn", &GooseFEM::Mesh::Stitch::conn, "Connectivity of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::conn`.") .def("dofs", &GooseFEM::Mesh::Stitch::dofs, "DOF numbers per node." "See :cpp:func:`GooseFEM::Mesh::Stitch::dofs`.") .def("nodemap", py::overload_cast<>( &GooseFEM::Mesh::Stitch::nodemap, py::const_), "Node-map for givall sub-meshes." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodemap`.") .def("elemmap", py::overload_cast<>( &GooseFEM::Mesh::Stitch::elemmap, py::const_), "Element-map for all sub-meshes." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemmap`.") .def("nodemap", py::overload_cast( &GooseFEM::Mesh::Stitch::nodemap, py::const_), "Node-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodemap`.", py::arg("mesh_index")) .def("elemmap", py::overload_cast( &GooseFEM::Mesh::Stitch::elemmap, py::const_), "Element-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemmap`.", py::arg("mesh_index")) .def("nodeset", - py::overload_cast&, size_t>( - &GooseFEM::Mesh::Stitch::nodeset, py::const_), + static_cast (GooseFEM::Mesh::Stitch::*)( + const xt::xtensor&, size_t) const>(&GooseFEM::Mesh::Stitch::nodeset), "Convert node-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.", py::arg("set"), py::arg("mesh_index")) .def("elemset", - py::overload_cast&, size_t>( - &GooseFEM::Mesh::Stitch::elemset, py::const_), + static_cast (GooseFEM::Mesh::Stitch::*)( + const xt::xtensor&, size_t) const>(&GooseFEM::Mesh::Stitch::elemset), "Convert element-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.", py::arg("set"), py::arg("mesh_index")) .def("nodeset", - py::overload_cast>&>( - &GooseFEM::Mesh::Stitch::nodeset, py::const_), + static_cast (GooseFEM::Mesh::Stitch::*)( + const std::vector>&) const>(&GooseFEM::Mesh::Stitch::nodeset), "Convert node-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.", py::arg("set")) .def("elemset", - py::overload_cast>&>( - &GooseFEM::Mesh::Stitch::elemset, py::const_), + static_cast (GooseFEM::Mesh::Stitch::*)( + const std::vector>&) const>(&GooseFEM::Mesh::Stitch::elemset), "Convert element-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.", py::arg("set")) .def("__repr__", [](const GooseFEM::Mesh::Stitch&) { return ""; }); py::class_(m, "Vstack") .def(py::init(), "Vstack meshes." "See :cpp:class:`GooseFEM::Mesh::Vstack`.", py::arg("check_overlap") = true, py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8) .def("push_back", - &GooseFEM::Mesh::Vstack::push_back, + &GooseFEM::Mesh::Vstack::push_back, xt::xtensor, xt::xtensor>, "Add mesh to be stitched." "See :cpp:func:`GooseFEM::Mesh::Vstack::push_back`.", py::arg("coor"), py::arg("conn"), py::arg("nodes_bottom"), py::arg("nodes_top")) .def("__repr__", [](const GooseFEM::Mesh::Vstack&) { return ""; }); py::class_(m, "Renumber") .def(py::init&>(), "Renumber to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Renumber`.", py::arg("dofs")) .def("get", &GooseFEM::Mesh::Renumber::apply>, "Get renumbered DOFs." "See :cpp:func:`GooseFEM::Mesh::Renumber::get`.") .def("apply", &GooseFEM::Mesh::Renumber::apply>, "Get renumbered list." "See :cpp:func:`GooseFEM::Mesh::Renumber::apply`.") .def("index", &GooseFEM::Mesh::Renumber::index, "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``." "See :cpp:func:`GooseFEM::Mesh::Renumber::index`.") .def("__repr__", [](const GooseFEM::Mesh::Renumber&) { return ""; }); py::class_(m, "Reorder") .def(py::init([](xt::xtensor& a) { return new GooseFEM::Mesh::Reorder({a}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def(py::init([](xt::xtensor& a, xt::xtensor& b) { return new GooseFEM::Mesh::Reorder({a, b}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def(py::init( [](xt::xtensor& a, xt::xtensor& b, xt::xtensor& c) { return new GooseFEM::Mesh::Reorder({a, b, c}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def(py::init([](xt::xtensor& a, xt::xtensor& b, xt::xtensor& c, xt::xtensor& d) { return new GooseFEM::Mesh::Reorder({a, b, c, d}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def("get", &GooseFEM::Mesh::Reorder::apply>, "Reorder matrix (e.g. ``dofs``)." "See :cpp:func:`GooseFEM::Mesh::Reorder::get`.", py::arg("dofs")) .def("apply", &GooseFEM::Mesh::Reorder::apply>, "Get reordered list." "See :cpp:func:`GooseFEM::Mesh::Reorder::apply`.") .def("index", &GooseFEM::Mesh::Reorder::index, "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``." "See :cpp:func:`GooseFEM::Mesh::Reorder::index`.") .def("__repr__", [](const GooseFEM::Mesh::Reorder&) { return ""; }); m.def("dofs", - &GooseFEM::Mesh::dofs, - "List with DOF-numbers in sequential order." - "See :cpp:func:`GooseFEM::Mesh::dofs`.", - py::arg("nnode"), - py::arg("ndim")); + &GooseFEM::Mesh::dofs, + "List with DOF-numbers in sequential order." + "See :cpp:func:`GooseFEM::Mesh::dofs`.", + py::arg("nnode"), + py::arg("ndim")); m.def("renumber", - &GooseFEM::Mesh::renumber, - "Renumber to lowest possible indices." - "See :cpp:func:`GooseFEM::Mesh::renumber`.", - py::arg("dofs")); + &GooseFEM::Mesh::renumber>, + "Renumber to lowest possible indices." + "See :cpp:func:`GooseFEM::Mesh::renumber`.", + py::arg("dofs")); m.def("coordination", - &GooseFEM::Mesh::coordination, - "Coordination number of each node." - "See :cpp:func:`GooseFEM::Mesh::coordination`.", - py::arg("conn")); + &GooseFEM::Mesh::coordination>, + "Coordination number of each node." + "See :cpp:func:`GooseFEM::Mesh::coordination`.", + py::arg("conn")); m.def("elem2node", - &GooseFEM::Mesh::elem2node, - "Element-numbers connected to each node." - "See :cpp:func:`GooseFEM::Mesh::elem2node`.", - py::arg("conn"), - py::arg("sorted") = true); + &GooseFEM::Mesh::elem2node, + "Element-numbers connected to each node." + "See :cpp:func:`GooseFEM::Mesh::elem2node`.", + py::arg("conn"), + py::arg("sorted") = true); m.def("edgesize", py::overload_cast&, const xt::xtensor&>( &GooseFEM::Mesh::edgesize), "Get the edge size of all elements." "See :cpp:func:`GooseFEM::Mesh::edgesize`.", py::arg("coor"), py::arg("conn")); m.def("edgesize", py::overload_cast< const xt::xtensor&, const xt::xtensor&, GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize), "Get the edge size of all elements." "See :cpp:func:`GooseFEM::Mesh::edgesize`.", py::arg("coor"), py::arg("conn"), py::arg("type")); m.def("centers", py::overload_cast&, const xt::xtensor&>( &GooseFEM::Mesh::centers), "Coordinates of the center of each element." "See :cpp:func:`GooseFEM::Mesh::centers`.", py::arg("coor"), py::arg("conn")); m.def("centers", py::overload_cast< const xt::xtensor&, const xt::xtensor&, GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::centers), "Coordinates of the center of each element." "See :cpp:func:`GooseFEM::Mesh::centers`.", py::arg("coor"), py::arg("conn"), py::arg("type")); m.def("elemmap2nodemap", py::overload_cast< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(&GooseFEM::Mesh::elemmap2nodemap), "Convert an element-map to a node-map." "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.", py::arg("elem_map"), py::arg("coor"), py::arg("conn")); m.def("elemmap2nodemap", py::overload_cast< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&, GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::elemmap2nodemap), "Convert an element-map to a node-map." "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.", py::arg("elem_map"), py::arg("coor"), py::arg("conn"), py::arg("type")); } diff --git a/test/basic/Mesh.cpp b/test/basic/Mesh.cpp index 1d58a2a..9eda479 100644 --- a/test/basic/Mesh.cpp +++ b/test/basic/Mesh.cpp @@ -1,424 +1,424 @@ #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Mesh", "Mesh.h") { SECTION("overlapping") { GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0); auto coor_a = mesh.coor(); auto overlap_a = mesh.nodesTopEdge(); auto coor_b = mesh.coor(); auto overlap_b = mesh.nodesBottomEdge(); xt::view(coor_b, xt::all(), 1) += 5.0; auto overlap = GooseFEM::Mesh::overlapping(coor_a, coor_b); REQUIRE(xt::all(xt::equal(xt::view(overlap, 0, xt::all()), overlap_a))); REQUIRE(xt::all(xt::equal(xt::view(overlap, 1, xt::all()), overlap_b))); } SECTION("ManualStitch") { GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0); auto coor_a = mesh.coor(); auto conn_a = mesh.conn(); auto overlap_a = mesh.nodesTopEdge(); auto coor_b = mesh.coor(); auto conn_b = mesh.conn(); auto overlap_b = mesh.nodesBottomEdge(); xt::view(coor_b, xt::all(), 1) += 5.0; GooseFEM::Mesh::ManualStitch stitch(coor_a, conn_a, overlap_a, coor_b, conn_b, overlap_b); GooseFEM::Mesh::Quad4::Regular res(5, 10, 1.0); REQUIRE(xt::allclose(stitch.coor(), res.coor())); REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn()))); REQUIRE(stitch.nodemap().size() == 2); REQUIRE(stitch.elemmap().size() == 2); REQUIRE(xt::all(xt::equal(stitch.nodemap(0), xt::arange(6 * 6)))); REQUIRE(xt::all(xt::equal(stitch.nodemap(1), xt::arange(6 * 6) + 5 * 6))); REQUIRE(xt::all(xt::equal(stitch.elemmap(0), xt::arange(5 * 5)))); REQUIRE(xt::all(xt::equal(stitch.elemmap(1), xt::arange(5 * 5) + 5 * 5))); - REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange(6 * 6), 0), xt::arange(6 * 6)))); - REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange(6 * 6), 1), xt::arange(6 * 6) + 5 * 6))); - REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange(5 * 5), 0), xt::arange(5 * 5)))); - REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange(5 * 5), 1), xt::arange(5 * 5) + 5 * 5))); + REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::eval(xt::arange(6 * 6)), 0), xt::arange(6 * 6)))); + REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::eval(xt::arange(6 * 6)), 1), xt::arange(6 * 6) + 5 * 6))); + REQUIRE(xt::all(xt::equal(stitch.elemset(xt::eval(xt::arange(5 * 5)), 0), xt::arange(5 * 5)))); + REQUIRE(xt::all(xt::equal(stitch.elemset(xt::eval(xt::arange(5 * 5)), 1), xt::arange(5 * 5) + 5 * 5))); REQUIRE(xt::all(xt::equal(stitch.nodeset(overlap_a, 0), stitch.nodeset(overlap_b, 1)))); } SECTION("Stitch") { GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0); auto coor_a = mesh.coor(); auto conn_a = mesh.conn(); auto overlap_a = mesh.nodesTopEdge(); auto nset = mesh.nodesLeftEdge(); auto eset = xt::eval(xt::arange(mesh.nelem())); auto coor_b = mesh.coor(); auto conn_b = mesh.conn(); auto overlap_b = mesh.nodesBottomEdge(); xt::view(coor_b, xt::all(), 1) += 5.0; GooseFEM::Mesh::Stitch stitch; stitch.push_back(coor_a, conn_a); stitch.push_back(coor_b, conn_b); GooseFEM::Mesh::Quad4::Regular res(5, 10, 1.0); REQUIRE(xt::allclose(stitch.coor(), res.coor())); REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn()))); REQUIRE(stitch.nodemap().size() == 2); REQUIRE(stitch.elemmap().size() == 2); REQUIRE(xt::all(xt::equal(stitch.nodemap(0), xt::arange(6 * 6)))); REQUIRE(xt::all(xt::equal(stitch.nodemap(1), xt::arange(6 * 6) + 5 * 6))); REQUIRE(xt::all(xt::equal(stitch.elemmap(0), xt::arange(5 * 5)))); REQUIRE(xt::all(xt::equal(stitch.elemmap(1), xt::arange(5 * 5) + 5 * 5))); - REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange(6 * 6), 0), xt::arange(6 * 6)))); - REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange(6 * 6), 1), xt::arange(6 * 6) + 5 * 6))); - REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange(5 * 5), 0), xt::arange(5 * 5)))); - REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange(5 * 5), 1), xt::arange(5 * 5) + 5 * 5))); + REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::eval(xt::arange(6 * 6)), 0), xt::arange(6 * 6)))); + REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::eval(xt::arange(6 * 6)), 1), xt::arange(6 * 6) + 5 * 6))); + REQUIRE(xt::all(xt::equal(stitch.elemset(xt::eval(xt::arange(5 * 5)), 0), xt::arange(5 * 5)))); + REQUIRE(xt::all(xt::equal(stitch.elemset(xt::eval(xt::arange(5 * 5)), 1), xt::arange(5 * 5) + 5 * 5))); REQUIRE(xt::all(xt::equal(stitch.nodeset(overlap_a, 0), stitch.nodeset(overlap_b, 1)))); REQUIRE(xt::all(xt::equal(stitch.nodeset({nset, nset}), xt::arange(0, 6 * 6 + 5 * 6, 6)))); REQUIRE(xt::all(xt::equal(stitch.elemset({eset, eset}), xt::arange(2 * 5 * 5)))); } SECTION("Vstack - equivalent test as 'Stitch'") { GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0); GooseFEM::Mesh::Vstack stitch; stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); auto nset = mesh.nodesLeftEdge(); auto eset = xt::eval(xt::arange(mesh.nelem())); GooseFEM::Mesh::Quad4::Regular res(5, 10, 1.0); REQUIRE(xt::allclose(stitch.coor(), res.coor())); REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn()))); REQUIRE(stitch.nodemap().size() == 2); REQUIRE(stitch.elemmap().size() == 2); REQUIRE(xt::all(xt::equal(stitch.nodemap(0), xt::arange(6 * 6)))); REQUIRE(xt::all(xt::equal(stitch.nodemap(1), xt::arange(6 * 6) + 5 * 6))); REQUIRE(xt::all(xt::equal(stitch.elemmap(0), xt::arange(5 * 5)))); REQUIRE(xt::all(xt::equal(stitch.elemmap(1), xt::arange(5 * 5) + 5 * 5))); - REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange(6 * 6), 0), xt::arange(6 * 6)))); - REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange(6 * 6), 1), xt::arange(6 * 6) + 5 * 6))); - REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange(5 * 5), 0), xt::arange(5 * 5)))); - REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange(5 * 5), 1), xt::arange(5 * 5) + 5 * 5))); + REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::eval(xt::arange(6 * 6)), 0), xt::arange(6 * 6)))); + REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::eval(xt::arange(6 * 6)), 1), xt::arange(6 * 6) + 5 * 6))); + REQUIRE(xt::all(xt::equal(stitch.elemset(xt::eval(xt::arange(5 * 5)), 0), xt::arange(5 * 5)))); + REQUIRE(xt::all(xt::equal(stitch.elemset(xt::eval(xt::arange(5 * 5)), 1), xt::arange(5 * 5) + 5 * 5))); REQUIRE(xt::all(xt::equal(stitch.nodeset({nset, nset}), xt::arange(0, 6 * 6 + 5 * 6, 6)))); REQUIRE(xt::all(xt::equal(stitch.elemset({eset, eset}), xt::arange(2 * 5 * 5)))); } SECTION("Vstack - several layers") { GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0); GooseFEM::Mesh::Vstack stitch; stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); stitch.push_back(mesh.coor(), mesh.conn(), mesh.nodesBottomEdge(), mesh.nodesTopEdge()); GooseFEM::Mesh::Quad4::Regular res(5, 5 * 5, 1.0); REQUIRE(xt::allclose(stitch.coor(), res.coor())); REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn()))); REQUIRE(stitch.nodemap().size() == 5); REQUIRE(stitch.elemmap().size() == 5); } SECTION("edgesize") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 10.0); auto s = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn()); auto t = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn(), mesh.getElementType()); REQUIRE(xt::allclose(s, 10.0)); REQUIRE(xt::allclose(t, 10.0)); } SECTION("coordination") { GooseFEM::Mesh::Quad4::Regular mesh(3, 3); auto N = GooseFEM::Mesh::coordination(mesh.conn()); xt::xtensor ret = {1, 2, 2, 1, 2, 4, 4, 2, 2, 4, 4, 2, 1, 2, 2, 1}; REQUIRE(xt::all(xt::equal(N, ret))); } SECTION("centers") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 2.0); xt::xtensor c = { {1.0, 1.0}, {3.0, 1.0}, {1.0, 3.0}, {3.0, 3.0}}; REQUIRE(xt::allclose(GooseFEM::Mesh::centers(mesh.coor(), mesh.conn()), c)); } SECTION("elem2node") { GooseFEM::Mesh::Quad4::Regular mesh(3, 3); auto tonode = GooseFEM::Mesh::elem2node(mesh.conn()); REQUIRE(tonode.size() == 16); REQUIRE(tonode[0] == std::vector{0}); REQUIRE(tonode[1] == std::vector{0, 1}); REQUIRE(tonode[2] == std::vector{1, 2}); REQUIRE(tonode[3] == std::vector{2}); REQUIRE(tonode[4] == std::vector{0, 3}); REQUIRE(tonode[5] == std::vector{0, 1, 3, 4}); REQUIRE(tonode[6] == std::vector{1, 2, 4, 5}); REQUIRE(tonode[7] == std::vector{2, 5}); REQUIRE(tonode[8] == std::vector{3, 6}); REQUIRE(tonode[9] == std::vector{3, 4, 6, 7}); REQUIRE(tonode[10] == std::vector{4, 5, 7, 8}); REQUIRE(tonode[11] == std::vector{5, 8}); REQUIRE(tonode[12] == std::vector{6}); REQUIRE(tonode[13] == std::vector{6, 7}); REQUIRE(tonode[14] == std::vector{7, 8}); REQUIRE(tonode[15] == std::vector{8}); } SECTION("elemmap2nodemap") { GooseFEM::Mesh::Quad4::Regular mesh(3, 3); xt::xtensor elmap0 = { 0, 1, 2, 3, 4, 5, 6, 7, 8 }; xt::xtensor elmap1 = { 2, 0, 1, 5, 3, 4, 8, 6, 7 }; xt::xtensor elmap2 = { 1, 2, 0, 4, 5, 3, 7, 8, 6 }; xt::xtensor nodemap0 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 }; xt::xtensor nodemap1 = { 2, 0, 1, 2, 6, 4, 5, 6, 10, 8, 9, 10, 14, 15, 13, 14 }; xt::xtensor nodemap2 = { 1, 2, 0, 1, 5, 6, 4, 5, 9, 10, 8, 9, 13, 14, 15, 13 }; REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap0, mesh.coor(), mesh.conn()), nodemap0))); REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap1, mesh.coor(), mesh.conn()), nodemap1))); REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap2, mesh.coor(), mesh.conn()), nodemap2))); } SECTION("elemmap2nodemap - example 1") { GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3); xt::xtensor elemval = { 1, 0, 0, 1, 0, 0, 1, 0, 0 }; xt::xtensor elemval_r1 = { 0, 1, 0, 0, 1, 0, 0, 1, 0 }; xt::xtensor elemval_r2 = { 0, 0, 1, 0, 0, 1, 0, 0, 1 }; xt::xtensor nodeval = { 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1 }; xt::xtensor nodeval_r1 = { 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0 }; xt::xtensor nodeval_r2 = { 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0 }; { auto elemmap = mesh.roll(0); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(1); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(2); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(3); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(4); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(5); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap))))); } } SECTION("elemmap2nodemap - example 2") { GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3); xt::xtensor elemval = { 1, 0, 0, 0, 1, 0, 0, 0, 1 }; xt::xtensor elemval_r1 = { 0, 1, 0, 0, 0, 1, 1, 0, 0 }; xt::xtensor elemval_r2 = { 0, 0, 1, 1, 0, 0, 0, 1, 0 }; xt::xtensor nodeval = { 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1 }; xt::xtensor nodeval_r1 = { 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0 }; xt::xtensor nodeval_r2 = { 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0 }; { auto elemmap = mesh.roll(0); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(1); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(2); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(3); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(4); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap))))); } { auto elemmap = mesh.roll(5); auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn()); REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap))))); REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap))))); } } }