diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h index 46c3790..b52fd16 100644 --- a/include/GooseFEM/Mesh.h +++ b/include/GooseFEM/Mesh.h @@ -1,590 +1,772 @@ /** Generic mesh operations. \file Mesh.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESH_H #define GOOSEFEM_MESH_H #include "config.h" namespace GooseFEM { namespace Mesh { - /** Enumerator for element-types */ enum class ElementType { + Unknown, ///< Unknown element-type Quad4, ///< Quadrilateral: 4-noded element in 2-d Hex8, ///< Hexahedron: 8-noded element in 3-d Tri3 ///< Triangle: 3-noded element in 2-d }; /** Extract the element type based on the connectivity. \param coor Nodal coordinates. \param conn Connectivity. \return ElementType(). */ inline ElementType defaultElementType( const xt::xtensor& coor, const xt::xtensor& conn); +/** +Base-class for regular meshes in 2-d. +This class does not have a specific element-type in mind, it is used mostly internally +to derive from such that common methods do not have to be reimplementation. +*/ +class RegularBase2d { +public: + + RegularBase2d() = default; + + /** + \return Number of elements. + */ + size_t nelem() const; + + /** + \return Number of nodes. + */ + size_t nnode() const; + + /** + \return Number of nodes-per-element == 4. + */ + size_t nne() const; + + /** + \return Number of dimensions == 2. + */ + size_t ndim() const; + + /** + \return Number of elements in x-direction == width of the mesh in units of #h. + */ + size_t nelx() const; + + /** + \return Number of elements in y-direction == height of the mesh, in units of #h, + */ + size_t nely() const; + + /** + \return Linear edge size of one 'block'. + */ + double h() const; + + /** + \return The ElementType(). + */ + virtual ElementType getElementType() const; + + /** + \return Nodal coordinates [#nnode, #ndim]. + */ + virtual xt::xtensor coor() const; + + /** + \return Connectivity [#nelem, #nne]. + */ + virtual xt::xtensor conn() const; + + /** + \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. + */ + virtual xt::xtensor dofs() const; + + /** + Nodes along the bottom edge (y = 0), in order of increasing x. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottomEdge() const; + + /** + Nodes along the top edge (y = #nely * #h), in order of increasing x. + + \return List of node numbers. + */ + virtual xt::xtensor nodesTopEdge() const; + + /** + Nodes along the left edge (x = 0), in order of increasing y. + + \return List of node numbers. + */ + virtual xt::xtensor nodesLeftEdge() const; + + /** + Nodes along the right edge (x = #nelx * #h), in order of increasing y. + + \return List of node numbers. + */ + virtual xt::xtensor nodesRightEdge() const; + + /** + Nodes along the bottom edge (y = 0), without the corners (at x = 0 and x = #nelx * #h). + Same as: nodesBottomEdge()[1: -1]. + + \return List of node numbers. + */ + virtual xt::xtensor nodesBottomOpenEdge() const; + + /** + Nodes along the top edge (y = #nely * #h), without the corners (at x = 0 and x = #nelx * #h). + Same as: nodesTopEdge()[1: -1]. + + \return List of node numbers. + */ + virtual xt::xtensor nodesTopOpenEdge() const; + + /** + Nodes along the left edge (x = 0), without the corners (at y = 0 and y = #nely * #h). + Same as: nodesLeftEdge()[1: -1]. + + \return List of node numbers. + */ + virtual xt::xtensor nodesLeftOpenEdge() const; + + /** + Nodes along the right edge (x = #nelx * #h), without the corners (at y = 0 and y = #nely * #h). + Same as: nodesRightEdge()[1: -1]. + + \return List of node numbers. + */ + virtual xt::xtensor nodesRightOpenEdge() const; + + /** + The bottom-left corner node (at x = 0, y = 0). + Same as nodesBottomEdge()[0] and nodesLeftEdge()[0]. + + \return Node number. + */ + virtual size_t nodesBottomLeftCorner() const; + + /** + The bottom-right corner node (at x = #nelx * #h, y = 0). + Same as nodesBottomEdge()[-1] and nodesRightEdge()[0]. + + \return Node number. + */ + virtual size_t nodesBottomRightCorner() const; + + /** + The top-left corner node (at x = 0, y = #nely * #h). + Same as nodesTopEdge()[0] and nodesRightEdge()[-1]. + + \return Node number. + */ + virtual size_t nodesTopLeftCorner() const; + + /** + The top-right corner node (at x = #nelx * #h, y = #nely * #h). + Same as nodesTopEdge()[-1] and nodesRightEdge()[-1]. + + \return Node number. + */ + virtual size_t nodesTopRightCorner() const; + + /** + \return Alias of nodesBottomLeftCorner(). + */ + size_t nodesLeftBottomCorner() const; + + /** + \return Alias of nodesTopLeftCorner(). + */ + size_t nodesLeftTopCorner() const; + + /** + \return Alias of nodesBottomRightCorner(). + */ + size_t nodesRightBottomCorner() const; + + /** + \return Alias of nodesTopRightCorner(). + */ + size_t nodesRightTopCorner() const; + + /** + DOF-numbers for the case that the periodicity if fully eliminated. Such that: + + dofs[nodesRightOpenEdge(), :] = dofs[nodesLeftOpenEdge(), :] + dofs[nodesTopOpenEdge(), :] = dofs[nodesBottomOpenEdge(), :] + dofs[nodesBottomRightCorner(), :] = dofs[nodesBottomLeftCorner(), :] + dofs[nodesTopRightCorner(), :] = dofs[nodesBottomLeftCorner(), :] + dofs[nodesTopLeftCorner(), :] = dofs[nodesBottomLeftCorner(), :] + + \return DOF numbers for each node [#nnode, #ndim]. + */ + xt::xtensor dofsPeriodic() const; + + /** + Periodic node pairs, in two columns: (independent, dependent) + + \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() +}; + /** Find overlapping nodes. The output has the following structure: [[nodes_from_mesh_a], [nodes_from_mesh_b]] \param coor_a Nodal coordinates of mesh "a". \param coor_b Nodal coordinates of mesh "b". \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. \return Overlapping nodes. */ inline xt::xtensor overlapping( const xt::xtensor& coor_a, const xt::xtensor& coor_b, double rtol = 1e-5, double atol = 1e-8); /** Stitch two mesh objects, specifying overlapping nodes by hand. */ class ManualStitch { public: ManualStitch() = default; /** \param coor_a Nodal coordinates of mesh "a". \param conn_a Connectivity of mesh "a". \param overlapping_nodes_a Node-numbers of mesh "a" that overlap with mesh "b". \param coor_b Nodal coordinates of mesh "b". \param conn_b Connectivity of mesh "b". \param overlapping_nodes_b Node-numbers of mesh "b" that overlap with mesh "a". \param check_position If ``true`` the nodes are checked for position overlap. \param rtol Relative tolerance for check on position overlap. \param atol Absolute tolerance for check on position overlap. */ ManualStitch( const xt::xtensor& coor_a, const xt::xtensor& conn_a, const xt::xtensor& overlapping_nodes_a, const xt::xtensor& coor_b, const xt::xtensor& conn_b, const xt::xtensor& overlapping_nodes_b, bool check_position = true, double rtol = 1e-5, double atol = 1e-8); /** - Number of sub meshes. - - \return 2. + \return Number of sub meshes == 2. */ size_t nmesh() const; /** - Number of elements. - - \return unsigned int. + \return Number of elements. */ size_t nelem() const; /** - Number of nodes. - - \return unsigned int. + \return Number of nodes. */ size_t nnode() const; /** - Number of nodes-per-element. - - \return unsigned int. + \return Number of nodes-per-element. */ size_t nne() const; /** - Number of dimensions. - - \return unsigned int. + \return Number of dimensions. */ size_t ndim() const; /** - Nodal coordinates. - - \return [#nnode, #ndim]. + \return Nodal coordinates [#nnode, #ndim]. */ xt::xtensor coor() const; /** - Connectivity. - - \return [#nelem, #nne]. + \return Connectivity [#nelem, #nne]. */ xt::xtensor conn() const; /** - DOF numbers for each node (numbered sequentially). - - \return [#nnode, #ndim]. + \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. */ xt::xtensor dofs() const; /** \return Node-map per sub-mesh. */ std::vector> nodemap() const; /** \return Element-map per sub-mesh. */ std::vector> elemmap() const; /** \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return Node-map for a given mesh. */ xt::xtensor nodemap(size_t mesh_index) const; /** \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return Element-map for a given mesh. */ xt::xtensor elemmap(size_t mesh_index) const; /** Convert set of node numbers for an original mesh to the stitched mesh. \param set List of node numbers. \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return List of node numbers for the stitched mesh. */ xt::xtensor nodeset(const xt::xtensor& set, size_t mesh_index) const; /** Convert set of element numbers for an original mesh to the stitched mesh. \param set List of element numbers. \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return List of element numbers for the stitched mesh. */ xt::xtensor elemset(const xt::xtensor& set, size_t mesh_index) const; private: xt::xtensor m_coor; xt::xtensor m_conn; xt::xtensor m_map_b; size_t m_nnd_a; size_t m_nel_a; size_t m_nel_b; }; /** Stitch mesh objects, automatically searching for overlapping nodes. */ class Stitch { public: /** \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. */ Stitch(double rtol = 1e-5, double atol = 1e-8); /** Add mesh to be stitched. \param coor Nodal coordinates. \param conn Connectivity. */ void push_back(const xt::xtensor& coor, const xt::xtensor& conn); /** - Number of sub meshes. - - \return unsigned int + \return Number of sub meshes. */ size_t nmesh() const; /** - Number of elements. - - \return unsigned int. + \return Number of elements. */ size_t nelem() const; /** - Number of nodes. - - \return unsigned int. + \return Number of nodes. */ size_t nnode() const; /** - Number of nodes-per-element. - - \return unsigned int. + \return Number of nodes-per-element. */ size_t nne() const; /** - Number of dimensions. - - \return unsigned int. + \return Number of dimensions. */ size_t ndim() const; /** - Nodal coordinates. - - \return [#nnode, #ndim]. + \return Nodal coordinates [#nnode, #ndim]. */ xt::xtensor coor() const; /** - Connectivity. - - \return [#nelem, #nne]. + \return Connectivity [#nelem, #nne]. */ xt::xtensor conn() const; /** - DOF numbers for each node (numbered sequentially). - - \return [#nnode, #ndim]. + \return DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. */ xt::xtensor dofs() const; /** \return Node-map per sub-mesh. */ std::vector> nodemap() const; /** \return Element-map per sub-mesh. */ std::vector> elemmap() const; /** The node numbers in the stitched mesh that are coming from a specific sub-mesh. \param mesh_index Index of the sub-mesh. \return List of node numbers. */ xt::xtensor nodemap(size_t mesh_index) const; /** The element numbers in the stitched mesh that are coming from a specific sub-mesh. \param mesh_index Index of the sub-mesh. \return List of element numbers. */ xt::xtensor elemmap(size_t mesh_index) const; /** Convert set of node-numbers for a sub-mesh to the stitched mesh. \param set List of node numbers. \param mesh_index Index of the sub-mesh. \return List of node numbers for the stitched mesh. */ xt::xtensor nodeset(const xt::xtensor& set, size_t mesh_index) const; /** Convert set of element-numbers for a sub-mesh to the stitched mesh. \param set List of element numbers. \param mesh_index Index of the sub-mesh. \return List of element numbers for the stitched mesh. */ xt::xtensor elemset(const xt::xtensor& set, size_t mesh_index) const; /** Combine set of node numbers for an original to the final mesh (removes duplicates). \param set List of node numbers per mesh. \return List of node numbers for the stitched mesh. */ xt::xtensor nodeset(const std::vector>& set) const; /** Combine set of element numbers for an original to the final mesh. \param set List of element numbers per mesh. \return List of element numbers for the stitched mesh. */ xt::xtensor elemset(const std::vector>& set) const; private: xt::xtensor m_coor; xt::xtensor m_conn; std::vector> m_map; std::vector m_nel; ///< Number of elements per sub-mesh. std::vector m_el_offset; double m_rtol; double m_atol; }; /** \rst Renumber indices to lowest possible index. For example: .. math:: \begin{bmatrix} 0 & 1 \\ 5 & 4 \end{bmatrix} is renumbered to .. math:: \begin{bmatrix} 0 & 1 \\ 3 & 2 \end{bmatrix} Or, in pseudo-code, the result of this function is that: .. code-block:: python dofs = renumber(dofs) sort(unique(dofs[:])) == range(max(dofs+1)) .. tip:: One can use the wrapper function :cpp:func:`GooseFEM::Mesh::renumber`. This class gives more advanced features. \endrst */ class Renumber { public: Renumber() = default; /** \param dofs DOF-numbers. */ template Renumber(const T& dofs); /** Get renumbered DOFs (same as ``Renumber::apply(dofs)``). \param dofs List of (DOF-)numbers. \return Renumbered list of (DOF-)numbers. */ [[deprecated]] xt::xtensor get(const xt::xtensor& dofs) const; /** Apply renumbering to other set. \param list List of (DOF-)numbers. \return Renumbered list of (DOF-)numbers. */ template T apply(const T& list) const; /** Get the list needed to renumber, e.g.: dofs_renumbered(i, j) = index(dofs(i, j)) \return Renumber-index. */ xt::xtensor index() const; private: xt::xtensor m_renum; }; /** Renumber to lowest possible index (see GooseFEM::Mesh::Renumber). \param dofs DOF-numbers. \return Renumbered DOF-numbers. */ inline xt::xtensor renumber(const xt::xtensor& dofs); /** Reorder to lowest possible index, in specific order. For example for ``Reorder({iiu, iip})`` after reordering: iiu = xt::range(nnu); iip = xt::range(nnp) + nnu; */ class Reorder { public: Reorder() = default; /** \param args List of (DOF-)numbers. */ Reorder(const std::initializer_list> args); /** Get reordered DOFs (same as ``Reorder::apply(dofs)``). \param dofs List of (DOF-)numbers. \return Reordered list of (DOF-)numbers. */ [[deprecated]] xt::xtensor get(const xt::xtensor& dofs) const; /** Apply reordering to other set. \param list List of (DOF-)numbers. \return Reordered list of (DOF-)numbers. */ template T apply(const T& list) const; /** Get the list needed to reorder, e.g.: dofs_reordered(i, j) = index(dofs(i, j)) \return Reorder-index. */ xt::xtensor index() const; private: xt::xtensor m_renum; }; /** List with DOF-numbers in sequential order. The output is a sequential list of DOF-numbers for each vector-component of each node. For example for 3 nodes in 2 dimensions the output is \rst .. math:: \begin{bmatrix} 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \endrst \param nnode Number of nodes. \param ndim Number of dimensions. \return DOF-numbers. */ inline xt::xtensor dofs(size_t nnode, size_t ndim); /** Number of elements connected to each node. \param conn Connectivity. \return Coordination per node. */ inline xt::xtensor coordination(const xt::xtensor& conn); /** Elements connected to each node. \param conn Connectivity. \param sorted If ``true`` the output is sorted. \return Elements per node. */ inline std::vector> elem2node( const xt::xtensor& conn, bool sorted = true); /** Return size of each element edge. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Edge-sizes per element. */ inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type); /** Return size of each element edge. The element-type is automatically determined, see defaultElementType(). \param coor Nodal coordinates. \param conn Connectivity. \return Edge-sizes per element. */ inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn); /** Coordinates of the center of each element. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Center of each element. */ inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type); /** Coordinates of the center of each element. The element-type is automatically determined, see defaultElementType(). \param coor Nodal coordinates. \param conn Connectivity. \return Center of each element. */ inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn); /** Convert an element-map to a node-map. \param elem_map Element-map such that ``new_elvar = elvar[elem_map]``. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Node-map such that ``new_nodevar = nodevar[node_map]`` */ inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn, ElementType type); /** Convert an element-map to a node-map. The element-type is automatically determined, see defaultElementType(). \param elem_map Element-map such that ``new_elvar = elvar[elem_map]``. \param coor Nodal coordinates. \param conn Connectivity. \return Node-map such that ``new_nodevar = nodevar[node_map]`` */ inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn); } // namespace Mesh } // namespace GooseFEM #include "Mesh.hpp" #endif diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp index e01110e..9259e51 100644 --- a/include/GooseFEM/Mesh.hpp +++ b/include/GooseFEM/Mesh.hpp @@ -1,656 +1,841 @@ /** Implementation of Mesh.h \file Mesh.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESH_HPP #define GOOSEFEM_MESH_HPP #include "Mesh.h" namespace GooseFEM { namespace Mesh { inline ElementType defaultElementType( const xt::xtensor& coor, const xt::xtensor& conn) { if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) { return ElementType::Tri3; } if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) { return ElementType::Quad4; } if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) { return ElementType::Hex8; } throw std::runtime_error("Element-type not implemented"); } +inline size_t RegularBase2d::nelem() const +{ + return m_nelem; +} + +inline size_t RegularBase2d::nnode() const +{ + return m_nnode; +} + +inline size_t RegularBase2d::nne() const +{ + return m_nne; +} + +inline size_t RegularBase2d::ndim() const +{ + return m_ndim; +} + +inline size_t RegularBase2d::nelx() const +{ + return m_nelx; +} + +inline size_t RegularBase2d::nely() const +{ + return m_nely; +} + +inline double RegularBase2d::h() const +{ + return m_h; +} + +inline ElementType RegularBase2d::getElementType() const +{ + return ElementType::Unknown; +} + +inline xt::xtensor RegularBase2d::coor() const +{ + return xt::empty({0, 0}); +} + +inline xt::xtensor RegularBase2d::conn() const +{ + return xt::empty({0, 0}); +} + +inline xt::xtensor RegularBase2d::nodesBottomEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesTopEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesLeftEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesRightEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesBottomOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesTopOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesLeftOpenEdge() const +{ + return xt::empty({0}); +} + +inline xt::xtensor RegularBase2d::nodesRightOpenEdge() const +{ + return xt::empty({0}); +} + +inline size_t RegularBase2d::nodesBottomLeftCorner() const +{ + return 0; +} + +inline size_t RegularBase2d::nodesBottomRightCorner() const +{ + return 0; +} + +inline size_t RegularBase2d::nodesTopLeftCorner() const +{ + return 0; +} + +inline size_t RegularBase2d::nodesTopRightCorner() const +{ + return 0; +} + +inline size_t RegularBase2d::nodesLeftBottomCorner() const +{ + return nodesBottomLeftCorner(); +} + +inline size_t RegularBase2d::nodesLeftTopCorner() const +{ + return nodesTopLeftCorner(); +} + +inline size_t RegularBase2d::nodesRightBottomCorner() const +{ + return nodesBottomRightCorner(); +} + +inline size_t RegularBase2d::nodesRightTopCorner() const +{ + return 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); +} + namespace detail { template inline T renum(const T& arg, const R& mapping) { T ret = T::from_shape(arg.shape()); auto jt = ret.begin(); for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) { *jt = mapping(*it); } return ret; } } // namespace detail inline xt::xtensor overlapping( const xt::xtensor& coor_a, const xt::xtensor& coor_b, double rtol, double atol) { GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); std::vector ret_a; std::vector ret_b; for (size_t i = 0; i < coor_a.shape(0); ++i) { auto idx = xt::flatten_indices(xt::argwhere(xt::prod(xt::isclose( coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1))); for (auto& j : idx) { ret_a.push_back(i); ret_b.push_back(j); } } xt::xtensor ret = xt::empty({size_t(2), ret_a.size()}); for (size_t i = 0; i < ret_a.size(); ++i) { ret(0, i) = ret_a[i]; ret(1, i) = ret_b[i]; } return ret; } inline ManualStitch::ManualStitch( const xt::xtensor& coor_a, const xt::xtensor& conn_a, const xt::xtensor& overlapping_nodes_a, const xt::xtensor& coor_b, const xt::xtensor& conn_b, const xt::xtensor& overlapping_nodes_b, bool check_position, double rtol, double atol) { UNUSED(rtol); UNUSED(atol); GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape())); GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1)); if (check_position) { GOOSEFEM_ASSERT(xt::allclose( xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()), xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()), rtol, atol)); } size_t nnda = coor_a.shape(0); size_t nndb = coor_b.shape(0); size_t ndim = coor_a.shape(1); size_t nelim = overlapping_nodes_a.size(); size_t nela = conn_a.shape(0); size_t nelb = conn_b.shape(0); size_t nne = conn_a.shape(1); m_nel_a = nela; m_nel_b = nelb; m_nnd_a = nnda; xt::xtensor keep_b = xt::setdiff1d(xt::arange(nndb), overlapping_nodes_b); m_map_b = xt::empty({nndb}); xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a; xt::view(m_map_b, xt::keep(keep_b)) = xt::arange(keep_b.size()) + nnda; m_conn = xt::empty({nela + nelb, nne}); xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a; xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b); m_coor = xt::empty({nnda + nndb - nelim, ndim}); xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a; xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) = xt::view(coor_b, xt::keep(keep_b), xt::all()); } inline xt::xtensor ManualStitch::coor() const { return m_coor; } inline xt::xtensor ManualStitch::conn() const { return m_conn; } inline size_t ManualStitch::nmesh() const { return 2; } inline size_t ManualStitch::nelem() const { return m_conn.shape(0); } inline size_t ManualStitch::nnode() const { return m_coor.shape(0); } inline size_t ManualStitch::nne() const { return m_conn.shape(1); } inline size_t ManualStitch::ndim() const { return m_coor.shape(1); } inline xt::xtensor ManualStitch::dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline std::vector> ManualStitch::nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } inline std::vector> ManualStitch::elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } inline xt::xtensor ManualStitch::nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nnd_a); } return m_map_b; } inline xt::xtensor ManualStitch::elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nel_a); } return xt::arange(m_nel_b) + m_nel_a; } inline xt::xtensor ManualStitch::nodeset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a); return set; } GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size()); return detail::renum(set, m_map_b); } inline xt::xtensor ManualStitch::elemset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a); return set; } GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b); return set + m_nel_a; } inline Stitch::Stitch(double rtol, double atol) { m_rtol = rtol; m_atol = atol; } inline void Stitch::push_back( const xt::xtensor& coor, const xt::xtensor& conn) { if (m_map.size() == 0) { m_coor = coor; m_conn = conn; m_map.push_back(xt::eval(xt::arange(coor.shape(0)))); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(0); return; } auto overlap = overlapping(m_coor, coor, m_rtol, m_atol); size_t index = m_map.size(); ManualStitch stich( m_coor, m_conn, xt::view(overlap, 0, xt::all()), coor, conn, xt::view(overlap, 1, xt::all()), false); m_coor = stich.coor(); m_conn = stich.conn(); m_map.push_back(stich.nodemap(1)); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]); } inline size_t Stitch::nmesh() const { return m_map.size(); } inline xt::xtensor Stitch::coor() const { return m_coor; } inline xt::xtensor Stitch::conn() const { return m_conn; } inline size_t Stitch::nelem() const { return m_conn.shape(0); } inline size_t Stitch::nnode() const { return m_coor.shape(0); } inline size_t Stitch::nne() const { return m_conn.shape(1); } inline size_t Stitch::ndim() const { return m_coor.shape(1); } inline xt::xtensor Stitch::dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline std::vector> Stitch::nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } inline std::vector> Stitch::elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } inline xt::xtensor Stitch::nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return m_map[mesh_index]; } inline xt::xtensor Stitch::elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return xt::arange(m_nel[mesh_index]) + m_el_offset[mesh_index]; } inline xt::xtensor Stitch::nodeset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size()); return detail::renum(set, m_map[mesh_index]); } inline xt::xtensor Stitch::elemset(const xt::xtensor& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]); return set + m_el_offset[mesh_index]; } inline xt::xtensor Stitch::nodeset(const std::vector>& set) const { GOOSEFEM_ASSERT(set.size() == m_map.size()); size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->nodeset(set[i], i); n += set[i].size(); } return xt::unique(ret); } inline xt::xtensor Stitch::elemset(const std::vector>& set) const { GOOSEFEM_ASSERT(set.size() == m_map.size()); size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i); n += set[i].size(); } return ret; } template inline Renumber::Renumber(const T& dofs) { size_t n = xt::amax(dofs)() + 1; size_t i = 0; xt::xtensor unique = xt::unique(dofs); m_renum = xt::empty({n}); for (auto& j : unique) { m_renum(j) = i; ++i; } } inline xt::xtensor Renumber::get(const xt::xtensor& dofs) const { GOOSEFEM_WARNING("Renumber::get is deprecated, use Renumber::apply"); return this->apply(dofs); } template inline T Renumber::apply(const T& list) const { return detail::renum(list, m_renum); } inline xt::xtensor Renumber::index() const { return m_renum; } inline xt::xtensor renumber(const xt::xtensor& dofs) { return Renumber(dofs).apply(dofs); } inline Reorder::Reorder(const std::initializer_list> args) { size_t n = 0; size_t i = 0; for (auto& arg : args) { if (arg.size() == 0) { continue; } n = std::max(n, xt::amax(arg)() + 1); } #ifdef GOOSEFEM_ENABLE_ASSERT for (auto& arg : args) { GOOSEFEM_ASSERT(xt::unique(arg) == xt::sort(arg)); } #endif m_renum = xt::empty({n}); for (auto& arg : args) { for (auto& j : arg) { m_renum(j) = i; ++i; } } } inline xt::xtensor Reorder::get(const xt::xtensor& dofs) const { GOOSEFEM_WARNING("Reorder::get is deprecated, use Reorder::apply"); return this->apply(dofs); } template inline T Reorder::apply(const T& list) const { T ret = T::from_shape(list.shape()); auto jt = ret.begin(); for (auto it = list.begin(); it != list.end(); ++it, ++jt) { *jt = m_renum(*it); } return ret; } inline xt::xtensor Reorder::index() const { return m_renum; } inline xt::xtensor dofs(size_t nnode, size_t ndim) { return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline xt::xtensor coordination(const xt::xtensor& conn) { size_t nnode = xt::amax(conn)() + 1; xt::xtensor N = xt::zeros({nnode}); for (auto it = conn.begin(); it != conn.end(); ++it) { N(*it) += 1; } return N; } inline std::vector> elem2node(const xt::xtensor& conn, bool sorted) { auto N = coordination(conn); auto nnode = N.size(); std::vector> ret(nnode); for (size_t i = 0; i < nnode; ++i) { ret[i].reserve(N(i)); } for (size_t e = 0; e < conn.shape(0); ++e) { for (size_t m = 0; m < conn.shape(1); ++m) { ret[conn(e, m)].push_back(e); } } if (sorted) { for (auto& row : ret) { std::sort(row.begin(), row.end()); } } return ret; } inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); if (type == ElementType::Quad4) { GOOSEFEM_ASSERT(coor.shape(1) == 2ul); GOOSEFEM_ASSERT(conn.shape(1) == 4ul); xt::xtensor n0 = xt::view(conn, xt::all(), 0); xt::xtensor n1 = xt::view(conn, xt::all(), 1); xt::xtensor n2 = xt::view(conn, xt::all(), 2); xt::xtensor n3 = xt::view(conn, xt::all(), 3); xt::xtensor x0 = xt::view(coor, xt::keep(n0), 0); xt::xtensor x1 = xt::view(coor, xt::keep(n1), 0); xt::xtensor x2 = xt::view(coor, xt::keep(n2), 0); xt::xtensor x3 = xt::view(coor, xt::keep(n3), 0); xt::xtensor y0 = xt::view(coor, xt::keep(n0), 1); xt::xtensor y1 = xt::view(coor, xt::keep(n1), 1); xt::xtensor y2 = xt::view(coor, xt::keep(n2), 1); xt::xtensor y3 = xt::view(coor, xt::keep(n3), 1); xt::xtensor ret = xt::empty(conn.shape()); xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0)); xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0)); xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0)); xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0)); return ret; } throw std::runtime_error("Element-type not implemented"); } inline xt::xtensor edgesize( const xt::xtensor& coor, const xt::xtensor& conn) { return edgesize(coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); xt::xtensor ret = xt::zeros({conn.shape(0), coor.shape(1)}); if (type == ElementType::Quad4) { GOOSEFEM_ASSERT(coor.shape(1) == 2); GOOSEFEM_ASSERT(conn.shape(1) == 4); for (size_t i = 0; i < 4; ++i) { auto n = xt::view(conn, xt::all(), i); ret += xt::view(coor, xt::keep(n), xt::all()); } ret /= 4.0; return ret; } throw std::runtime_error("Element-type not implemented"); } inline xt::xtensor centers( const xt::xtensor& coor, const xt::xtensor& conn) { return centers(coor, conn, defaultElementType(coor, conn)); } inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn, ElementType type) { GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0)); size_t N = coor.shape(0); xt::xtensor ret = N * xt::ones({N}); if (type == ElementType::Quad4) { GOOSEFEM_ASSERT(coor.shape(1) == 2); GOOSEFEM_ASSERT(conn.shape(1) == 4); for (size_t i = 0; i < 4; ++i) { xt::xtensor t = N * xt::ones({N}); auto old_nd = xt::view(conn, xt::all(), i); auto new_nd = xt::view(conn, xt::keep(elem_map), i); xt::view(t, xt::keep(old_nd)) = new_nd; ret = xt::where(xt::equal(ret, N), t, ret); } return ret; } throw std::runtime_error("Element-type not implemented"); } inline xt::xtensor elemmap2nodemap( const xt::xtensor& elem_map, const xt::xtensor& coor, const xt::xtensor& conn) { return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn)); } } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h index 7c5421d..8a6ce7a 100644 --- a/include/GooseFEM/MeshQuad4.h +++ b/include/GooseFEM/MeshQuad4.h @@ -1,626 +1,493 @@ /** Generate simple meshes of 4-noded quadrilateral elements in 2d (GooseFEM::Mesh::ElementType::Quad4). \file MeshQuad4.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHQUAD4_H #define GOOSEFEM_MESHQUAD4_H #include "config.h" +#include "Mesh.h" namespace GooseFEM { namespace Mesh { + +/** +Simple meshes of quadrilateral elements of type ElementType::Quad4. +*/ namespace Quad4 { // pre-allocation - namespace Map { class FineLayer2Regular; } /** Regular mesh: equi-sized elements. */ -class Regular { +class Regular : public RegularBase2d { public: Regular() = default; /** Constructor. \param nelx Number of elements in horizontal (x) direction. \param nely Number of elements in vertical (y) direction. \param h Edge size (width == height). */ Regular(size_t nelx, size_t nely, double h = 1.0); - /** - Number of elements. - - \return unsigned int. - */ - size_t nelem() const; + ElementType getElementType() const override; + xt::xtensor coor() const override; + xt::xtensor conn() const override; + xt::xtensor nodesBottomEdge() const override; + xt::xtensor nodesTopEdge() const override; + xt::xtensor nodesLeftEdge() const override; + xt::xtensor nodesRightEdge() const override; + xt::xtensor nodesBottomOpenEdge() const override; + xt::xtensor nodesTopOpenEdge() const override; + xt::xtensor nodesLeftOpenEdge() const override; + xt::xtensor nodesRightOpenEdge() const override; + size_t nodesBottomLeftCorner() const override; + size_t nodesBottomRightCorner() const override; + size_t nodesTopLeftCorner() const override; + size_t nodesTopRightCorner() const override; /** - Number of nodes. + Element numbers as 'matrix'. - \return unsigned int. + \return [#nely, #nelx]. */ - 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; - - /** - Number of elements in x-direction == width of the mesh in units of #h. - - \return unsigned int. - */ - size_t nelx() const; - - /** - Number of elements in y-direction == height of the mesh, in units of #h, - - \return unsigned int. - */ - size_t nely() const; - - /** - Edge size of one element. - - \return double. - */ - double h() const; - - /** - Element type. - - \return GooseFEM::Mesh::ElementType(). - */ - ElementType getElementType() const; - - /** - Nodal coordinates. - - \return [#nnode, #ndim]. - */ - xt::xtensor coor() const; - - /** - Connectivity. - - \return [#nelem, #nne]. - */ - xt::xtensor conn() const; - - /** - Nodes along the bottom edge (y = 0). - - \return List of node numbers. - */ - xt::xtensor nodesBottomEdge() const; - - /** - Nodes along the top edge (y = #nely * #h). - - \return List of node numbers. - */ - xt::xtensor nodesTopEdge() const; - - /** - Nodes along the left edge (x = 0). - - \return List of node numbers. - */ - xt::xtensor nodesLeftEdge() const; - - /** - Nodes along the right edge (x = #nelx * #h). - - \return List of node numbers. - */ - xt::xtensor nodesRightEdge() const; - - /** - Nodes along the bottom edge (y = 0), with the corners (at x = 0 and x = #nelx * #h). - - \return List of node numbers. - */ - xt::xtensor nodesBottomOpenEdge() const; - - /** - Nodes along the top edge (y = #nely * #h), with the corners (at x = 0 and x = #nelx * #h). - - \return List of node numbers. - */ - xt::xtensor nodesTopOpenEdge() const; - - /** - Nodes along the left edge (x = 0), with the corners (at y = 0 and y = #nely * #h). - - \return List of node numbers. - */ - xt::xtensor nodesLeftOpenEdge() const; - - /** - Nodes along the right edge (x = #nelx * #h), with the corners (at y = 0 and y = #nely * #h). - - \return List of node numbers. - */ - xt::xtensor nodesRightOpenEdge() const; - - // boundary nodes: corners (including aliases) - size_t nodesBottomLeftCorner() const; - size_t nodesBottomRightCorner() const; - size_t nodesTopLeftCorner() const; - size_t nodesTopRightCorner() const; - size_t nodesLeftBottomCorner() const; - size_t nodesLeftTopCorner() const; - size_t nodesRightBottomCorner() const; - size_t nodesRightTopCorner() const; - - // DOF-numbers for each component of each node (sequential) - xt::xtensor dofs() const; - - // DOF-numbers for the case that the periodicity if fully eliminated - xt::xtensor dofsPeriodic() const; - - // periodic node pairs [:,2]: (independent, dependent) - xt::xtensor nodesPeriodic() const; - - // front-bottom-left node, used as reference for periodicity - size_t nodesOrigin() const; - - // element numbers as matrix xt::xtensor elementgrid() const; - -private: - double m_h; // elementary element edge-size (in all directions) - size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") - size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") - size_t m_nelem; // number of elements - size_t m_nnode; // number of nodes - static const size_t m_nne = 4; // number of nodes-per-element - static const size_t m_ndim = 2; // number of dimensions }; // Mesh with fine middle layer, and coarser elements towards the top and bottom class FineLayer { public: FineLayer() = default; FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1); /** Reconstruct class for given coordinates / connectivity. \param coor Nodal coordinates ``[nnode, ndim]`` with ``ndim == 2``. \param conn Connectivity ``[nne, nne]`` with ``nne == 4``. \throw GOOSEFEM_CHECK() */ FineLayer(const xt::xtensor& coor, const xt::xtensor& conn); /** 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; /** Number of elements in x-direction along the middle layer. By definition equal to the width of the mesh in units of #h. \return unsigned int. */ size_t nelx() const; /** Height of the mesh, in units of #h. \return unsigned int. */ size_t nely() const; /** Edge size of the smallest elements (along the middle layer). \return double. */ double h() const; /** Edge size in x-direction of a block, in units of #h, per row of blocks. Note that a block is equal to an element except in refinement layers where it contains three elements. \returns List of size equal to the number of rows of blocks. */ xt::xtensor elemrow_nhx() const; /** Edge size in y-direction of a block, in units of #h, per row of blocks. Note that a block is equal to an element except in refinement layers where it contains three elements. \returns List of size equal to the number of rows of blocks. */ xt::xtensor elemrow_nhy() const; /** Number of elements per row of blocks. Note that a block is equal to an element except in refinement layers where it contains three elements. \returns List of size equal to the number of rows of blocks. */ xt::xtensor elemrow_nelem() const; /** Element type. \return GooseFEM::Mesh::ElementType(). */ ElementType getElementType() const; /** Nodal coordinates. \return ``[nnode, ndim]``. */ xt::xtensor coor() const; /** Connectivity. \return ``[nelem, nne]``. */ xt::xtensor conn() const; // elements in the middle (fine) layer xt::xtensor elementsMiddleLayer() const; // extract elements along a layer xt::xtensor elementsLayer(size_t layer) const; // select region of elements from 'matrix' of element numbers xt::xtensor elementgrid_ravel( std::vector rows_start_stop, std::vector cols_start_stop) const; /** Select region of elements from 'matrix' of element numbers around an element: square box with edge-size ``(2 * size + 1) * h``, around ``element``. \param element The element around which to select elements. \param size Edge size of the square box encapsulating the selected element. \param periodic Assume the mesh periodic. \returns List of elements. */ xt::xtensor elementgrid_around_ravel( size_t element, size_t size, bool periodic = true); /** Select region of elements from 'matrix' of element numbers around an element: left/right from ``element`` (on the same layer). \param element The element around which to select elements. \param left Number of elements to select to the left. \param right Number of elements to select to the right. \param periodic Assume the mesh periodic. \returns List of elements. */ // - xt::xtensor elementgrid_leftright( size_t element, size_t left, size_t right, bool periodic = true); // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor nodesRightOpenEdge() const; // boundary nodes: corners (including aliases) size_t nodesBottomLeftCorner() const; size_t nodesBottomRightCorner() const; size_t nodesTopLeftCorner() const; size_t nodesTopRightCorner() const; size_t nodesLeftBottomCorner() const; size_t nodesLeftTopCorner() const; size_t nodesRightBottomCorner() const; size_t nodesRightTopCorner() const; // DOF-numbers for each component of each node (sequential) xt::xtensor dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; // mapping to 'roll' periodically in the x-direction, // returns element mapping, such that: new_elemvar = elemvar[elem_map] xt::xtensor roll(size_t n); private: double m_h; // elementary element edge-size (in all directions) double m_Lx; // mesh size in "x" size_t m_nelem; // number of elements size_t m_nnode; // number of nodes static const size_t m_nne = 4; // number of nodes-per-element static const size_t m_ndim = 2; // number of dimensions xt::xtensor m_nelx; ///< See elemrow_nelem(). xt::xtensor m_nhx; ///< See elemrow_nhx(). xt::xtensor m_nhy; ///< See elemrow_nhy(). xt::xtensor m_nnd; // total number of nodes in the main node layer (**) xt::xtensor m_refine; // refine direction (-1:no refine, 0:"x" (*) xt::xtensor m_startElem; // start element (*) xt::xtensor m_startNode; // start node (**) // (*) per element layer in "y" // (**) per node layer in "y" void init(size_t nelx, size_t nely, double h, size_t nfine = 1); void map(const xt::xtensor& coor, const xt::xtensor& conn); friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular; }; // Mesh mappings namespace Map { // Return "FineLayer"-class responsible for generating a connectivity // Throws if conversion is not possible [[ deprecated ]] GooseFEM::Mesh::Quad4::FineLayer FineLayer( const xt::xtensor& coor, const xt::xtensor& conn); /** Refine a Regular mesh: subdivide elements in several smaller elements. */ class RefineRegular { public: RefineRegular() = default; /** Constructor. \param mesh the coarse mesh. \param nx for each coarse element: number of fine elements in x-direction. \param ny for each coarse element: number of fine elements in y-direction. */ RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny); /** For each coarse element: number of fine elements in x-direction. \return unsigned int (same as used in constructor) */ size_t nx() const; /** For each coarse element: number of fine elements in y-direction. \return unsigned int (same as used in constructor) */ size_t ny() const; /** Obtain the coarse mesh (copy of the mesh passed to the constructor). \return mesh */ GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const; /** Obtain the fine mesh. \return mesh */ GooseFEM::Mesh::Quad4::Regular getFineMesh() const; /** Get element-mapping: elements of the fine mesh per element of the coarse mesh. \return [nelem_coarse, nx() * ny()] */ xt::xtensor getMap() const; // map field [[ deprecated ]] xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per el [[ deprecated ]] xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per intpnt [[ deprecated ]] xt::xtensor mapToCoarse(const xt::xtensor& data) const; // tensor per intpnt /** Compute the mean of the quantity define on the fine mesh when mapped on the coarse mesh. \tparam T type of the data (e.g. ``double``). \tparam rank rank of the data. \param data the data [nelem_fine, ...] \return the average data of the coarse mesh [nelem_coarse, ...] */ template xt::xtensor meanToCoarse(const xt::xtensor& data) const; /** Compute the average of the quantity define on the fine mesh when mapped on the coarse mesh. \tparam T type of the data (e.g. ``double``). \tparam rank rank of the data. \tparam S type of the weights (e.g. ``double``). \param data the data [nelem_fine, ...] \param weights the weights [nelem_fine, ...] \return the average data of the coarse mesh [nelem_coarse, ...] */ template xt::xtensor averageToCoarse( const xt::xtensor& data, const xt::xtensor& weights) const; /** Map element quantities to the fine mesh. The mapping is a bit simplistic: no interpolation is involved. The mapping is such that:: ret[e_fine, ...] <- data[e_coarse, ...] \tparam T type of the data (e.g. ``double``). \tparam rank rank of the data. \param data the data. \return mapped data. */ template xt::xtensor mapToFine(const xt::xtensor& data) const; private: GooseFEM::Mesh::Quad4::Regular m_coarse; ///< the coarse mesh GooseFEM::Mesh::Quad4::Regular m_fine; ///< the fine mesh size_t m_nx; ///< see nx() size_t m_ny; ///< see ny() xt::xtensor m_coarse2fine; ///< see getMap() }; /** Map a FineLayer mesh to a Regular mesh. The element size of the Regular corresponds to the smallest elements of the FineLayer mesh (along the middle layer). */ class FineLayer2Regular { public: FineLayer2Regular() = default; /** Constructors. \param mesh The FineLayer mesh. */ FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh); /** Obtain the Regular mesh. \return mesh. */ GooseFEM::Mesh::Quad4::Regular getRegularMesh() const; /** Obtain the FineLayer mesh (copy of the mesh passed to the constructor). \return mesh. */ GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const; // elements of the Regular mesh per element of the FineLayer mesh // and the fraction by which the overlap is /** Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh. The number of Regular elements varies between elements of the FineLayer mesh. \return [nelem_finelayer, ?] */ std::vector> getMap() const; /** To overlap fraction for each item in the mapping in getMap(). \return [nelem_finelayer, ?] */ std::vector> getMapFraction() const; /** Map element quantities to Regular. The mapping is a bit simplistic: no interpolation is involved, the function just accounts the fraction of overlap between the FineLayer element and the Regular element. The mapping is such that:: ret[e_regular, ...] <- arg[e_finelayer, ...] \tparam T type of the data (e.g. ``double``). \tparam rank rank of the data. \param arg data. \return mapped data. */ template xt::xtensor mapToRegular(const xt::xtensor& arg) const; private: GooseFEM::Mesh::Quad4::FineLayer m_finelayer; ///< the FineLayer mesh to map GooseFEM::Mesh::Quad4::Regular m_regular; ///< the new Regular mesh to which to map std::vector> m_elem_regular; ///< see getMap() std::vector> m_frac_regular; ///< see getMapFraction() }; } // namespace Map } // namespace Quad4 } // namespace Mesh } // namespace GooseFEM #include "MeshQuad4.hpp" #endif diff --git a/include/GooseFEM/MeshQuad4.hpp b/include/GooseFEM/MeshQuad4.hpp index 45c6f38..9a1bf1a 100644 --- a/include/GooseFEM/MeshQuad4.hpp +++ b/include/GooseFEM/MeshQuad4.hpp @@ -1,1652 +1,1548 @@ /** Implementation of MeshQuad4.h \file MeshQuad4.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHQUAD4_HPP #define GOOSEFEM_MESHQUAD4_HPP #include "MeshQuad4.h" namespace GooseFEM { namespace Mesh { namespace Quad4 { -inline Regular::Regular(size_t nelx, size_t nely, double h) : m_h(h), m_nelx(nelx), m_nely(nely) +inline Regular::Regular(size_t nelx, size_t nely, double h) { - GOOSEFEM_ASSERT(m_nelx >= 1ul); - GOOSEFEM_ASSERT(m_nely >= 1ul); + m_h = h; + m_nelx = nelx; + m_nely = nely; + m_ndim = 2; + m_nne = 4; + + GOOSEFEM_ASSERT(m_nelx >= 1); + GOOSEFEM_ASSERT(m_nely >= 1); m_nnode = (m_nelx + 1) * (m_nely + 1); m_nelem = m_nelx * m_nely; } -inline size_t Regular::nelem() const -{ - return m_nelem; -} - -inline size_t Regular::nnode() const -{ - return m_nnode; -} - -inline size_t Regular::nne() const -{ - return m_nne; -} - -inline size_t Regular::ndim() const -{ - return m_ndim; -} - -inline size_t Regular::nelx() const -{ - return m_nelx; -} - -inline size_t Regular::nely() const -{ - return m_nely; -} - -inline double Regular::h() const -{ - return m_h; -} - inline ElementType Regular::getElementType() const { return ElementType::Quad4; } inline xt::xtensor Regular::coor() const { xt::xtensor ret = xt::empty({m_nnode, m_ndim}); xt::xtensor x = xt::linspace(0.0, m_h * static_cast(m_nelx), m_nelx + 1); xt::xtensor y = xt::linspace(0.0, m_h * static_cast(m_nely), m_nely + 1); size_t inode = 0; for (size_t iy = 0; iy < m_nely + 1; ++iy) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ++inode; } } return ret; } inline xt::xtensor Regular::conn() const { xt::xtensor ret = xt::empty({m_nelem, m_nne}); size_t ielem = 0; for (size_t iy = 0; iy < m_nely; ++iy) { for (size_t ix = 0; ix < m_nelx; ++ix) { ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix); ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1); ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix); ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1); ++ielem; } } return ret; } inline xt::xtensor Regular::nodesBottomEdge() const { return xt::arange(m_nelx + 1); } inline xt::xtensor Regular::nodesTopEdge() const { return xt::arange(m_nelx + 1) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1) + m_nelx; } inline xt::xtensor Regular::nodesBottomOpenEdge() const { return xt::arange(1, m_nelx); } inline xt::xtensor Regular::nodesTopOpenEdge() const { return xt::arange(1, m_nelx) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftOpenEdge() const { return xt::arange(1, m_nely) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightOpenEdge() const { return xt::arange(1, m_nely) * (m_nelx + 1) + m_nelx; } inline size_t Regular::nodesBottomLeftCorner() const { return 0; } inline size_t Regular::nodesBottomRightCorner() const { return m_nelx; } inline size_t Regular::nodesTopLeftCorner() const { return m_nely * (m_nelx + 1); } inline size_t Regular::nodesTopRightCorner() const { return m_nely * (m_nelx + 1) + m_nelx; } -inline size_t Regular::nodesLeftBottomCorner() const -{ - return nodesBottomLeftCorner(); -} - -inline size_t Regular::nodesLeftTopCorner() const -{ - return nodesTopLeftCorner(); -} - -inline size_t Regular::nodesRightBottomCorner() const -{ - return nodesBottomRightCorner(); -} - -inline size_t Regular::nodesRightTopCorner() const -{ - return nodesTopRightCorner(); -} - -inline xt::xtensor Regular::nodesPeriodic() const -{ - xt::xtensor bot = nodesBottomOpenEdge(); - xt::xtensor top = nodesTopOpenEdge(); - xt::xtensor lft = nodesLeftOpenEdge(); - xt::xtensor rgt = nodesRightOpenEdge(); - std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; - xt::xtensor ret = xt::empty(shape); - - ret(0, 0) = nodesBottomLeftCorner(); - ret(0, 1) = nodesBottomRightCorner(); - - ret(1, 0) = nodesBottomLeftCorner(); - ret(1, 1) = nodesTopRightCorner(); - - ret(2, 0) = nodesBottomLeftCorner(); - ret(2, 1) = nodesTopLeftCorner(); - - size_t i = 3; - - xt::view(ret, xt::range(i, i + bot.size()), 0) = bot; - xt::view(ret, xt::range(i, i + bot.size()), 1) = top; - - i += bot.size(); - - xt::view(ret, xt::range(i, i + lft.size()), 0) = lft; - xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt; - - return ret; -} - -inline size_t Regular::nodesOrigin() const -{ - return nodesBottomLeftCorner(); -} - -inline xt::xtensor Regular::dofs() const -{ - return GooseFEM::Mesh::dofs(m_nnode, m_ndim); -} - -inline xt::xtensor Regular::dofsPeriodic() const -{ - xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); - xt::xtensor nodePer = nodesPeriodic(); - xt::xtensor independent = xt::view(nodePer, xt::all(), 0); - xt::xtensor dependent = xt::view(nodePer, xt::all(), 1); - - for (size_t j = 0; j < m_ndim; ++j) { - xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j); - } - - return GooseFEM::Mesh::renumber(ret); -} - inline xt::xtensor Regular::elementgrid() const { return xt::arange(m_nelem).reshape({m_nely, m_nelx}); } inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine) { this->init(nelx, nely, h, nfine); } inline FineLayer::FineLayer(const xt::xtensor& coor, const xt::xtensor& conn) { this->map(coor, conn); } inline void FineLayer::init(size_t nelx, size_t nely, double h, size_t nfine) { GOOSEFEM_ASSERT(nelx >= 1ul); GOOSEFEM_ASSERT(nely >= 1ul); m_h = h; m_Lx = m_h * static_cast(nelx); // compute element size in y-direction (use symmetry, compute upper half) // temporary variables size_t nmin, ntot; xt::xtensor nhx = xt::ones({nely}); xt::xtensor nhy = xt::ones({nely}); xt::xtensor refine = -1 * xt::ones({nely}); // minimum height in y-direction (half of the height because of symmetry) if (nely % 2 == 0) { nmin = nely / 2; } else { nmin = (nely + 1) / 2; } // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half) if (nfine % 2 == 0) { nfine = nfine / 2 + 1; } else { nfine = (nfine + 1) / 2; } if (nfine < 1) { nfine = 1; } if (nfine > nmin) { nfine = nmin; } // loop over element layers in y-direction, try to coarsen using these rules: // (1) element size in y-direction <= distance to origin in y-direction // (2) element size in x-direction should fit the total number of elements in x-direction // (3) a certain number of layers have the minimum size "1" (are fine) for (size_t iy = nfine;;) { // initialize current size in y-direction if (iy == nfine) { ntot = nfine; } // check to stop if (iy >= nely || ntot >= nmin) { nely = iy; break; } // rules (1,2) satisfied: coarsen in x-direction if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) { refine(iy) = 0; nhy(iy) *= 2; auto vnhy = xt::view(nhy, xt::range(iy + 1, _)); auto vnhx = xt::view(nhx, xt::range(iy, _)); vnhy *= 3; vnhx *= 3; } // update the number of elements in y-direction ntot += nhy(iy); // proceed to next element layer in y-direction ++iy; // check to stop if (iy >= nely || ntot >= nmin) { nely = iy; break; } } // symmetrize, compute full information // allocate mesh constructor parameters m_nhx = xt::empty({nely * 2 - 1}); m_nhy = xt::empty({nely * 2 - 1}); m_refine = xt::empty({nely * 2 - 1}); m_nelx = xt::empty({nely * 2 - 1}); m_nnd = xt::empty({nely * 2}); m_startElem = xt::empty({nely * 2 - 1}); m_startNode = xt::empty({nely * 2}); // fill // - lower half for (size_t iy = 0; iy < nely; ++iy) { m_nhx(iy) = nhx(nely - iy - 1); m_nhy(iy) = nhy(nely - iy - 1); m_refine(iy) = refine(nely - iy - 1); } // - upper half for (size_t iy = 0; iy < nely - 1; ++iy) { m_nhx(iy + nely) = nhx(iy + 1); m_nhy(iy + nely) = nhy(iy + 1); m_refine(iy + nely) = refine(iy + 1); } // update size nely = m_nhx.size(); // compute the number of elements per element layer in y-direction for (size_t iy = 0; iy < nely; ++iy) { m_nelx(iy) = nelx / m_nhx(iy); } // compute the number of nodes per node layer in y-direction for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { m_nnd(iy) = m_nelx(iy) + 1; } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { m_nnd(iy + 1) = m_nelx(iy) + 1; } // compute mesh dimensions // initialize m_nnode = 0; m_nelem = 0; m_startNode(0) = 0; // loop over element layers (bottom -> middle, elements become finer) for (size_t i = 0; i < (nely - 1) / 2; ++i) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if (m_refine(i) == 0) { m_nnode += (3 * m_nelx(i) + 1); } else { m_nnode += (m_nelx(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_nelx(i)); } else { m_nelem += (m_nelx(i)); } // - store the starting node of the next layer m_startNode(i + 1) = m_nnode; } // loop over element layers (middle -> top, elements become coarser) for (size_t i = (nely - 1) / 2; i < nely; ++i) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if (m_refine(i) == 0) { m_nnode += (5 * m_nelx(i) + 1); } else { m_nnode += (m_nelx(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_nelx(i)); } else { m_nelem += (m_nelx(i)); } // - store the starting node of the next layer m_startNode(i + 1) = m_nnode; } // - add the top row of nodes m_nnode += m_nelx(nely - 1) + 1; } inline size_t FineLayer::nelem() const { return m_nelem; } inline size_t FineLayer::nnode() const { return m_nnode; } inline size_t FineLayer::nne() const { return m_nne; } inline size_t FineLayer::ndim() const { return m_ndim; } inline size_t FineLayer::nelx() const { return xt::amax(m_nelx)(); } inline size_t FineLayer::nely() const { return xt::sum(m_nhy)(); } inline double FineLayer::h() const { return m_h; } inline xt::xtensor FineLayer::elemrow_nhx() const { return m_nhx; } inline xt::xtensor FineLayer::elemrow_nhy() const { return m_nhy; } inline xt::xtensor FineLayer::elemrow_nelem() const { return m_nelx; } inline ElementType FineLayer::getElementType() const { return ElementType::Quad4; } inline xt::xtensor FineLayer::coor() const { // allocate output xt::xtensor ret = xt::empty({m_nnode, m_ndim}); // current node, number of element layers size_t inode = 0; size_t nely = static_cast(m_nhy.size()); // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) // - allocate xt::xtensor y = xt::empty({nely + 1}); // - initialize y(0) = 0.0; // - compute for (size_t iy = 1; iy < nely + 1; ++iy) { y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h; } // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes for (size_t iy = 0;; ++iy) { // get positions along the x- and z-axis xt::xtensor x = xt::linspace(0.0, m_Lx, m_nelx(iy) + 1); // add nodes of the bottom layer of this element for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ++inode; } // stop at middle layer if (iy == (nely - 1) / 2) { break; } // add extra nodes of the intermediate layer, for refinement in x-direction if (m_refine(iy) == 0) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t ix = 0; ix < m_nelx(iy); ++ix) { for (size_t j = 0; j < 2; ++j) { ret(inode, 0) = x(ix) + dx * static_cast(j + 1); ret(inode, 1) = y(iy) + dy; ++inode; } } } } // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { // get positions along the x- and z-axis xt::xtensor x = xt::linspace(0.0, m_Lx, m_nelx(iy) + 1); // add extra nodes of the intermediate layer, for refinement in x-direction if (m_refine(iy) == 0) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t ix = 0; ix < m_nelx(iy); ++ix) { for (size_t j = 0; j < 2; ++j) { ret(inode, 0) = x(ix) + dx * static_cast(j + 1); ret(inode, 1) = y(iy) + dy; ++inode; } } } // add nodes of the top layer of this element for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy + 1); ++inode; } } return ret; } inline xt::xtensor FineLayer::conn() const { // allocate output xt::xtensor ret = xt::empty({m_nelem, m_nne}); // current element, number of element layers, starting nodes of each node layer size_t ielem = 0; size_t nely = static_cast(m_nhy.size()); size_t bot, mid, top; // loop over all element layers for (size_t iy = 0; iy < nely; ++iy) { // - get: starting nodes of bottom(, middle) and top layer bot = m_startNode(iy); mid = m_startNode(iy) + m_nnd(iy); top = m_startNode(iy + 1); // - define connectivity: no coarsening/refinement if (m_refine(iy) == -1) { for (size_t ix = 0; ix < m_nelx(iy); ++ix) { ret(ielem, 0) = bot + (ix); ret(ielem, 1) = bot + (ix + 1); ret(ielem, 2) = top + (ix + 1); ret(ielem, 3) = top + (ix); ielem++; } } // - define connectivity: refinement along the x-direction (below the middle layer) else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) { for (size_t ix = 0; ix < m_nelx(iy); ++ix) { // -- bottom element ret(ielem, 0) = bot + (ix); ret(ielem, 1) = bot + (ix + 1); ret(ielem, 2) = mid + (2 * ix + 1); ret(ielem, 3) = mid + (2 * ix); ielem++; // -- top-right element ret(ielem, 0) = bot + (ix + 1); ret(ielem, 1) = top + (3 * ix + 3); ret(ielem, 2) = top + (3 * ix + 2); ret(ielem, 3) = mid + (2 * ix + 1); ielem++; // -- top-center element ret(ielem, 0) = mid + (2 * ix); ret(ielem, 1) = mid + (2 * ix + 1); ret(ielem, 2) = top + (3 * ix + 2); ret(ielem, 3) = top + (3 * ix + 1); ielem++; // -- top-left element ret(ielem, 0) = bot + (ix); ret(ielem, 1) = mid + (2 * ix); ret(ielem, 2) = top + (3 * ix + 1); ret(ielem, 3) = top + (3 * ix); ielem++; } } // - define connectivity: coarsening along the x-direction (above the middle layer) else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) { for (size_t ix = 0; ix < m_nelx(iy); ++ix) { // -- lower-left element ret(ielem, 0) = bot + (3 * ix); ret(ielem, 1) = bot + (3 * ix + 1); ret(ielem, 2) = mid + (2 * ix); ret(ielem, 3) = top + (ix); ielem++; // -- lower-center element ret(ielem, 0) = bot + (3 * ix + 1); ret(ielem, 1) = bot + (3 * ix + 2); ret(ielem, 2) = mid + (2 * ix + 1); ret(ielem, 3) = mid + (2 * ix); ielem++; // -- lower-right element ret(ielem, 0) = bot + (3 * ix + 2); ret(ielem, 1) = bot + (3 * ix + 3); ret(ielem, 2) = top + (ix + 1); ret(ielem, 3) = mid + (2 * ix + 1); ielem++; // -- upper element ret(ielem, 0) = mid + (2 * ix); ret(ielem, 1) = mid + (2 * ix + 1); ret(ielem, 2) = top + (ix + 1); ret(ielem, 3) = top + (ix); ielem++; } } } return ret; } inline xt::xtensor FineLayer::elementsMiddleLayer() const { size_t nely = m_nhy.size(); size_t iy = (nely - 1) / 2; return m_startElem(iy) + xt::arange(m_nelx(iy)); } inline xt::xtensor FineLayer::elementsLayer(size_t iy) const { GOOSEFEM_ASSERT(iy < m_nelx.size()); size_t n = m_nelx(iy); if (m_refine(iy) != -1) { n *= 4; } return m_startElem(iy) + xt::arange(n); } inline xt::xtensor FineLayer::elementgrid_ravel( std::vector start_stop_rows, std::vector start_stop_cols) const { GOOSEFEM_ASSERT(start_stop_rows.size() == 0 || start_stop_rows.size() == 2); GOOSEFEM_ASSERT(start_stop_cols.size() == 0 || start_stop_cols.size() == 2); std::array rows; std::array cols; if (start_stop_rows.size() == 2) { std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin()); GOOSEFEM_ASSERT(rows[0] <= this->nely()); GOOSEFEM_ASSERT(rows[1] <= this->nely()); } else { rows[0] = 0; rows[1] = this->nely(); } if (start_stop_cols.size() == 2) { std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin()); GOOSEFEM_ASSERT(cols[0] <= this->nelx()); GOOSEFEM_ASSERT(cols[1] <= this->nelx()); } else { cols[0] = 0; cols[1] = this->nelx(); } if (rows[0] == rows[1] || cols[0] == cols[1]) { xt::xtensor ret = xt::empty({0}); return ret; } // Compute dimensions auto H = xt::cumsum(m_nhy); size_t yl = 0; if (rows[0] > 0) { yl = xt::argmax(H > rows[0])(); } size_t yu = xt::argmax(H >= rows[1])(); size_t hx = std::max(m_nhx(yl), m_nhx(yu)); size_t xl = (cols[0] - cols[0] % hx) / hx; size_t xu = (cols[1] - cols[1] % hx) / hx; // Allocate output size_t N = 0; for (size_t i = yl; i <= yu; ++i) { // no refinement size_t n = (xu - xl) * hx / m_nhx(i); // refinement if (m_refine(i) != -1) { n *= 4; } N += n; } xt::xtensor ret = xt::empty({N}); // Write output N = 0; for (size_t i = yl; i <= yu; ++i) { // no refinement size_t n = (xu - xl) * hx / m_nhx(i); size_t h = hx; // refinement if (m_refine(i) != -1) { n *= 4; h *= 4; } xt::xtensor e = m_startElem(i) + xl * h / m_nhx(i) + xt::arange(n); xt::view(ret, xt::range(N, N + n)) = e; N += n; } return ret; } inline xt::xtensor FineLayer::elementgrid_around_ravel( size_t e, size_t size, bool periodic) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_nelx(iy); GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2); auto H = xt::cumsum(m_nhy); if (2 * size >= H(H.size() - 1)) { return xt::arange(this->nelem()); } size_t hy = H(iy); size_t l = xt::argmax(H > (hy - size - 1))(); size_t u = xt::argmax(H >= (hy + size))(); size_t lh = 0; if (l > 0) { lh = H(l - 1); } size_t uh = H(u); size_t step = xt::amax(m_nhx)(); size_t relx = (e - m_startElem(iy)) % step; size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx; size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step; size_t dx = m_nhx(u); size_t xl = 0; size_t xu = nel; if (mid >= size) { xl = mid - size; } if (mid + size < nel) { xu = mid + size + 1; } xl = xl - xl % dx; xu = xu - xu % dx; if (mid - xl < size) { if (xl < dx) { xl = 0; } else { xl -= dx; } } if (xu - mid < size) { if (xu > nel - dx) { xu = nel; } else { xu += dx; } } auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu}); auto map = this->roll(nroll); return xt::view(map, xt::keep(ret)); } inline xt::xtensor FineLayer::elementgrid_leftright( size_t e, size_t left, size_t right, bool periodic) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_nelx(iy); GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2); size_t step = xt::amax(m_nhx)(); size_t relx = (e - m_startElem(iy)) % step; size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx; size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step; size_t dx = m_nhx(iy); size_t xl = 0; size_t xu = nel; if (mid >= left) { xl = mid - left; } if (mid + right < nel) { xu = mid + right + 1; } xl = xl - xl % dx; xu = xu - xu % dx; if (mid - xl < left) { if (xl < dx) { xl = 0; } else { xl -= dx; } } if (xu - mid < right) { if (xu > nel - dx) { xu = nel; } else { xu += dx; } } auto H = xt::cumsum(m_nhy); size_t yl = 0; if (iy > 0) { yl = H(iy - 1); } auto ret = this->elementgrid_ravel({yl, H(iy)}, {xl, xu}); auto map = this->roll(nroll); return xt::view(map, xt::keep(ret)); } inline xt::xtensor FineLayer::nodesBottomEdge() const { return m_startNode(0) + xt::arange(m_nelx(0) + 1); } inline xt::xtensor FineLayer::nodesTopEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(m_nelx(nely - 1) + 1); } inline xt::xtensor FineLayer::nodesLeftEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely + 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j)); xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1)); return ret; } inline xt::xtensor FineLayer::nodesRightEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely + 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j)) + xt::view(m_nelx, xt::range(i, j)); xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1)) + xt::view(m_nelx, xt::range(k, l)); return ret; } inline xt::xtensor FineLayer::nodesBottomOpenEdge() const { return m_startNode(0) + xt::arange(1, m_nelx(0)); } inline xt::xtensor FineLayer::nodesTopOpenEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(1, m_nelx(nely - 1)); } inline xt::xtensor FineLayer::nodesLeftOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely - 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j)); xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l)); return ret; } inline xt::xtensor FineLayer::nodesRightOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely - 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j)) + xt::view(m_nelx, xt::range(i + 1, j)); xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l)) + xt::view(m_nelx, xt::range(k, l - 1)); return ret; } inline size_t FineLayer::nodesBottomLeftCorner() const { return m_startNode(0); } inline size_t FineLayer::nodesBottomRightCorner() const { return m_startNode(0) + m_nelx(0); } inline size_t FineLayer::nodesTopLeftCorner() const { size_t nely = m_nhy.size(); return m_startNode(nely); } inline size_t FineLayer::nodesTopRightCorner() const { size_t nely = m_nhy.size(); return m_startNode(nely) + m_nelx(nely - 1); } inline size_t FineLayer::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } inline size_t FineLayer::nodesRightBottomCorner() const { return nodesBottomRightCorner(); } inline size_t FineLayer::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } inline size_t FineLayer::nodesRightTopCorner() const { return nodesTopRightCorner(); } inline xt::xtensor FineLayer::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(shape); ret(0, 0) = nodesBottomLeftCorner(); ret(0, 1) = nodesBottomRightCorner(); ret(1, 0) = nodesBottomLeftCorner(); ret(1, 1) = nodesTopRightCorner(); ret(2, 0) = nodesBottomLeftCorner(); ret(2, 1) = nodesTopLeftCorner(); size_t i = 3; xt::view(ret, xt::range(i, i + bot.size()), 0) = bot; xt::view(ret, xt::range(i, i + bot.size()), 1) = top; i += bot.size(); xt::view(ret, xt::range(i, i + lft.size()), 0) = lft; xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt; return ret; } inline size_t FineLayer::nodesOrigin() const { return nodesBottomLeftCorner(); } inline xt::xtensor FineLayer::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor FineLayer::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor dependent = xt::view(nodePer, xt::all(), 1); for (size_t j = 0; j < m_ndim; ++j) { xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j); } return GooseFEM::Mesh::renumber(ret); } inline xt::xtensor FineLayer::roll(size_t n) { auto conn = this->conn(); size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_nelem}); // loop over all element layers for (size_t iy = 0; iy < nely; ++iy) { // no refinement size_t shift = n * (m_nelx(iy) / m_nelx(0)); size_t nel = m_nelx(iy); // refinement if (m_refine(iy) != -1) { shift = n * (m_nelx(iy) / m_nelx(0)) * 4; nel = m_nelx(iy) * 4; } // element numbers of the layer, and roll them auto e = m_startElem(iy) + xt::arange(nel); xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift); } return ret; } inline void FineLayer::map(const xt::xtensor& coor, const xt::xtensor& conn) { GOOSEFEM_ASSERT(coor.shape(1) == 2); GOOSEFEM_ASSERT(conn.shape(1) == 4); GOOSEFEM_ASSERT(conn.shape(0) > 0); GOOSEFEM_ASSERT(coor.shape(0) >= 4); GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); if (conn.shape(0) == 1) { this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0)); GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn))); GOOSEFEM_CHECK(xt::allclose(this->coor(), coor)); return; } // Identify the middle layer size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2; size_t eleft = emid; size_t eright = emid; while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) { eleft--; } while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) { eright++; } GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0)); // Get element sizes along the middle layer auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0); auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1); auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 2); auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0); auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1); auto hx = xt::amin(dx)(); auto hy = xt::amin(dy)(); GOOSEFEM_CHECK(xt::allclose(hx, hy)); GOOSEFEM_CHECK(xt::allclose(dx, hx)); GOOSEFEM_CHECK(xt::allclose(dy, hy)); // Extract shape and initialise size_t nelx = eright - eleft + 1; size_t nely = static_cast((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx); this->init(nelx, nely, hx); GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn))); GOOSEFEM_CHECK(xt::allclose(this->coor(), coor)); GOOSEFEM_CHECK(xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange(nelx)))); } namespace Map { inline RefineRegular::RefineRegular( const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny) : m_coarse(mesh), m_nx(nx), m_ny(ny) { m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h()); xt::xtensor elmat_coarse = m_coarse.elementgrid(); xt::xtensor elmat_fine = m_fine.elementgrid(); m_coarse2fine = xt::empty({m_coarse.nelem(), nx * ny}); for (size_t i = 0; i < elmat_coarse.shape(0); ++i) { for (size_t j = 0; j < elmat_coarse.shape(1); ++j) { xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view( elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx))); } } } inline size_t RefineRegular::nx() const { return m_nx; } inline size_t RefineRegular::ny() const { return m_ny; } inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getCoarseMesh() const { return m_coarse; } inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getFineMesh() const { return m_fine; } inline xt::xtensor RefineRegular::getMap() const { return m_coarse2fine; } template inline xt::xtensor RefineRegular::meanToCoarse(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_coarse2fine.shape(0); xt::xtensor ret = xt::empty(shape); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); auto d = xt::view(data, xt::keep(e)); xt::view(ret, i) = xt::mean(d, 0); } return ret; } template inline xt::xtensor RefineRegular::averageToCoarse( const xt::xtensor& data, const xt::xtensor& weights) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_coarse2fine.shape(0); xt::xtensor ret = xt::empty(shape); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); xt::xtensor d = xt::view(data, xt::keep(e)); xt::xtensor w = xt::view(weights, xt::keep(e)); xt::view(ret, i) = xt::average(d, w, {0}); } return ret; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& data) const { GOOSEFEM_WARNING("mapToCoarse is deprecated, use meanToCoarse or averageToCoarse"); GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); size_t m = m_coarse2fine.shape(0); size_t n = m_coarse2fine.shape(1); xt::xtensor ret = xt::empty({m, n}); for (size_t i = 0; i < m; ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); xt::view(ret, i) = xt::view(data, xt::keep(e)); } return ret; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& data) const { GOOSEFEM_WARNING("mapToCoarse is deprecated, use meanToCoarse or averageToCoarse"); GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); size_t m = m_coarse2fine.shape(0); size_t n = m_coarse2fine.shape(1); size_t N = data.shape(1); xt::xtensor ret = xt::empty({m, n * N}); for (size_t i = 0; i < m; ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); for (size_t q = 0; q < data.shape(1); ++q) { xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q); } } return ret; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& data) const { GOOSEFEM_WARNING("mapToCoarse is deprecated, use meanToCoarse or averageToCoarse"); GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); size_t m = m_coarse2fine.shape(0); size_t n = m_coarse2fine.shape(1); size_t N = data.shape(1); xt::xtensor ret = xt::empty({m, n * N, data.shape(2), data.shape(3)}); for (size_t i = 0; i < m; ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); for (size_t q = 0; q < data.shape(1); ++q) { xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q); } } return ret; } template inline xt::xtensor RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_coarse2fine.size(); xt::xtensor ret = xt::empty(shape); for (size_t e = 0; e < m_coarse2fine.shape(0); ++e) { for (size_t i = 0; i < m_coarse2fine.shape(1); ++i) { xt::view(ret, m_coarse2fine(e, i)) = xt::view(data, e); } } return ret; } inline FineLayer2Regular::FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh) : m_finelayer(mesh) { // ------------ // Regular-mesh // ------------ m_regular = GooseFEM::Mesh::Quad4::Regular( xt::amax(m_finelayer.m_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h); // ------- // mapping // ------- // allocate mapping m_elem_regular.resize(m_finelayer.m_nelem); m_frac_regular.resize(m_finelayer.m_nelem); // alias xt::xtensor nhx = m_finelayer.m_nhx; xt::xtensor nhy = m_finelayer.m_nhy; xt::xtensor nelx = m_finelayer.m_nelx; xt::xtensor start = m_finelayer.m_startElem; // 'matrix' of element numbers of the Regular-mesh xt::xtensor elementgrid = m_regular.elementgrid(); // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh xt::xtensor cum_nhy = xt::concatenate(xt::xtuple(xt::zeros({1}), xt::cumsum(nhy))); // number of element layers in y-direction of the FineLayer-mesh size_t nely = nhy.size(); // loop over layers of the FineLayer-mesh for (size_t iy = 0; iy < nely; ++iy) { // element numbers of the Regular-mesh along this layer of the FineLayer-mesh auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all()); // no coarsening/refinement // ------------------------ if (m_finelayer.m_refine(iy) == -1) { // element numbers of the FineLayer-mesh along this layer xt::xtensor el_old = start(iy) + xt::arange(nelx(iy)); // loop along this layer of the FineLayer-mesh for (size_t ix = 0; ix < nelx(iy); ++ix) { // get the element numbers of the Regular-mesh for this element of the // FineLayer-mesh auto block = xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy))); // write to mapping for (auto& i : block) { m_elem_regular[el_old(ix)].push_back(i); m_frac_regular[el_old(ix)].push_back(1.0); } } } // refinement along the x-direction (below the middle layer) // --------------------------------------------------------- else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) { // element numbers of the FineLayer-mesh along this layer // rows: coarse block, columns element numbers per block xt::xtensor el_old = start(iy) + xt::arange(nelx(iy) * 4ul).reshape({-1, 4}); // loop along this layer of the FineLayer-mesh for (size_t ix = 0; ix < nelx(iy); ++ix) { // get the element numbers of the Regular-mesh for this block of the FineLayer-mesh auto block = xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy))); // bottom: wide-to-narrow { for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, j, xt::range(j, nhx(iy) - j)); m_elem_regular[el_old(ix, 0)].push_back(e(0)); m_frac_regular[el_old(ix, 0)].push_back(0.5); for (size_t k = 1; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 0)].push_back(e(k)); m_frac_regular[el_old(ix, 0)].push_back(1.0); } m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 0)].push_back(0.5); } } // top: regular small element { auto e = xt::view( block, xt::range(nhy(iy) / 2, nhy(iy)), xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 2)].push_back(i); m_frac_regular[el_old(ix, 2)].push_back(1.0); } } // left { // left-bottom: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, j, xt::range(0, j + 1)); for (size_t k = 0; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 3)].push_back(e(k)); m_frac_regular[el_old(ix, 3)].push_back(1.0); } m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 3)].push_back(0.5); } // left-top: regular { auto e = xt::view( block, xt::range(nhy(iy) / 2, nhy(iy)), xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 3)].push_back(i); m_frac_regular[el_old(ix, 3)].push_back(1.0); } } } // right { // right-bottom: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy))); m_elem_regular[el_old(ix, 1)].push_back(e(0)); m_frac_regular[el_old(ix, 1)].push_back(0.5); for (size_t k = 1; k < e.size(); ++k) { m_elem_regular[el_old(ix, 1)].push_back(e(k)); m_frac_regular[el_old(ix, 1)].push_back(1.0); } } // right-top: regular { auto e = xt::view( block, xt::range(nhy(iy) / 2, nhy(iy)), xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 1)].push_back(i); m_frac_regular[el_old(ix, 1)].push_back(1.0); } } } } } // coarsening along the x-direction (above the middle layer) else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) { // element numbers of the FineLayer-mesh along this layer // rows: coarse block, columns element numbers per block xt::xtensor el_old = start(iy) + xt::arange(nelx(iy) * 4ul).reshape({-1, 4}); // loop along this layer of the FineLayer-mesh for (size_t ix = 0; ix < nelx(iy); ++ix) { // get the element numbers of the Regular-mesh for this block of the FineLayer-mesh auto block = xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy))); // top: narrow-to-wide { for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view( block, nhy(iy) / 2 + j, xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1)); m_elem_regular[el_old(ix, 3)].push_back(e(0)); m_frac_regular[el_old(ix, 3)].push_back(0.5); for (size_t k = 1; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 3)].push_back(e(k)); m_frac_regular[el_old(ix, 3)].push_back(1.0); } m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 3)].push_back(0.5); } } // bottom: regular small element { auto e = xt::view( block, xt::range(0, nhy(iy) / 2), xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 1)].push_back(i); m_frac_regular[el_old(ix, 1)].push_back(1.0); } } // left { // left-bottom: regular { auto e = xt::view( block, xt::range(0, nhy(iy) / 2), xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 0)].push_back(i); m_frac_regular[el_old(ix, 0)].push_back(1.0); } } // left-top: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j)); for (size_t k = 0; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 0)].push_back(e(k)); m_frac_regular[el_old(ix, 0)].push_back(1.0); } m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 0)].push_back(0.5); } } // right { // right-bottom: regular { auto e = xt::view( block, xt::range(0, nhy(iy) / 2), xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 2)].push_back(i); m_frac_regular[el_old(ix, 2)].push_back(1.0); } } // right-top: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view( block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy))); m_elem_regular[el_old(ix, 2)].push_back(e(0)); m_frac_regular[el_old(ix, 2)].push_back(0.5); for (size_t k = 1; k < e.size(); ++k) { m_elem_regular[el_old(ix, 2)].push_back(e(k)); m_frac_regular[el_old(ix, 2)].push_back(1.0); } } } } } } } inline GooseFEM::Mesh::Quad4::Regular FineLayer2Regular::getRegularMesh() const { return m_regular; } inline GooseFEM::Mesh::Quad4::FineLayer FineLayer2Regular::getFineLayerMesh() const { return m_finelayer; } inline std::vector> FineLayer2Regular::getMap() const { return m_elem_regular; } inline std::vector> FineLayer2Regular::getMapFraction() const { return m_frac_regular; } template inline xt::xtensor FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_regular.nelem(); xt::xtensor ret = xt::zeros(shape); for (size_t e = 0; e < m_finelayer.nelem(); ++e) { for (size_t i = 0; i < m_elem_regular[e].size(); ++i) { xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e); } } return ret; } } // namespace Map } // namespace Quad4 } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/MeshTri3.h b/include/GooseFEM/MeshTri3.h index dfb6387..142c6af 100644 --- a/include/GooseFEM/MeshTri3.h +++ b/include/GooseFEM/MeshTri3.h @@ -1,98 +1,74 @@ /** Generate simple meshes of 3-noded triangular elements in 2d (GooseFEM::Mesh::ElementType::Tri3). \file MeshTri3.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHTRI3_H #define GOOSEFEM_MESHTRI3_H #include "config.h" +#include "Mesh.h" namespace GooseFEM { namespace Mesh { + +/** +Simple meshes of triangular elements of type ElementType::Tri3. +*/ namespace Tri3 { -class Regular { +/** +Regular grid of squares, with each square cut into two triangular elements. +*/ +class Regular : public RegularBase2d { public: - Regular(size_t nelx, size_t nely, double h = 1.); - - // size - size_t nelem() const; // number of elements - size_t nnode() const; // number of nodes - size_t nne() const; // number of nodes-per-element - size_t ndim() const; // number of dimensions - - // type - ElementType getElementType() const; - // mesh - xt::xtensor coor() const; // nodal positions [nnode, ndim] - xt::xtensor conn() const; // connectivity [nelem, nne] + /** + Constructor. - // boundary nodes: edges - xt::xtensor nodesBottomEdge() const; - xt::xtensor nodesTopEdge() const; - xt::xtensor nodesLeftEdge() const; - xt::xtensor nodesRightEdge() const; - - // boundary nodes: edges, without corners - xt::xtensor nodesBottomOpenEdge() const; - xt::xtensor nodesTopOpenEdge() const; - xt::xtensor nodesLeftOpenEdge() const; - xt::xtensor nodesRightOpenEdge() const; - - // boundary nodes: corners (including aliases) - size_t nodesBottomLeftCorner() const; - size_t nodesBottomRightCorner() const; - size_t nodesTopLeftCorner() const; - size_t nodesTopRightCorner() const; - size_t nodesLeftBottomCorner() const; - size_t nodesLeftTopCorner() const; - size_t nodesRightBottomCorner() const; - size_t nodesRightTopCorner() const; - - // DOF-numbers for each component of each node (sequential) - xt::xtensor dofs() const; - - // DOF-numbers for the case that the periodicity if fully eliminated - xt::xtensor dofsPeriodic() const; - - // periodic node pairs [:,2]: (independent, dependent) - xt::xtensor nodesPeriodic() const; - - // front-bottom-left node, used as reference for periodicity - size_t nodesOrigin() const; + \param nelx Number of elements in x-direction. + \param nely Number of elements in y-direction. + \param h Edge-size (of the squares, and thus of two of three element-edges). + */ + Regular(size_t nelx, size_t nely, double h = 1.); -private: - double m_h; // elementary element edge-size (in all directions) - size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") - size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") - size_t m_nelem; // number of elements - size_t m_nnode; // number of nodes - static const size_t m_nne = 3; // number of nodes-per-element - static const size_t m_ndim = 2; // number of dimensions + ElementType getElementType() const override; + xt::xtensor coor() const override; + xt::xtensor conn() const override; + xt::xtensor nodesBottomEdge() const override; + xt::xtensor nodesTopEdge() const override; + xt::xtensor nodesLeftEdge() const override; + xt::xtensor nodesRightEdge() const override; + xt::xtensor nodesBottomOpenEdge() const override; + xt::xtensor nodesTopOpenEdge() const override; + xt::xtensor nodesLeftOpenEdge() const override; + xt::xtensor nodesRightOpenEdge() const override; + size_t nodesBottomLeftCorner() const override; + size_t nodesBottomRightCorner() const override; + size_t nodesTopLeftCorner() const override; + size_t nodesTopRightCorner() const override; }; // read / set the orientation (-1/+1) of all triangles inline xt::xtensor getOrientation( const xt::xtensor& coor, const xt::xtensor& conn); inline xt::xtensor setOrientation( const xt::xtensor& coor, const xt::xtensor& conn, int orientation = -1); inline xt::xtensor setOrientation( const xt::xtensor& coor, const xt::xtensor& conn, const xt::xtensor& current, // (output of "getOrientation") int orientation = -1); } // namespace Tri3 } // namespace Mesh } // namespace GooseFEM #include "MeshTri3.hpp" #endif diff --git a/include/GooseFEM/MeshTri3.hpp b/include/GooseFEM/MeshTri3.hpp index bf859fb..a19b556 100644 --- a/include/GooseFEM/MeshTri3.hpp +++ b/include/GooseFEM/MeshTri3.hpp @@ -1,356 +1,261 @@ /** Implementation of MeshTri3.h \file MeshTri3.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHTRI3_HPP #define GOOSEFEM_MESHTRI3_HPP #include "MeshTri3.h" namespace GooseFEM { namespace Mesh { namespace Tri3 { -inline Regular::Regular(size_t nelx, size_t nely, double h) : m_h(h), m_nelx(nelx), m_nely(nely) +inline Regular::Regular(size_t nelx, size_t nely, double h) { - GOOSEFEM_ASSERT(m_nelx >= 1ul); - GOOSEFEM_ASSERT(m_nely >= 1ul); + m_h = h; + m_nelx = nelx; + m_nely = nely; + m_ndim = 2; + m_nne = 3; + + GOOSEFEM_ASSERT(m_nelx >= 1); + GOOSEFEM_ASSERT(m_nely >= 1); m_nnode = (m_nelx + 1) * (m_nely + 1); m_nelem = m_nelx * m_nely * 2; } -inline size_t Regular::nelem() const -{ - return m_nelem; -} - -inline size_t Regular::nnode() const -{ - return m_nnode; -} - -inline size_t Regular::nne() const -{ - return m_nne; -} - -inline size_t Regular::ndim() const -{ - return m_ndim; -} - inline ElementType Regular::getElementType() const { return ElementType::Tri3; } inline xt::xtensor Regular::coor() const { xt::xtensor ret = xt::empty({m_nnode, m_ndim}); xt::xtensor x = xt::linspace(0.0, m_h * static_cast(m_nelx), m_nelx + 1); xt::xtensor y = xt::linspace(0.0, m_h * static_cast(m_nely), m_nely + 1); size_t inode = 0; for (size_t iy = 0; iy < m_nely + 1; ++iy) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ++inode; } } return ret; } inline xt::xtensor Regular::conn() const { xt::xtensor ret = xt::empty({m_nelem, m_nne}); size_t ielem = 0; for (size_t iy = 0; iy < m_nely; ++iy) { for (size_t ix = 0; ix < m_nelx; ++ix) { ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix); ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1); ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix); ++ielem; ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix + 1); ret(ielem, 1) = (iy + 1) * (m_nelx + 1) + (ix + 1); ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix); ++ielem; } } return ret; } inline xt::xtensor Regular::nodesBottomEdge() const { xt::xtensor ret = xt::empty({m_nelx + 1}); for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(ix) = ix; } return ret; } inline xt::xtensor Regular::nodesTopEdge() const { xt::xtensor ret = xt::empty({m_nelx + 1}); for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(ix) = ix + m_nely * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesLeftEdge() const { xt::xtensor ret = xt::empty({m_nely + 1}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iy) = iy * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesRightEdge() const { xt::xtensor ret = xt::empty({m_nely + 1}); for (size_t iy = 0; iy < m_nely + 1; ++iy) { ret(iy) = iy * (m_nelx + 1) + m_nelx; } return ret; } inline xt::xtensor Regular::nodesBottomOpenEdge() const { xt::xtensor ret = xt::empty({m_nelx - 1}); for (size_t ix = 1; ix < m_nelx; ++ix) { ret(ix - 1) = ix; } return ret; } inline xt::xtensor Regular::nodesTopOpenEdge() const { xt::xtensor ret = xt::empty({m_nelx - 1}); for (size_t ix = 1; ix < m_nelx; ++ix) { ret(ix - 1) = ix + m_nely * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesLeftOpenEdge() const { xt::xtensor ret = xt::empty({m_nely - 1}); for (size_t iy = 1; iy < m_nely; ++iy) { ret(iy - 1) = iy * (m_nelx + 1); } return ret; } inline xt::xtensor Regular::nodesRightOpenEdge() const { xt::xtensor ret = xt::empty({m_nely - 1}); for (size_t iy = 1; iy < m_nely; ++iy) { ret(iy - 1) = iy * (m_nelx + 1) + m_nelx; } return ret; } inline size_t Regular::nodesBottomLeftCorner() const { return 0; } inline size_t Regular::nodesBottomRightCorner() const { return m_nelx; } inline size_t Regular::nodesTopLeftCorner() const { return m_nely * (m_nelx + 1); } inline size_t Regular::nodesTopRightCorner() const { return m_nely * (m_nelx + 1) + m_nelx; } -inline size_t Regular::nodesLeftBottomCorner() const -{ - return nodesBottomLeftCorner(); -} - -inline size_t Regular::nodesLeftTopCorner() const -{ - return nodesTopLeftCorner(); -} - -inline size_t Regular::nodesRightBottomCorner() const -{ - return nodesBottomRightCorner(); -} - -inline size_t Regular::nodesRightTopCorner() const -{ - return nodesTopRightCorner(); -} - -inline xt::xtensor Regular::nodesPeriodic() const -{ - xt::xtensor bot = nodesBottomOpenEdge(); - xt::xtensor top = nodesTopOpenEdge(); - xt::xtensor lft = nodesLeftOpenEdge(); - xt::xtensor rgt = nodesRightOpenEdge(); - - size_t tedge = bot.size() + lft.size(); - size_t tnode = 3; - xt::xtensor ret = xt::empty({tedge + tnode, std::size_t(2)}); - - size_t i = 0; - - ret(i, 0) = nodesBottomLeftCorner(); - ret(i, 1) = nodesBottomRightCorner(); - ++i; - ret(i, 0) = nodesBottomLeftCorner(); - ret(i, 1) = nodesTopRightCorner(); - ++i; - ret(i, 0) = nodesBottomLeftCorner(); - ret(i, 1) = nodesTopLeftCorner(); - ++i; - - for (size_t j = 0; j < bot.size(); ++j) { - ret(i, 0) = bot(j); - ret(i, 1) = top(j); - ++i; - } - for (size_t j = 0; j < lft.size(); ++j) { - ret(i, 0) = lft(j); - ret(i, 1) = rgt(j); - ++i; - } - - return ret; -} - -inline size_t Regular::nodesOrigin() const -{ - return nodesBottomLeftCorner(); -} - -inline xt::xtensor Regular::dofs() const -{ - return GooseFEM::Mesh::dofs(m_nnode, m_ndim); -} - -inline xt::xtensor Regular::dofsPeriodic() const -{ - xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); - xt::xtensor nodePer = nodesPeriodic(); - - for (size_t i = 0; i < nodePer.shape(0); ++i) { - for (size_t j = 0; j < m_ndim; ++j) { - ret(nodePer(i, 1), j) = ret(nodePer(i, 0), j); - } - } - - return GooseFEM::Mesh::renumber(ret); -} - inline xt::xtensor getOrientation(const xt::xtensor& coor, const xt::xtensor& conn) { GOOSEFEM_ASSERT(conn.shape(1) == 3ul); GOOSEFEM_ASSERT(coor.shape(1) == 2ul); double k; size_t nelem = conn.shape(0); xt::xtensor ret = xt::empty({nelem}); for (size_t ielem = 0; ielem < nelem; ++ielem) { auto v1 = xt::view(coor, conn(ielem, 0), xt::all()) - xt::view(coor, conn(ielem, 1), xt::all()); auto v2 = xt::view(coor, conn(ielem, 2), xt::all()) - xt::view(coor, conn(ielem, 1), xt::all()); k = v1(0) * v2(1) - v2(0) * v1(1); if (k < 0) { ret(ielem) = -1; } else { ret(ielem) = +1; } } return ret; } inline xt::xtensor setOrientation( const xt::xtensor& coor, const xt::xtensor& conn, int orientation) { GOOSEFEM_ASSERT(conn.shape(1) == 3ul); GOOSEFEM_ASSERT(coor.shape(1) == 2ul); GOOSEFEM_ASSERT(orientation == -1 || orientation == +1); xt::xtensor val = getOrientation(coor, conn); return setOrientation(coor, conn, val, orientation); } inline xt::xtensor setOrientation( const xt::xtensor& coor, const xt::xtensor& conn, const xt::xtensor& val, int orientation) { GOOSEFEM_ASSERT(conn.shape(1) == 3ul); GOOSEFEM_ASSERT(coor.shape(1) == 2ul); GOOSEFEM_ASSERT(conn.shape(0) == val.size()); GOOSEFEM_ASSERT(orientation == -1 || orientation == +1); UNUSED(coor); size_t nelem = conn.shape(0); xt::xtensor ret = conn; for (size_t ielem = 0; ielem < nelem; ++ielem) { if ((orientation == -1 && val(ielem) > 0) || (orientation == +1 && val(ielem) < 0)) { std::swap(ret(ielem, 2), ret(ielem, 1)); } } return ret; } } // namespace Tri3 } // namespace Mesh } // namespace GooseFEM #endif