diff --git a/include/GooseFEM/MatrixDiagonal.h b/include/GooseFEM/MatrixDiagonal.h index 3dea2b3..7858f9f 100644 --- a/include/GooseFEM/MatrixDiagonal.h +++ b/include/GooseFEM/MatrixDiagonal.h @@ -1,177 +1,180 @@ /** Diagonal matrix. \file MatrixDiagonal.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MATRIXDIAGONAL_H #define GOOSEFEM_MATRIXDIAGONAL_H #include "Element.h" #include "config.h" namespace GooseFEM { /** Diagonal matrix. See Vector() for bookkeeping definitions. */ class MatrixDiagonal { public: MatrixDiagonal() = default; virtual ~MatrixDiagonal() = default; /** Constructor. + \tparam C e.g. `xt::xtensor` + \tparam D e.g. `xt::xtensor` \param conn connectivity [#nelem, #nne]. \param dofs DOFs per node [#nnode, #ndim]. */ - MatrixDiagonal(const xt::xtensor& conn, const xt::xtensor& dofs); + template + MatrixDiagonal(const C& conn, const D& dofs); /** \return Number of elements. */ size_t nelem() const; /** \return Number of nodes per element. */ size_t nne() const; /** \return Number of nodes. */ size_t nnode() const; /** \return Number of dimensions. */ size_t ndim() const; /** \return Number of DOFs. */ size_t ndof() const; /** \return DOFs per node [#nnode, #ndim] */ xt::xtensor dofs() const; /** Assemble from matrices stored per element. \warning Ignores any off-diagonal terms. \param elemmat [#nelem, #nne * #ndim, #nne * #ndim]. */ virtual void assemble(const xt::xtensor& elemmat); /** Set all (diagonal) matrix components. \param A The matrix [#ndof]. */ virtual void set(const xt::xtensor& A); /** Return matrix as diagonal matrix. \param [#ndof]. */ virtual xt::xtensor Todiagonal() const; /** Dot-product \f$ b_i = A_{ij} x_j \f$. \param x nodevec [#nelem, #ndim]. \return b nodevec overwritten [#nelem, #ndim]. */ xt::xtensor Dot(const xt::xtensor& x) const; /** Same as Dot(const xt::xtensor&, xt::xtensor& b) but writing to preallocated data. \param x nodevec [#nelem, #ndim]. \param b nodevec overwritten [#nelem, #ndim]. */ virtual void dot(const xt::xtensor& x, xt::xtensor& b) const; /** Dot-product \f$ b_i = A_{ij} x_j \f$. \param x dofval [#ndof]. \return b dofval overwritten [#ndof]. */ xt::xtensor Dot(const xt::xtensor& x) const; /** Same as Dot(const xt::xtensor&, xt::xtensor& b) but writing to preallocated data. \param x dofval [#ndof]. \param b dofval overwritten [#ndof]. */ virtual void dot(const xt::xtensor& x, xt::xtensor& b) const; /** Solve \f$ x = A^{-1} b \f$. \param b nodevec [nelem, ndim]. \return x nodevec [nelem, ndim]. */ xt::xtensor Solve(const xt::xtensor& b); /** Same as Solve(const xt::xtensor&) but writing to preallocated data. \param b nodevec [nelem, ndim]. \param x nodevec overwritten [nelem, ndim]. */ virtual void solve(const xt::xtensor& b, xt::xtensor& x); /** Same as Solve(const xt::xtensor&) but for "dofval" input and output. \param b dofval [ndof]. \return x dofval [ndof]. */ xt::xtensor Solve(const xt::xtensor& b); /** Same as Solve(const xt::xtensor&) but writing to preallocated data. \param b dofval [ndof]. \param x dofval overwritten [ndof]. */ virtual void solve(const xt::xtensor& b, xt::xtensor& x); protected: xt::xtensor m_conn; ///< Connectivity [#nelem, #nne]. xt::xtensor m_dofs; ///< DOF-numbers per node [#nnode, #ndim]. size_t m_nelem; ///< See nelem(). size_t m_nne; ///< See nne(). size_t m_nnode; ///< See nnode(). size_t m_ndim; ///< See ndim(). size_t m_ndof; ///< See ndof(). bool m_changed = true; ///< Signal changes to data. private: xt::xtensor m_A; ///< The matrix. xt::xtensor m_inv; /// Inverse of the matrix. void factorize(); ///< Compute inverse (automatically evaluated by "solve"). }; } // namespace GooseFEM #include "MatrixDiagonal.hpp" #endif diff --git a/include/GooseFEM/MatrixDiagonal.hpp b/include/GooseFEM/MatrixDiagonal.hpp index b15e3a8..27c21fa 100644 --- a/include/GooseFEM/MatrixDiagonal.hpp +++ b/include/GooseFEM/MatrixDiagonal.hpp @@ -1,181 +1,182 @@ /** Implementation of MatrixDiagonal.h \file MatrixDiagonal.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MATRIXDIAGONAL_HPP #define GOOSEFEM_MATRIXDIAGONAL_HPP #include "MatrixDiagonal.h" namespace GooseFEM { +template inline MatrixDiagonal::MatrixDiagonal( - const xt::xtensor& conn, - const xt::xtensor& dofs) + const C& conn, + const D& dofs) : m_conn(conn), m_dofs(dofs) { m_nelem = m_conn.shape(0); m_nne = m_conn.shape(1); m_nnode = m_dofs.shape(0); m_ndim = m_dofs.shape(1); m_ndof = xt::amax(m_dofs)() + 1; m_A = xt::empty({m_ndof}); m_inv = xt::empty({m_ndof}); GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); } inline size_t MatrixDiagonal::nelem() const { return m_nelem; } inline size_t MatrixDiagonal::nne() const { return m_nne; } inline size_t MatrixDiagonal::nnode() const { return m_nnode; } inline size_t MatrixDiagonal::ndim() const { return m_ndim; } inline size_t MatrixDiagonal::ndof() const { return m_ndof; } inline xt::xtensor MatrixDiagonal::dofs() const { return m_dofs; } inline void MatrixDiagonal::factorize() { if (!m_changed) { return; } #pragma omp parallel for for (size_t d = 0; d < m_ndof; ++d) { m_inv(d) = 1.0 / m_A(d); } m_changed = false; } inline void MatrixDiagonal::set(const xt::xtensor& A) { GOOSEFEM_ASSERT(A.size() == m_ndof); std::copy(A.begin(), A.end(), m_A.begin()); m_changed = true; } inline void MatrixDiagonal::assemble(const xt::xtensor& elemmat) { GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); GOOSEFEM_ASSERT(Element::isDiagonal(elemmat)); m_A.fill(0.0); for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { m_A(m_dofs(m_conn(e, m), i)) += elemmat(e, m * m_ndim + i, m * m_ndim + i); } } } m_changed = true; } inline void MatrixDiagonal::dot(const xt::xtensor& x, xt::xtensor& b) const { GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { b(m, i) = m_A(m_dofs(m, i)) * x(m, i); } } } inline void MatrixDiagonal::dot(const xt::xtensor& x, xt::xtensor& b) const { GOOSEFEM_ASSERT(x.size() == m_ndof); GOOSEFEM_ASSERT(b.size() == m_ndof); xt::noalias(b) = m_A * x; } inline void MatrixDiagonal::solve(const xt::xtensor& b, xt::xtensor& x) { GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim})); this->factorize(); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { x(m, i) = m_inv(m_dofs(m, i)) * b(m, i); } } } inline void MatrixDiagonal::solve(const xt::xtensor& b, xt::xtensor& x) { GOOSEFEM_ASSERT(b.size() == m_ndof); GOOSEFEM_ASSERT(x.size() == m_ndof); this->factorize(); xt::noalias(x) = m_inv * b; } inline xt::xtensor MatrixDiagonal::Todiagonal() const { return m_A; } inline xt::xtensor MatrixDiagonal::Dot(const xt::xtensor& x) const { xt::xtensor b = xt::empty({m_nnode, m_ndim}); this->dot(x, b); return b; } inline xt::xtensor MatrixDiagonal::Dot(const xt::xtensor& x) const { xt::xtensor b = xt::empty({m_ndof}); this->dot(x, b); return b; } inline xt::xtensor MatrixDiagonal::Solve(const xt::xtensor& b) { xt::xtensor x = xt::empty({m_nnode, m_ndim}); this->solve(b, x); return x; } inline xt::xtensor MatrixDiagonal::Solve(const xt::xtensor& b) { xt::xtensor x = xt::empty({m_ndof}); this->solve(b, x); return x; } } // namespace GooseFEM #endif diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h index f09744b..3a49073 100644 --- a/include/GooseFEM/Mesh.h +++ b/include/GooseFEM/Mesh.h @@ -1,1497 +1,1524 @@ /** 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" +#include "Vector.h" +#include "MatrixDiagonal.h" namespace GooseFEM { /** Generic mesh operations, and simple mesh definitions. */ namespace Mesh { /** Enumerator for element-types */ enum class ElementType { Unknown, ///< Unknown element-type Quad4, ///< Quadrilateral: 4-noded element in 2-d Hex8, ///< Hexahedron: 8-noded element in 3-d Tri3 ///< Triangle: 3-noded element in 2-d }; /** Extract the element type based on the connectivity. \param coor Nodal coordinates [nnode, ndim]. \param conn Connectivity [nelem, nne]. \return ElementType(). */ template inline ElementType defaultElementType(const S& coor, const T& conn); /** CRTP base class for regular meshes. */ template class RegularBase { public: /** Underlying type. */ using derived_type = D; /** Number of elements. \return unsigned int */ auto nelem() const; /** Number of nodes. \return unsigned int */ auto nnode() const; /** Number of nodes-per-element == 4. \return unsigned int */ auto nne() const; /** Number of dimensions == 2. \return unsigned int */ auto ndim() const; /** Number of elements in x-direction == width of the mesh in units of #h. \return unsigned int */ auto nelx() const; /** Number of elements in y-direction == height of the mesh, in units of #h, \return unsigned int */ auto nely() const; /** Linear edge size of one 'block'. \return double */ auto h() const; /** The ElementType(). \return element type */ auto getElementType() const; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ auto coor() const; /** Connectivity [#nelem, #nne]. \return nodes per element */ auto conn() const; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ auto dofs() const; /** DOF-numbers for the case that the periodicity if fully eliminated. \return DOF numbers for each node [#nnode, #ndim]. */ auto dofsPeriodic() const; /** Periodic node pairs, in two columns: (independent, dependent). \return [ntyings, #ndim]. */ auto nodesPeriodic() const; /** Reference node to use for periodicity, because all corners are tied to it. \return Node number. */ auto nodesOrigin() const; private: auto derived_cast() -> derived_type&; auto derived_cast() const -> const derived_type&; }; /** CRTP base class for regular meshes in 2d. */ template class RegularBase2d : public RegularBase { public: /** Underlying type. */ using derived_type = D; /** Nodes along the bottom edge (y = 0), in order of increasing x. \return List of node numbers. */ auto nodesBottomEdge() const; /** Nodes along the top edge (y = #nely * #h), in order of increasing x. \return List of node numbers. */ auto nodesTopEdge() const; /** Nodes along the left edge (x = 0), in order of increasing y. \return List of node numbers. */ auto nodesLeftEdge() const; /** Nodes along the right edge (x = #nelx * #h), in order of increasing y. \return List of node numbers. */ auto 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. */ auto 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. */ auto 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. */ auto 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. */ auto nodesRightOpenEdge() const; /** The bottom-left corner node (at x = 0, y = 0). Same as nodesBottomEdge()[0] and nodesLeftEdge()[0]. \return Node number. */ auto nodesBottomLeftCorner() const; /** The bottom-right corner node (at x = #nelx * #h, y = 0). Same as nodesBottomEdge()[-1] and nodesRightEdge()[0]. \return Node number. */ auto nodesBottomRightCorner() const; /** The top-left corner node (at x = 0, y = #nely * #h). Same as nodesTopEdge()[0] and nodesRightEdge()[-1]. \return Node number. */ auto nodesTopLeftCorner() const; /** The top-right corner node (at x = #nelx * #h, y = #nely * #h). Same as nodesTopEdge()[-1] and nodesRightEdge()[-1]. \return Node number. */ auto nodesTopRightCorner() const; /** Alias of nodesBottomLeftCorner(). \return Node number. */ auto nodesLeftBottomCorner() const; /** Alias of nodesTopLeftCorner(). \return Node number. */ auto nodesLeftTopCorner() const; /** Alias of nodesBottomRightCorner(). \return Node number. */ auto nodesRightBottomCorner() const; /** Alias of nodesTopRightCorner(). \return Node number. */ auto nodesRightTopCorner() const; private: auto derived_cast() -> derived_type&; auto derived_cast() const -> const derived_type&; friend class RegularBase; xt::xtensor nodesPeriodic_impl() const; auto nodesOrigin_impl() const; }; /** CRTP base class for regular meshes in 3d. */ template class RegularBase3d : public RegularBase { public: /** Underlying type. */ using derived_type = D; /** Number of elements in y-direction == height of the mesh, in units of #h, \return unsigned int */ auto nelz() const; /** Nodes along the bottom face (y = 0). \return List of node numbers. */ auto nodesBottom() const; /** Nodes along the top face (y = #nely * #h). \return List of node numbers. */ auto nodesTop() const; /** Nodes along the left face (x = 0). \return List of node numbers. */ auto nodesLeft() const; /** Nodes along the right face (x = #nelx * #h). \return List of node numbers. */ auto nodesRight() const; /** Nodes along the front face (z = 0). \return List of node numbers. */ auto nodesFront() const; /** Nodes along the back face (z = #nelz * #h). \return List of node numbers. */ auto nodesBack() const; /** Nodes along the edge at the intersection of the front and bottom faces (z = 0 and y = 0). \return List of node numbers. */ auto nodesFrontBottomEdge() const; /** Nodes along the edge at the intersection of the front and top faces (z = 0 and y = #nely * #h). \return List of node numbers. */ auto nodesFrontTopEdge() const; /** Nodes along the edge at the intersection of the front and left faces (z = 0 and x = 0). \return List of node numbers. */ auto nodesFrontLeftEdge() const; /** Nodes along the edge at the intersection of the front and right faces (z = 0 and x = #nelx * #h). \return List of node numbers. */ auto nodesFrontRightEdge() const; /** Nodes along the edge at the intersection of the back and bottom faces (z = #nelz * #h and y = #nely * #h). \return List of node numbers. */ auto nodesBackBottomEdge() const; /** Nodes along the edge at the intersection of the back and top faces (z = #nelz * #h and x = 0). \return List of node numbers. */ auto nodesBackTopEdge() const; /** Nodes along the edge at the intersection of the back and left faces (z = #nelz * #h and x = #nelx * #h). \return List of node numbers. */ auto nodesBackLeftEdge() const; /** Nodes along the edge at the intersection of the back and right faces (? = #nelz * #h and ? = ?). \return List of node numbers. */ auto nodesBackRightEdge() const; /** Nodes along the edge at the intersection of the bottom and left faces (y = 0 and x = 0). \return List of node numbers. */ auto nodesBottomLeftEdge() const; /** Nodes along the edge at the intersection of the bottom and right faces (y = 0 and x = #nelx * #h). \return List of node numbers. */ auto nodesBottomRightEdge() const; /** Nodes along the edge at the intersection of the top and left faces (y = 0 and x = #nelx * #h). \return List of node numbers. */ auto nodesTopLeftEdge() const; /** Nodes along the edge at the intersection of the top and right faces (y = #nely * #h and x = #nelx * #h). \return List of node numbers. */ auto nodesTopRightEdge() const; /** Alias of nodesFrontBottomEdge() \return List of node numbers. */ auto nodesBottomFrontEdge() const; /** Alias of nodesBackBottomEdge() \return List of node numbers. */ auto nodesBottomBackEdge() const; /** Alias of nodesFrontTopEdge() \return List of node numbers. */ auto nodesTopFrontEdge() const; /** Alias of nodesBackTopEdge() \return List of node numbers. */ auto nodesTopBackEdge() const; /** Alias of nodesBottomLeftEdge() \return List of node numbers. */ auto nodesLeftBottomEdge() const; /** Alias of nodesFrontLeftEdge() \return List of node numbers. */ auto nodesLeftFrontEdge() const; /** Alias of nodesBackLeftEdge() \return List of node numbers. */ auto nodesLeftBackEdge() const; /** Alias of nodesTopLeftEdge() \return List of node numbers. */ auto nodesLeftTopEdge() const; /** Alias of nodesBottomRightEdge() \return List of node numbers. */ auto nodesRightBottomEdge() const; /** Alias of nodesTopRightEdge() \return List of node numbers. */ auto nodesRightTopEdge() const; /** Alias of nodesFrontRightEdge() \return List of node numbers. */ auto nodesRightFrontEdge() const; /** Alias of nodesBackRightEdge() \return List of node numbers. */ auto nodesRightBackEdge() const; /** Nodes along the front face excluding edges. Same as different between nodesFront() and [nodesFrontBottomEdge(), nodesFrontTopEdge(), nodesFrontLeftEdge(), nodesFrontRightEdge()] \return list of node numbers. */ auto nodesFrontFace() const; /** Nodes along the back face excluding edges. Same as different between nodesBack() and [nodesBackBottomEdge(), nodesBackTopEdge(), nodesBackLeftEdge(), nodesBackRightEdge()] \return list of node numbers. */ auto nodesBackFace() const; /** Nodes along the left face excluding edges. Same as different between nodesLeft() and [nodesFrontLeftEdge(), nodesBackLeftEdge(), nodesBottomLeftEdge(), nodesTopLeftEdge()] \return list of node numbers. */ auto nodesLeftFace() const; /** Nodes along the right face excluding edges. Same as different between nodesRight() and [nodesFrontRightEdge(), nodesBackRightEdge(), nodesBottomRightEdge(), nodesTopRightEdge()] \return list of node numbers. */ auto nodesRightFace() const; /** Nodes along the bottom face excluding edges. Same as different between nodesBottom() and [nodesBackBottomEdge(), nodesBackTopEdge(), nodesBackLeftEdge(), nodesBackRightEdge()] \return list of node numbers. */ auto nodesBottomFace() const; /** Nodes along the top face excluding edges. Same as different between nodesTop() and [nodesFrontBottomEdge(), nodesFrontTopEdge(), nodesFrontLeftEdge(), nodesFrontRightEdge()] \return list of node numbers. */ auto nodesTopFace() const; /** Same as nodesFrontBottomEdge() but without corners. \return List of node numbers. */ auto nodesFrontBottomOpenEdge() const; /** Same as nodesFrontTopEdge() but without corners. \return List of node numbers. */ auto nodesFrontTopOpenEdge() const; /** Same as nodesFrontLeftEdge() but without corners. \return List of node numbers. */ auto nodesFrontLeftOpenEdge() const; /** Same as nodesFrontRightEdge() but without corners. \return List of node numbers. */ auto nodesFrontRightOpenEdge() const; /** Same as nodesBackBottomEdge() but without corners. \return List of node numbers. */ auto nodesBackBottomOpenEdge() const; /** Same as nodesBackTopEdge() but without corners. \return List of node numbers. */ auto nodesBackTopOpenEdge() const; /** Same as nodesBackLeftEdge() but without corners. \return List of node numbers. */ auto nodesBackLeftOpenEdge() const; /** Same as nodesBackRightEdge() but without corners. \return List of node numbers. */ auto nodesBackRightOpenEdge() const; /** Same as nodesBottomLeftEdge() but without corners. \return List of node numbers. */ auto nodesBottomLeftOpenEdge() const; /** Same as nodesBottomRightEdge() but without corners. \return List of node numbers. */ auto nodesBottomRightOpenEdge() const; /** Same as nodesTopLeftEdge() but without corners. \return List of node numbers. */ auto nodesTopLeftOpenEdge() const; /** Same as nodesTopRightEdge() but without corners. \return List of node numbers. */ auto nodesTopRightOpenEdge() const; /** Alias of nodesFrontBottomOpenEdge(). \return List of node numbers. */ auto nodesBottomFrontOpenEdge() const; /** Alias of nodesBackBottomOpenEdge(). \return List of node numbers. */ auto nodesBottomBackOpenEdge() const; /** Alias of nodesFrontTopOpenEdge(). \return List of node numbers. */ auto nodesTopFrontOpenEdge() const; /** Alias of nodesBackTopOpenEdge(). \return List of node numbers. */ auto nodesTopBackOpenEdge() const; /** Alias of nodesBottomLeftOpenEdge(). \return List of node numbers. */ auto nodesLeftBottomOpenEdge() const; /** Alias of nodesFrontLeftOpenEdge(). \return List of node numbers. */ auto nodesLeftFrontOpenEdge() const; /** Alias of nodesBackLeftOpenEdge(). \return List of node numbers. */ auto nodesLeftBackOpenEdge() const; /** Alias of nodesTopLeftOpenEdge(). \return List of node numbers. */ auto nodesLeftTopOpenEdge() const; /** Alias of nodesBottomRightOpenEdge(). \return List of node numbers. */ auto nodesRightBottomOpenEdge() const; /** Alias of nodesTopRightOpenEdge(). \return List of node numbers. */ auto nodesRightTopOpenEdge() const; /** Alias of nodesFrontRightOpenEdge(). \return List of node numbers. */ auto nodesRightFrontOpenEdge() const; /** Alias of nodesBackRightOpenEdge(). \return List of node numbers. */ auto nodesRightBackOpenEdge() const; /** Front-Bottom-Left corner node. \return Node number. */ auto nodesFrontBottomLeftCorner() const; /** Front-Bottom-Right corner node. \return Node number. */ auto nodesFrontBottomRightCorner() const; /** Front-Top-Left corner node. \return Node number. */ auto nodesFrontTopLeftCorner() const; /** Front-Top-Right corner node. \return Node number. */ auto nodesFrontTopRightCorner() const; /** Back-Bottom-Left corner node. \return Node number. */ auto nodesBackBottomLeftCorner() const; /** Back-Bottom-Right corner node. \return Node number. */ auto nodesBackBottomRightCorner() const; /** Back-Top-Left corner node. \return Node number. */ auto nodesBackTopLeftCorner() const; /** Back-Top-Right corner node. \return Node number. */ auto nodesBackTopRightCorner() const; /** Alias of nodesFrontBottomLeftCorner(). \return Node number. */ auto nodesFrontLeftBottomCorner() const; /** Alias of nodesFrontBottomLeftCorner(). \return Node number. */ auto nodesBottomFrontLeftCorner() const; /** Alias of nodesFrontBottomLeftCorner(). \return Node number. */ auto nodesBottomLeftFrontCorner() const; /** Alias of nodesFrontBottomLeftCorner(). \return Node number. */ auto nodesLeftFrontBottomCorner() const; /** Alias of nodesFrontBottomLeftCorner(). \return Node number. */ auto nodesLeftBottomFrontCorner() const; /** Alias of nodesFrontBottomRightCorner(). \return Node number. */ auto nodesFrontRightBottomCorner() const; /** Alias of nodesFrontBottomRightCorner(). \return Node number. */ auto nodesBottomFrontRightCorner() const; /** Alias of nodesFrontBottomRightCorner(). \return Node number. */ auto nodesBottomRightFrontCorner() const; /** Alias of nodesFrontBottomRightCorner(). \return Node number. */ auto nodesRightFrontBottomCorner() const; /** Alias of nodesFrontBottomRightCorner(). \return Node number. */ auto nodesRightBottomFrontCorner() const; /** Alias of nodesFrontTopLeftCorner(). \return Node number. */ auto nodesFrontLeftTopCorner() const; /** Alias of nodesFrontTopLeftCorner(). \return Node number. */ auto nodesTopFrontLeftCorner() const; /** Alias of nodesFrontTopLeftCorner(). \return Node number. */ auto nodesTopLeftFrontCorner() const; /** Alias of nodesFrontTopLeftCorner(). \return Node number. */ auto nodesLeftFrontTopCorner() const; /** Alias of nodesFrontTopLeftCorner(). \return Node number. */ auto nodesLeftTopFrontCorner() const; /** Alias of nodesFrontTopRightCorner(). \return Node number. */ auto nodesFrontRightTopCorner() const; /** Alias of nodesFrontTopRightCorner(). \return Node number. */ auto nodesTopFrontRightCorner() const; /** Alias of nodesFrontTopRightCorner(). \return Node number. */ auto nodesTopRightFrontCorner() const; /** Alias of nodesFrontTopRightCorner(). \return Node number. */ auto nodesRightFrontTopCorner() const; /** Alias of nodesFrontTopRightCorner(). \return Node number. */ auto nodesRightTopFrontCorner() const; /** Alias of nodesBackBottomLeftCorner(). \return Node number. */ auto nodesBackLeftBottomCorner() const; /** Alias of nodesBackBottomLeftCorner(). \return Node number. */ auto nodesBottomBackLeftCorner() const; /** Alias of nodesBackBottomLeftCorner(). \return Node number. */ auto nodesBottomLeftBackCorner() const; /** Alias of nodesBackBottomLeftCorner(). \return Node number. */ auto nodesLeftBackBottomCorner() const; /** Alias of nodesBackBottomLeftCorner(). \return Node number. */ auto nodesLeftBottomBackCorner() const; /** Alias of nodesBackBottomRightCorner(). \return Node number. */ auto nodesBackRightBottomCorner() const; /** Alias of nodesBackBottomRightCorner(). \return Node number. */ auto nodesBottomBackRightCorner() const; /** Alias of nodesBackBottomRightCorner(). \return Node number. */ auto nodesBottomRightBackCorner() const; /** Alias of nodesBackBottomRightCorner(). \return Node number. */ auto nodesRightBackBottomCorner() const; /** Alias of nodesBackBottomRightCorner(). \return Node number. */ auto nodesRightBottomBackCorner() const; /** Alias of nodesBackTopLeftCorner(). \return Node number. */ auto nodesBackLeftTopCorner() const; /** Alias of nodesBackTopLeftCorner(). \return Node number. */ auto nodesTopBackLeftCorner() const; /** Alias of nodesBackTopLeftCorner(). \return Node number. */ auto nodesTopLeftBackCorner() const; /** Alias of nodesBackTopLeftCorner(). \return Node number. */ auto nodesLeftBackTopCorner() const; /** Alias of nodesBackTopLeftCorner(). \return Node number. */ auto nodesLeftTopBackCorner() const; /** Alias of nodesBackTopRightCorner(). \return Node number. */ auto nodesBackRightTopCorner() const; /** Alias of nodesBackTopRightCorner(). \return Node number. */ auto nodesTopBackRightCorner() const; /** Alias of nodesBackTopRightCorner(). \return Node number. */ auto nodesTopRightBackCorner() const; /** Alias of nodesBackTopRightCorner(). \return Node number. */ auto nodesRightBackTopCorner() const; /** Alias of nodesBackTopRightCorner(). \return Node number. */ auto nodesRightTopBackCorner() const; private: auto derived_cast() -> derived_type&; auto derived_cast() const -> const derived_type&; friend class RegularBase; xt::xtensor nodesPeriodic_impl() const; auto nodesOrigin_impl() const; }; /** 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" [nnode, ndim]. \param coor_b Nodal coordinates of mesh "b" [nnode, ndim]. \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. \return Overlapping nodes. */ template inline xt::xtensor overlapping(const S& coor_a, const T& coor_b, double rtol = 1e-5, double atol = 1e-8); /** Stitch two mesh objects, specifying overlapping nodes by hand. */ class ManualStitch { public: ManualStitch() = default; /** \param coor_a Nodal coordinates of mesh "a" [nnode, ndim]. \param conn_a Connectivity of mesh "a" [nelem, nne]. \param overlapping_nodes_a Node-numbers of mesh "a" that overlap with mesh "b" [n]. \param coor_b Nodal coordinates of mesh "b" [nnode, ndim]. \param conn_b Connectivity of mesh "b" [nelem, nne]. \param overlapping_nodes_b Node-numbers of mesh "b" that overlap with mesh "a" [n]. \param check_position If ``true`` the nodes are checked for position overlap. \param rtol Relative tolerance for check on position overlap. \param atol Absolute tolerance for check on position overlap. */ template ManualStitch( const CA& coor_a, const EA& conn_a, const NA& overlapping_nodes_a, const CB& coor_b, const EB& conn_b, const NB& overlapping_nodes_b, bool check_position = true, double rtol = 1e-5, double atol = 1e-8); /** Number of sub meshes == 2. \return unsigned int */ size_t nmesh() const; /** Number of elements. \return unsigned int */ size_t nelem() const; /** Number of nodes. \return unsigned int */ size_t nnode() const; /** Number of nodes-per-element. \return unsigned int */ size_t nne() const; /** Number of dimensions. \return unsigned int */ size_t ndim() const; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ xt::xtensor coor() const; /** Connectivity [#nelem, #nne]. \return nodes per element */ xt::xtensor conn() const; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ xt::xtensor dofs() const; /** Node-map per sub-mesh. \return nodes per mesh */ std::vector> nodemap() const; /** Element-map per sub-mesh. \return elements per mesh */ std::vector> elemmap() const; /** \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return Node-map for a given mesh. */ xt::xtensor nodemap(size_t mesh_index) const; /** \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return Element-map for a given mesh. */ xt::xtensor elemmap(size_t mesh_index) const; /** Convert set of node numbers for an original mesh to the stitched mesh. \param set List of node numbers. \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return List of node numbers for the stitched mesh. */ template T nodeset(const T& set, size_t mesh_index) const; /** Convert set of element numbers for an original mesh to the stitched mesh. \param set List of element numbers. \param mesh_index Index of the mesh ("a" = 1, "b" = 1). \return List of element numbers for the stitched mesh. */ template T elemset(const T& set, size_t mesh_index) const; private: xt::xtensor m_coor; xt::xtensor m_conn; xt::xtensor m_map_b; size_t m_nnd_a; size_t m_nel_a; size_t m_nel_b; }; /** Stitch mesh objects, automatically searching for overlapping nodes. */ class Stitch { public: /** \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. */ Stitch(double rtol = 1e-5, double atol = 1e-8); /** Add mesh to be stitched. \param coor Nodal coordinates [nnode, ndim]. \param conn Connectivity [nelem, nne]. */ template void push_back(const C& coor, const E& conn); /** Number of sub meshes. \return unsigned int */ size_t nmesh() const; /** Number of elements. \return unsigned int */ size_t nelem() const; /** Number of nodes. \return unsigned int */ size_t nnode() const; /** Number of nodes-per-element. \return unsigned int */ size_t nne() const; /** Number of dimensions. \return unsigned int */ size_t ndim() const; /** Nodal coordinates [#nnode, #ndim]. \return coordinates per node */ xt::xtensor coor() const; /** Connectivity [#nelem, #nne]. \return nodes per element */ xt::xtensor conn() const; /** DOF numbers for each node (numbered sequentially) [#nnode, #ndim]. \return DOFs per node */ xt::xtensor dofs() const; /** Node-map per sub-mesh. \return nodes per mesh */ std::vector> nodemap() const; /** Element-map per sub-mesh. \return elements per mesh */ std::vector> elemmap() const; /** The node numbers in the stitched mesh that are coming from a specific sub-mesh. \param mesh_index Index of the sub-mesh. \return List of node numbers. */ xt::xtensor nodemap(size_t mesh_index) const; /** The element numbers in the stitched mesh that are coming from a specific sub-mesh. \param mesh_index Index of the sub-mesh. \return List of element numbers. */ xt::xtensor elemmap(size_t mesh_index) const; /** Convert set of node-numbers for a sub-mesh to the stitched mesh. \param set List of node numbers. \param mesh_index Index of the sub-mesh. \return List of node numbers for the stitched mesh. */ template T nodeset(const T& set, size_t mesh_index) const; /** Convert set of element-numbers for a sub-mesh to the stitched mesh. \param set List of element numbers. \param mesh_index Index of the sub-mesh. \return List of element numbers for the stitched mesh. */ template T elemset(const T& set, size_t mesh_index) const; /** Combine set of node numbers for an original to the final mesh (removes duplicates). \param set List of node numbers per mesh. \return List of node numbers for the stitched mesh. */ template T nodeset(const std::vector& set) const; /** \copydoc nodeset(const std::vector&) const */ template T nodeset(std::initializer_list set) const; /** Combine set of element numbers for an original to the final mesh. \param set List of element numbers per mesh. \return List of element numbers for the stitched mesh. */ template T elemset(const std::vector& set) const; /** \copydoc elemset(const std::vector&) const */ template T elemset(std::initializer_list set) const; protected: xt::xtensor m_coor; ///< Nodal coordinates [#nnode, #ndim] xt::xtensor m_conn; ///< Connectivity [#nelem, #nne] std::vector> m_map; ///< See nodemap(size_t) std::vector m_nel; ///< Number of elements per sub-mesh. std::vector m_el_offset; ///< First element of every sub-mesh. double m_rtol; ///< Relative tolerance to find overlapping nodes. double m_atol; ///< Absolute tolerance to find overlapping nodes. }; /** Vertically stack meshes. */ class Vstack : public Stitch { public: /** \param check_overlap Check if nodes are overlapping when adding a mesh. \param rtol Relative tolerance for position match. \param atol Absolute tolerance for position match. */ Vstack(bool check_overlap = true, double rtol = 1e-5, double atol = 1e-8); /** Add a mesh to the top of the current stack. Each time the current `nodes_bottom` are stitched with the then highest `nodes_top`. \param coor Nodal coordinates [nnode, ndim]. \param conn Connectivity [nelem, nne]. \param nodes_bottom Nodes along the bottom edge [n]. \param nodes_top Nodes along the top edge [n]. */ template void push_back(const C& coor, const E& conn, const N& nodes_bottom, const N& nodes_top); private: std::vector> m_nodes_bot; ///< Bottom nodes of each mesh (renumbered). std::vector> m_nodes_top; ///< Top nodes of each mesh (renumbered). bool m_check_overlap; ///< Check if nodes are overlapping when adding a mesh. }; /** Renumber indices to lowest possible index. For example: \f$ \begin{bmatrix} 0 & 1 \\ 5 & 4 \end{bmatrix} \f$ is renumbered to \f$ \begin{bmatrix} 0 & 1 \\ 3 & 2 \end{bmatrix} \f$ Or, in pseudo-code, the result of this function is that: dofs = renumber(dofs) sort(unique(dofs[:])) == range(max(dofs+1)) \note One can use the wrapper function renumber(). This class gives more advanced features. */ class Renumber { public: Renumber() = default; /** \param dofs DOF-numbers. */ template Renumber(const T& dofs); /** 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 [nnode, ndim]. \return Renumbered DOF-numbers. */ template inline T renumber(const T& dofs); /** Reorder to lowest possible index, in specific order. For example for ``Reorder({iiu, iip})`` after reordering: iiu = xt::range(nnu); iip = xt::range(nnp) + nnu; */ class Reorder { public: Reorder() = default; /** \param args List of (DOF-)numbers. */ template Reorder(const std::initializer_list args); /** Apply reordering to other set. \param list List of (DOF-)numbers. \return Reordered list of (DOF-)numbers. */ template T apply(const T& list) const; /** Get the list needed to reorder, e.g.: dofs_reordered(i, j) = index(dofs(i, j)) \return Reorder-index. */ xt::xtensor index() const; private: xt::xtensor m_renum; }; /** List with DOF-numbers in sequential order. The output is a sequential list of DOF-numbers for each vector-component of each node. For example for 3 nodes in 2 dimensions the output is \f$ \begin{bmatrix} 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \f$ \param nnode Number of nodes. \param ndim Number of dimensions. \return DOF-numbers. */ inline xt::xtensor dofs(size_t nnode, size_t ndim); /** Number of elements connected to each node. \param conn Connectivity [nelem, nne]. \return Coordination per node. */ template inline xt::xtensor coordination(const E& conn); /** Elements connected to each node. \param conn Connectivity. \param sorted If ``true`` the output is sorted. \return Elements per node. */ template inline std::vector> elem2node(const E& 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. */ template inline xt::xtensor edgesize(const C& coor, const E& 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. */ template inline xt::xtensor edgesize(const C& coor, const E& conn); /** Coordinates of the center of each element. \param coor Nodal coordinates. \param conn Connectivity. \param type ElementType. \return Center of each element. */ template inline xt::xtensor centers(const C& coor, const E& 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. */ template inline xt::xtensor centers(const C& coor, const E& 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]`` */ template inline xt::xtensor elemmap2nodemap(const T& elem_map, const C& coor, const E& 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]`` */ template inline xt::xtensor elemmap2nodemap(const T& elem_map, const C& coor, const E& conn); +/** +Compute the center of gravity of a mesh. + +\tparam C e.g. `xt::xtensor` +\tparam E e.g. `xt::xtensor` +\param coor Nodal coordinates `[nnode, ndim]`. +\param conn Connectivity `[nelem, nne]`. +\param type ElementType. +\return Center of gravity `[ndim]`. +*/ +template +inline xt::xtensor center_of_gravity(const C& coor, const E& conn, ElementType type); + +/** +Compute the center of gravity of a mesh. + +\tparam C e.g. `xt::xtensor` +\tparam E e.g. `xt::xtensor` +\param coor Nodal coordinates `[nnode, ndim]`. +\param conn Connectivity `[nelem, nne]`. +\return Center of gravity `[ndim]`. +*/ +template +inline xt::xtensor center_of_gravity(const C& coor, const E& conn); + } // namespace Mesh } // namespace GooseFEM #include "Mesh.hpp" #endif diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp index f65a64c..0feb034 100644 --- a/include/GooseFEM/Mesh.hpp +++ b/include/GooseFEM/Mesh.hpp @@ -1,1761 +1,1789 @@ /** 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" #include "assertions.h" namespace GooseFEM { namespace Mesh { template inline ElementType defaultElementType(const S& coor, const T& conn) { GOOSEFEM_ASSERT(coor.dimension() == 2); GOOSEFEM_ASSERT(conn.dimension() == 2); if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) { return ElementType::Tri3; } if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) { return ElementType::Quad4; } if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) { return ElementType::Hex8; } throw std::runtime_error("Element-type not implemented"); } template inline auto RegularBase::nelem() const { return derived_cast().m_nelem; } template inline auto RegularBase::nnode() const { return derived_cast().m_nnode; } template inline auto RegularBase::nne() const { return derived_cast().m_nne; } template inline auto RegularBase::ndim() const { return derived_cast().m_ndim; } template inline auto RegularBase::nelx() const { return derived_cast().nelx_impl(); } template inline auto RegularBase::nely() const { return derived_cast().nely_impl(); } template inline auto RegularBase::h() const { return derived_cast().m_h; } template inline auto RegularBase::getElementType() const { return derived_cast().getElementType_impl(); } template inline auto RegularBase::coor() const { return derived_cast().coor_impl(); } template inline auto RegularBase::conn() const { return derived_cast().conn_impl(); } template inline auto RegularBase::dofs() const { return GooseFEM::Mesh::dofs(this->nnode(), this->ndim()); } template inline auto RegularBase::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); } template inline auto RegularBase::derived_cast() -> derived_type& { return *static_cast(this); } template inline auto RegularBase::derived_cast() const -> const derived_type& { return *static_cast(this); } template inline auto RegularBase::nodesPeriodic() const { return derived_cast().nodesPeriodic_impl(); } template inline auto RegularBase::nodesOrigin() const { return derived_cast().nodesOrigin_impl(); } template inline auto RegularBase2d::nodesBottomEdge() const { return derived_cast().nodesBottomEdge_impl(); } template inline auto RegularBase2d::nodesTopEdge() const { return derived_cast().nodesTopEdge_impl(); } template inline auto RegularBase2d::nodesLeftEdge() const { return derived_cast().nodesLeftEdge_impl(); } template inline auto RegularBase2d::nodesRightEdge() const { return derived_cast().nodesRightEdge_impl(); } template inline auto RegularBase2d::nodesBottomOpenEdge() const { return derived_cast().nodesBottomOpenEdge_impl(); } template inline auto RegularBase2d::nodesTopOpenEdge() const { return derived_cast().nodesTopOpenEdge_impl(); } template inline auto RegularBase2d::nodesLeftOpenEdge() const { return derived_cast().nodesLeftOpenEdge_impl(); } template inline auto RegularBase2d::nodesRightOpenEdge() const { return derived_cast().nodesRightOpenEdge_impl(); } template inline auto RegularBase2d::nodesBottomLeftCorner() const { return derived_cast().nodesBottomLeftCorner_impl(); } template inline auto RegularBase2d::nodesTopLeftCorner() const { return derived_cast().nodesTopLeftCorner_impl(); } template inline auto RegularBase2d::nodesBottomRightCorner() const { return derived_cast().nodesBottomRightCorner_impl(); } template inline auto RegularBase2d::nodesTopRightCorner() const { return derived_cast().nodesTopRightCorner_impl(); } template inline auto RegularBase2d::nodesLeftBottomCorner() const { return derived_cast().nodesBottomLeftCorner_impl(); } template inline auto RegularBase2d::nodesLeftTopCorner() const { return derived_cast().nodesTopLeftCorner_impl(); } template inline auto RegularBase2d::nodesRightBottomCorner() const { return derived_cast().nodesBottomRightCorner_impl(); } template inline auto RegularBase2d::nodesRightTopCorner() const { return derived_cast().nodesTopRightCorner_impl(); } template inline xt::xtensor RegularBase2d::nodesPeriodic_impl() const { xt::xtensor bot = derived_cast().nodesBottomOpenEdge_impl(); xt::xtensor top = derived_cast().nodesTopOpenEdge_impl(); xt::xtensor lft = derived_cast().nodesLeftOpenEdge_impl(); xt::xtensor rgt = derived_cast().nodesRightOpenEdge_impl(); std::array shape = {bot.size() + lft.size() + size_t(3), size_t(2)}; auto ret = xt::xtensor::from_shape(shape); ret(0, 0) = derived_cast().nodesBottomLeftCorner_impl(); ret(0, 1) = derived_cast().nodesBottomRightCorner_impl(); ret(1, 0) = derived_cast().nodesBottomLeftCorner_impl(); ret(1, 1) = derived_cast().nodesTopRightCorner_impl(); ret(2, 0) = derived_cast().nodesBottomLeftCorner_impl(); ret(2, 1) = derived_cast().nodesTopLeftCorner_impl(); 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; } template inline auto RegularBase2d::nodesOrigin_impl() const { return derived_cast().nodesBottomLeftCorner_impl(); } template inline auto RegularBase2d::derived_cast() -> derived_type& { return *static_cast(this); } template inline auto RegularBase2d::derived_cast() const -> const derived_type& { return *static_cast(this); } template inline auto RegularBase3d::nelz() const { return derived_cast().nelz_impl(); } template inline auto RegularBase3d::nodesBottom() const { return derived_cast().nodesBottom_impl(); } template inline auto RegularBase3d::nodesTop() const { return derived_cast().nodesTop_impl(); } template inline auto RegularBase3d::nodesLeft() const { return derived_cast().nodesLeft_impl(); } template inline auto RegularBase3d::nodesRight() const { return derived_cast().nodesRight_impl(); } template inline auto RegularBase3d::nodesFront() const { return derived_cast().nodesFront_impl(); } template inline auto RegularBase3d::nodesBack() const { return derived_cast().nodesBack_impl(); } template inline auto RegularBase3d::nodesFrontBottomEdge() const { return derived_cast().nodesFrontBottomEdge_impl(); } template inline auto RegularBase3d::nodesFrontTopEdge() const { return derived_cast().nodesFrontTopEdge_impl(); } template inline auto RegularBase3d::nodesFrontLeftEdge() const { return derived_cast().nodesFrontLeftEdge_impl(); } template inline auto RegularBase3d::nodesFrontRightEdge() const { return derived_cast().nodesFrontRightEdge_impl(); } template inline auto RegularBase3d::nodesBackBottomEdge() const { return derived_cast().nodesBackBottomEdge_impl(); } template inline auto RegularBase3d::nodesBackTopEdge() const { return derived_cast().nodesBackTopEdge_impl(); } template inline auto RegularBase3d::nodesBackLeftEdge() const { return derived_cast().nodesBackLeftEdge_impl(); } template inline auto RegularBase3d::nodesBackRightEdge() const { return derived_cast().nodesBackRightEdge_impl(); } template inline auto RegularBase3d::nodesBottomLeftEdge() const { return derived_cast().nodesBottomLeftEdge_impl(); } template inline auto RegularBase3d::nodesBottomRightEdge() const { return derived_cast().nodesBottomRightEdge_impl(); } template inline auto RegularBase3d::nodesTopLeftEdge() const { return derived_cast().nodesTopLeftEdge_impl(); } template inline auto RegularBase3d::nodesTopRightEdge() const { return derived_cast().nodesTopRightEdge_impl(); } template inline auto RegularBase3d::nodesBottomFrontEdge() const { return derived_cast().nodesFrontBottomEdge_impl(); } template inline auto RegularBase3d::nodesBottomBackEdge() const { return derived_cast().nodesBackBottomEdge_impl(); } template inline auto RegularBase3d::nodesTopFrontEdge() const { return derived_cast().nodesFrontTopEdge_impl(); } template inline auto RegularBase3d::nodesTopBackEdge() const { return derived_cast().nodesBackTopEdge_impl(); } template inline auto RegularBase3d::nodesLeftBottomEdge() const { return derived_cast().nodesBottomLeftEdge_impl(); } template inline auto RegularBase3d::nodesLeftFrontEdge() const { return derived_cast().nodesFrontLeftEdge_impl(); } template inline auto RegularBase3d::nodesLeftBackEdge() const { return derived_cast().nodesBackLeftEdge_impl(); } template inline auto RegularBase3d::nodesLeftTopEdge() const { return derived_cast().nodesTopLeftEdge_impl(); } template inline auto RegularBase3d::nodesRightBottomEdge() const { return derived_cast().nodesBottomRightEdge_impl(); } template inline auto RegularBase3d::nodesRightTopEdge() const { return derived_cast().nodesTopRightEdge_impl(); } template inline auto RegularBase3d::nodesRightFrontEdge() const { return derived_cast().nodesFrontRightEdge_impl(); } template inline auto RegularBase3d::nodesRightBackEdge() const { return derived_cast().nodesBackRightEdge_impl(); } template inline auto RegularBase3d::nodesFrontFace() const { return derived_cast().nodesFrontFace_impl(); } template inline auto RegularBase3d::nodesBackFace() const { return derived_cast().nodesBackFace_impl(); } template inline auto RegularBase3d::nodesLeftFace() const { return derived_cast().nodesLeftFace_impl(); } template inline auto RegularBase3d::nodesRightFace() const { return derived_cast().nodesRightFace_impl(); } template inline auto RegularBase3d::nodesBottomFace() const { return derived_cast().nodesBottomFace_impl(); } template inline auto RegularBase3d::nodesTopFace() const { return derived_cast().nodesTopFace_impl(); } template inline auto RegularBase3d::nodesFrontBottomOpenEdge() const { return derived_cast().nodesFrontBottomOpenEdge_impl(); } template inline auto RegularBase3d::nodesFrontTopOpenEdge() const { return derived_cast().nodesFrontTopOpenEdge_impl(); } template inline auto RegularBase3d::nodesFrontLeftOpenEdge() const { return derived_cast().nodesFrontLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesFrontRightOpenEdge() const { return derived_cast().nodesFrontRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesBackBottomOpenEdge() const { return derived_cast().nodesBackBottomOpenEdge_impl(); } template inline auto RegularBase3d::nodesBackTopOpenEdge() const { return derived_cast().nodesBackTopOpenEdge_impl(); } template inline auto RegularBase3d::nodesBackLeftOpenEdge() const { return derived_cast().nodesBackLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesBackRightOpenEdge() const { return derived_cast().nodesBackRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesBottomLeftOpenEdge() const { return derived_cast().nodesBottomLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesBottomRightOpenEdge() const { return derived_cast().nodesBottomRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesTopLeftOpenEdge() const { return derived_cast().nodesTopLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesTopRightOpenEdge() const { return derived_cast().nodesTopRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesBottomFrontOpenEdge() const { return derived_cast().nodesFrontBottomOpenEdge_impl(); } template inline auto RegularBase3d::nodesBottomBackOpenEdge() const { return derived_cast().nodesBackBottomOpenEdge_impl(); } template inline auto RegularBase3d::nodesTopFrontOpenEdge() const { return derived_cast().nodesFrontTopOpenEdge_impl(); } template inline auto RegularBase3d::nodesTopBackOpenEdge() const { return derived_cast().nodesBackTopOpenEdge_impl(); } template inline auto RegularBase3d::nodesLeftBottomOpenEdge() const { return derived_cast().nodesBottomLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesLeftFrontOpenEdge() const { return derived_cast().nodesFrontLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesLeftBackOpenEdge() const { return derived_cast().nodesBackLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesLeftTopOpenEdge() const { return derived_cast().nodesTopLeftOpenEdge_impl(); } template inline auto RegularBase3d::nodesRightBottomOpenEdge() const { return derived_cast().nodesBottomRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesRightTopOpenEdge() const { return derived_cast().nodesTopRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesRightFrontOpenEdge() const { return derived_cast().nodesFrontRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesRightBackOpenEdge() const { return derived_cast().nodesBackRightOpenEdge_impl(); } template inline auto RegularBase3d::nodesFrontBottomLeftCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesFrontBottomRightCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesFrontTopLeftCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesFrontTopRightCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } template inline auto RegularBase3d::nodesBackBottomLeftCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesBackBottomRightCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesBackTopLeftCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesBackTopRightCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } template inline auto RegularBase3d::nodesFrontLeftBottomCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesBottomFrontLeftCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesBottomLeftFrontCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftFrontBottomCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftBottomFrontCorner() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesFrontRightBottomCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesBottomFrontRightCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesBottomRightFrontCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesRightFrontBottomCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesRightBottomFrontCorner() const { return derived_cast().nodesFrontBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesFrontLeftTopCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesTopFrontLeftCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesTopLeftFrontCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftFrontTopCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftTopFrontCorner() const { return derived_cast().nodesFrontTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesFrontRightTopCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } template inline auto RegularBase3d::nodesTopFrontRightCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } template inline auto RegularBase3d::nodesTopRightFrontCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } template inline auto RegularBase3d::nodesRightFrontTopCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } template inline auto RegularBase3d::nodesRightTopFrontCorner() const { return derived_cast().nodesFrontTopRightCorner_impl(); } template inline auto RegularBase3d::nodesBackLeftBottomCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesBottomBackLeftCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesBottomLeftBackCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftBackBottomCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftBottomBackCorner() const { return derived_cast().nodesBackBottomLeftCorner_impl(); } template inline auto RegularBase3d::nodesBackRightBottomCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesBottomBackRightCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesBottomRightBackCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesRightBackBottomCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesRightBottomBackCorner() const { return derived_cast().nodesBackBottomRightCorner_impl(); } template inline auto RegularBase3d::nodesBackLeftTopCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesTopBackLeftCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesTopLeftBackCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftBackTopCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesLeftTopBackCorner() const { return derived_cast().nodesBackTopLeftCorner_impl(); } template inline auto RegularBase3d::nodesBackRightTopCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } template inline auto RegularBase3d::nodesTopBackRightCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } template inline auto RegularBase3d::nodesTopRightBackCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } template inline auto RegularBase3d::nodesRightBackTopCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } template inline auto RegularBase3d::nodesRightTopBackCorner() const { return derived_cast().nodesBackTopRightCorner_impl(); } template inline xt::xtensor RegularBase3d::nodesPeriodic_impl() const { xt::xtensor fro = derived_cast().nodesFrontFace_impl(); xt::xtensor bck = derived_cast().nodesBackFace_impl(); xt::xtensor lft = derived_cast().nodesLeftFace_impl(); xt::xtensor rgt = derived_cast().nodesRightFace_impl(); xt::xtensor bot = derived_cast().nodesBottomFace_impl(); xt::xtensor top = derived_cast().nodesTopFace_impl(); xt::xtensor froBot = derived_cast().nodesFrontBottomOpenEdge_impl(); xt::xtensor froTop = derived_cast().nodesFrontTopOpenEdge_impl(); xt::xtensor froLft = derived_cast().nodesFrontLeftOpenEdge_impl(); xt::xtensor froRgt = derived_cast().nodesFrontRightOpenEdge_impl(); xt::xtensor bckBot = derived_cast().nodesBackBottomOpenEdge_impl(); xt::xtensor bckTop = derived_cast().nodesBackTopOpenEdge_impl(); xt::xtensor bckLft = derived_cast().nodesBackLeftOpenEdge_impl(); xt::xtensor bckRgt = derived_cast().nodesBackRightOpenEdge_impl(); xt::xtensor botLft = derived_cast().nodesBottomLeftOpenEdge_impl(); xt::xtensor botRgt = derived_cast().nodesBottomRightOpenEdge_impl(); xt::xtensor topLft = derived_cast().nodesTopLeftOpenEdge_impl(); xt::xtensor topRgt = derived_cast().nodesTopRightOpenEdge_impl(); size_t tface = fro.size() + lft.size() + bot.size(); size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size(); size_t tnode = 7; xt::xtensor ret = xt::empty({tface + tedge + tnode, std::size_t(2)}); size_t i = 0; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesFrontBottomRightCorner_impl(); ++i; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesBackBottomRightCorner_impl(); ++i; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesBackBottomLeftCorner_impl(); ++i; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesFrontTopLeftCorner_impl(); ++i; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesFrontTopRightCorner_impl(); ++i; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesBackTopRightCorner_impl(); ++i; ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl(); ret(i, 1) = derived_cast().nodesBackTopLeftCorner_impl(); ++i; for (size_t j = 0; j < froBot.size(); ++j) { ret(i, 0) = froBot(j); ret(i, 1) = bckBot(j); ++i; } for (size_t j = 0; j < froBot.size(); ++j) { ret(i, 0) = froBot(j); ret(i, 1) = bckTop(j); ++i; } for (size_t j = 0; j < froBot.size(); ++j) { ret(i, 0) = froBot(j); ret(i, 1) = froTop(j); ++i; } for (size_t j = 0; j < botLft.size(); ++j) { ret(i, 0) = botLft(j); ret(i, 1) = botRgt(j); ++i; } for (size_t j = 0; j < botLft.size(); ++j) { ret(i, 0) = botLft(j); ret(i, 1) = topRgt(j); ++i; } for (size_t j = 0; j < botLft.size(); ++j) { ret(i, 0) = botLft(j); ret(i, 1) = topLft(j); ++i; } for (size_t j = 0; j < froLft.size(); ++j) { ret(i, 0) = froLft(j); ret(i, 1) = froRgt(j); ++i; } for (size_t j = 0; j < froLft.size(); ++j) { ret(i, 0) = froLft(j); ret(i, 1) = bckRgt(j); ++i; } for (size_t j = 0; j < froLft.size(); ++j) { ret(i, 0) = froLft(j); ret(i, 1) = bckLft(j); ++i; } for (size_t j = 0; j < fro.size(); ++j) { ret(i, 0) = fro(j); ret(i, 1) = bck(j); ++i; } for (size_t j = 0; j < lft.size(); ++j) { ret(i, 0) = lft(j); ret(i, 1) = rgt(j); ++i; } for (size_t j = 0; j < bot.size(); ++j) { ret(i, 0) = bot(j); ret(i, 1) = top(j); ++i; } return ret; } template inline auto RegularBase3d::nodesOrigin_impl() const { return derived_cast().nodesFrontBottomLeftCorner_impl(); } template inline auto RegularBase3d::derived_cast() -> derived_type& { return *static_cast(this); } template inline auto RegularBase3d::derived_cast() const -> const derived_type& { return *static_cast(this); } namespace detail { template inline T renum(const T& arg, const R& mapping) { T ret = T::from_shape(arg.shape()); auto jt = ret.begin(); for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) { *jt = mapping(*it); } return ret; } } // namespace detail template inline xt::xtensor overlapping(const S& coor_a, const T& coor_b, double rtol, double atol) { GOOSEFEM_ASSERT(coor_a.dimension() == 2); GOOSEFEM_ASSERT(coor_b.dimension() == 2); GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); std::vector ret_a; std::vector ret_b; for (size_t i = 0; i < coor_a.shape(0); ++i) { auto idx = xt::flatten_indices(xt::argwhere( xt::prod(xt::isclose(coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1))); for (auto& j : idx) { ret_a.push_back(i); ret_b.push_back(j); } } xt::xtensor ret = xt::empty({size_t(2), ret_a.size()}); for (size_t i = 0; i < ret_a.size(); ++i) { ret(0, i) = ret_a[i]; ret(1, i) = ret_b[i]; } return ret; } template inline ManualStitch::ManualStitch( const CA& coor_a, const EA& conn_a, const NA& overlapping_nodes_a, const CB& coor_b, const EB& conn_b, const NB& overlapping_nodes_b, bool check_position, double rtol, double atol) { UNUSED(rtol); UNUSED(atol); GOOSEFEM_ASSERT(coor_a.dimension() == 2); GOOSEFEM_ASSERT(conn_a.dimension() == 2); GOOSEFEM_ASSERT(overlapping_nodes_a.dimension() == 1); GOOSEFEM_ASSERT(coor_b.dimension() == 2); GOOSEFEM_ASSERT(conn_b.dimension() == 2); GOOSEFEM_ASSERT(overlapping_nodes_b.dimension() == 1); GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape())); GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1)); GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1)); if (check_position) { GOOSEFEM_CHECK(xt::allclose( xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()), xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()), rtol, atol)); } size_t nnda = coor_a.shape(0); size_t nndb = coor_b.shape(0); size_t ndim = coor_a.shape(1); size_t nelim = overlapping_nodes_a.size(); size_t nela = conn_a.shape(0); size_t nelb = conn_b.shape(0); size_t nne = conn_a.shape(1); m_nel_a = nela; m_nel_b = nelb; m_nnd_a = nnda; xt::xtensor keep_b = xt::setdiff1d(xt::arange(nndb), overlapping_nodes_b); m_map_b = xt::empty({nndb}); xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a; xt::view(m_map_b, xt::keep(keep_b)) = xt::arange(keep_b.size()) + nnda; m_conn = xt::empty({nela + nelb, nne}); xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a; xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b); m_coor = xt::empty({nnda + nndb - nelim, ndim}); xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a; xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) = xt::view(coor_b, xt::keep(keep_b), xt::all()); } inline xt::xtensor ManualStitch::coor() const { return m_coor; } inline xt::xtensor ManualStitch::conn() const { return m_conn; } inline size_t ManualStitch::nmesh() const { return 2; } inline size_t ManualStitch::nelem() const { return m_conn.shape(0); } inline size_t ManualStitch::nnode() const { return m_coor.shape(0); } inline size_t ManualStitch::nne() const { return m_conn.shape(1); } inline size_t ManualStitch::ndim() const { return m_coor.shape(1); } inline xt::xtensor ManualStitch::dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline std::vector> ManualStitch::nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } inline std::vector> ManualStitch::elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } inline xt::xtensor ManualStitch::nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nnd_a); } return m_map_b; } inline xt::xtensor ManualStitch::elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { return xt::arange(m_nel_a); } return xt::arange(m_nel_b) + m_nel_a; } template inline T ManualStitch::nodeset(const T& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a); return set; } GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size()); return detail::renum(set, m_map_b); } template inline T ManualStitch::elemset(const T& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index <= 1); if (mesh_index == 0) { GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a); return set; } GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b); return set + m_nel_a; } inline Stitch::Stitch(double rtol, double atol) { m_rtol = rtol; m_atol = atol; } template inline void Stitch::push_back(const C& coor, const E& conn) { GOOSEFEM_ASSERT(coor.dimension() == 2); GOOSEFEM_ASSERT(conn.dimension() == 2); if (m_map.size() == 0) { m_coor = coor; m_conn = conn; m_map.push_back(xt::eval(xt::arange(coor.shape(0)))); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(0); return; } auto overlap = overlapping(m_coor, coor, m_rtol, m_atol); size_t index = m_map.size(); ManualStitch stitch( m_coor, m_conn, xt::eval(xt::view(overlap, 0, xt::all())), coor, conn, xt::eval(xt::view(overlap, 1, xt::all())), false); m_coor = stitch.coor(); m_conn = stitch.conn(); m_map.push_back(stitch.nodemap(1)); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]); } inline size_t Stitch::nmesh() const { return m_map.size(); } inline xt::xtensor Stitch::coor() const { return m_coor; } inline xt::xtensor Stitch::conn() const { return m_conn; } inline size_t Stitch::nelem() const { return m_conn.shape(0); } inline size_t Stitch::nnode() const { return m_coor.shape(0); } inline size_t Stitch::nne() const { return m_conn.shape(1); } inline size_t Stitch::ndim() const { return m_coor.shape(1); } inline xt::xtensor Stitch::dofs() const { size_t nnode = this->nnode(); size_t ndim = this->ndim(); return xt::reshape_view(xt::arange(nnode * ndim), {nnode, ndim}); } inline std::vector> Stitch::nodemap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->nodemap(i); } return ret; } inline std::vector> Stitch::elemmap() const { std::vector> ret(this->nmesh()); for (size_t i = 0; i < this->nmesh(); ++i) { ret[i] = this->elemmap(i); } return ret; } inline xt::xtensor Stitch::nodemap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return m_map[mesh_index]; } inline xt::xtensor Stitch::elemmap(size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); return xt::arange(m_nel[mesh_index]) + m_el_offset[mesh_index]; } template inline T Stitch::nodeset(const T& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size()); return detail::renum(set, m_map[mesh_index]); } template inline T Stitch::elemset(const T& set, size_t mesh_index) const { GOOSEFEM_ASSERT(mesh_index < m_map.size()); GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]); return set + m_el_offset[mesh_index]; } template inline T Stitch::nodeset(const std::vector& set) const { GOOSEFEM_ASSERT(set.size() == m_map.size()); size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->nodeset(set[i], i); n += set[i].size(); } return xt::unique(ret); } template inline T Stitch::elemset(const std::vector& set) const { GOOSEFEM_ASSERT(set.size() == m_map.size()); size_t n = 0; for (size_t i = 0; i < set.size(); ++i) { n += set[i].size(); } xt::xtensor ret = xt::empty({n}); n = 0; for (size_t i = 0; i < set.size(); ++i) { xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i); n += set[i].size(); } return ret; } template inline T Stitch::nodeset(std::initializer_list set) const { return this->nodeset(std::vector(set)); } template inline T Stitch::elemset(std::initializer_list set) const { return this->elemset(std::vector(set)); } inline Vstack::Vstack(bool check_overlap, double rtol, double atol) { m_check_overlap = check_overlap; m_rtol = rtol; m_atol = atol; } template inline void Vstack::push_back(const C& coor, const E& conn, const N& nodes_bot, const N& nodes_top) { if (m_map.size() == 0) { m_coor = coor; m_conn = conn; m_map.push_back(xt::eval(xt::arange(coor.shape(0)))); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(0); m_nodes_bot.push_back(nodes_bot); m_nodes_top.push_back(nodes_top); return; } GOOSEFEM_ASSERT(nodes_bot.size() == m_nodes_top.back().size()); size_t index = m_map.size(); double shift = xt::amax(xt::view(m_coor, xt::all(), 1))(); auto x = coor; xt::view(x, xt::all(), 1) += shift; ManualStitch stitch( m_coor, m_conn, m_nodes_top.back(), x, conn, nodes_bot, m_check_overlap, m_rtol, m_atol); m_nodes_bot.push_back(stitch.nodeset(nodes_bot, 1)); m_nodes_top.push_back(stitch.nodeset(nodes_top, 1)); m_coor = stitch.coor(); m_conn = stitch.conn(); m_map.push_back(stitch.nodemap(1)); m_nel.push_back(conn.shape(0)); m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]); } template inline Renumber::Renumber(const T& dofs) { size_t n = xt::amax(dofs)() + 1; size_t i = 0; xt::xtensor unique = xt::unique(dofs); m_renum = xt::empty({n}); for (auto& j : unique) { m_renum(j) = i; ++i; } } template inline T Renumber::apply(const T& list) const { return detail::renum(list, m_renum); } inline xt::xtensor Renumber::index() const { return m_renum; } template inline T renumber(const T& dofs) { return Renumber(dofs).apply(dofs); } template inline Reorder::Reorder(const std::initializer_list args) { size_t n = 0; size_t i = 0; for (auto& arg : args) { if (arg.size() == 0) { continue; } n = std::max(n, xt::amax(arg)() + 1); } #ifdef GOOSEFEM_ENABLE_ASSERT for (auto& arg : args) { GOOSEFEM_ASSERT(is_unique(arg)); } #endif m_renum = xt::empty({n}); for (auto& arg : args) { for (auto& j : arg) { m_renum(j) = i; ++i; } } } 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}); } template inline xt::xtensor coordination(const E& conn) { GOOSEFEM_ASSERT(conn.dimension() == 2); size_t nnode = xt::amax(conn)() + 1; xt::xtensor N = xt::zeros({nnode}); for (auto it = conn.begin(); it != conn.end(); ++it) { N(*it) += 1; } return N; } template inline std::vector> elem2node(const E& 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; } template inline xt::xtensor edgesize(const C& coor, const E& conn, ElementType type) { GOOSEFEM_ASSERT(coor.dimension() == 2); GOOSEFEM_ASSERT(conn.dimension() == 2); 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"); } template inline xt::xtensor edgesize(const C& coor, const E& conn) { return edgesize(coor, conn, defaultElementType(coor, conn)); } template inline xt::xtensor centers(const C& coor, const E& conn, ElementType type) { GOOSEFEM_ASSERT(coor.dimension() == 2); GOOSEFEM_ASSERT(conn.dimension() == 2); 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"); } template inline xt::xtensor centers(const C& coor, const E& conn) { return centers(coor, conn, defaultElementType(coor, conn)); } template inline xt::xtensor elemmap2nodemap(const T& elem_map, const C& coor, const E& conn, ElementType type) { GOOSEFEM_ASSERT(elem_map.dimension() == 1); GOOSEFEM_ASSERT(coor.dimension() == 2); GOOSEFEM_ASSERT(conn.dimension() == 2); 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"); } template inline xt::xtensor elemmap2nodemap(const T& elem_map, const C& coor, const E& conn) { return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn)); } +template +inline xt::xtensor center_of_gravity(const C& coor, const E& conn, ElementType type) +{ + auto dof = dofs(coor.shape(0), coor.shape(1)); + + GooseFEM::MatrixDiagonal M(conn, dof); + + GooseFEM::Vector vector(conn, dof); + + GooseFEM::Element::Quad4::Quadrature quad( + vector.AsElement(coor), + GooseFEM::Element::Quad4::Nodal::xi(), + GooseFEM::Element::Quad4::Nodal::w()); + + xt::xtensor rho = xt::ones({conn.shape(0), quad.nip()}); + + M.assemble(quad.Int_N_scalar_NT_dV(rho)); + + auto m = vector.AsNode(M.Todiagonal()); + return xt::average(coor, m, 0); +} + +template +inline xt::xtensor center_of_gravity(const C& coor, const E& conn) +{ + return center_of_gravity(coor, conn, defaultElementType(coor, conn)); +} + } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/Vector.hpp b/include/GooseFEM/Vector.hpp index 526ca7d..a54e174 100644 --- a/include/GooseFEM/Vector.hpp +++ b/include/GooseFEM/Vector.hpp @@ -1,553 +1,550 @@ /** Implementation of Vector.h \file Vector.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_VECTOR_HPP #define GOOSEFEM_VECTOR_HPP #include "Vector.h" namespace GooseFEM { template -inline Vector::Vector(const S& conn, const T& dofs) +inline Vector::Vector(const S& conn, const T& dofs) : m_conn(conn), m_dofs(dofs) { GOOSEFEM_ASSERT(conn.dimension() == 2); GOOSEFEM_ASSERT(dofs.dimension() == 2); - m_conn = conn; - m_dofs = dofs; - m_nelem = m_conn.shape(0); m_nne = m_conn.shape(1); m_nnode = m_dofs.shape(0); m_ndim = m_dofs.shape(1); m_ndof = xt::amax(m_dofs)() + 1; GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); } inline size_t Vector::nelem() const { return m_nelem; } inline size_t Vector::nne() const { return m_nne; } inline size_t Vector::nnode() const { return m_nnode; } inline size_t Vector::ndim() const { return m_ndim; } inline size_t Vector::ndof() const { return m_ndof; } inline xt::xtensor Vector::conn() const { return m_conn; } inline xt::xtensor Vector::dofs() const { return m_dofs; } template inline T Vector::Copy(const T& nodevec_src, const T& nodevec_dest) const { T ret = T::from_shape(nodevec_dest.shape()); this->copy(nodevec_src, ret); return ret; } template inline void Vector::copy(const T& nodevec_src, T& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, this->shape_nodevec())); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, this->shape_nodevec())); xt::noalias(nodevec_dest) = nodevec_src; } // asDofs template inline xt::xtensor Vector::AsDofs(const T& arg) const { xt::xtensor ret = xt::empty(this->shape_dofval()); this->asDofs_impl(arg, ret); return ret; } template void Vector::asDofs(const T& arg, R& ret) const { this->asDofs_impl(arg, ret); } // asDofs : implementation distribution template ::value, int>> inline void Vector::asDofs_impl(const T& arg, R& ret) const { if (arg.dimension() == 2) { this->asDofs_impl_nodevec(arg, ret); } else if (arg.dimension() == 3) { this->asDofs_impl_elemvec(arg, ret); } else { throw std::runtime_error("Vector::asDofs unknown dimension for conversion"); } } template ::value == 2, int>> inline void Vector::asDofs_impl(const T& arg, R& ret) const { this->asDofs_impl_nodevec(arg, ret); } template ::value == 3, int>> inline void Vector::asDofs_impl(const T& arg, R& ret) const { this->asDofs_impl_elemvec(arg, ret); } // asDofs : implementation template inline void Vector::asDofs_impl_nodevec(const T& arg, R& dofval) const { static_assert( xt::get_rank::value == 1 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(arg, this->shape_nodevec())); GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval())); dofval.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m, i)) = arg(m, i); } } } template inline void Vector::asDofs_impl_elemvec(const T& arg, R& dofval) const { static_assert( xt::get_rank::value == 1 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(arg, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval())); dofval.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m_conn(e, m), i)) = arg(e, m, i); } } } } // asNode template inline xt::xtensor Vector::AsNode(const T& arg) const { xt::xtensor ret = xt::empty(this->shape_nodevec()); this->asNode_impl(arg, ret); return ret; } template void Vector::asNode(const T& arg, R& ret) const { this->asNode_impl(arg, ret); } // asNode : implementation distribution template ::value, int>> inline void Vector::asNode_impl(const T& arg, R& ret) const { if (arg.dimension() == 1) { this->asNode_impl_dofval(arg, ret); } else if (arg.dimension() == 3) { this->asNode_impl_elemvec(arg, ret); } else { throw std::runtime_error("Vector::asNode unknown dimension for conversion"); } } template ::value == 1, int>> inline void Vector::asNode_impl(const T& arg, R& ret) const { this->asNode_impl_dofval(arg, ret); } template ::value == 3, int>> inline void Vector::asNode_impl(const T& arg, R& ret) const { this->asNode_impl_elemvec(arg, ret); } // asNode : implementation template inline void Vector::asNode_impl_dofval(const T& dofval, R& nodevec) const { static_assert( xt::get_rank::value == 2 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval())); GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec())); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { nodevec(m, i) = dofval(m_dofs(m, i)); } } } template inline void Vector::asNode_impl_elemvec(const T& elemvec, R& nodevec) const { static_assert( xt::get_rank::value == 2 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec())); nodevec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { nodevec(m_conn(e, m), i) = elemvec(e, m, i); } } } } // asElement template inline xt::xtensor Vector::AsElement(const T& arg) const { xt::xtensor ret = xt::empty(this->shape_elemvec()); this->asElement_impl(arg, ret); return ret; } template void Vector::asElement(const T& arg, R& ret) const { this->asElement_impl(arg, ret); } // asElement : implementation distribution template ::value, int>> inline void Vector::asElement_impl(const T& arg, R& ret) const { if (arg.dimension() == 1) { this->asElement_impl_dofval(arg, ret); } else if (arg.dimension() == 2) { this->asElement_impl_nodevec(arg, ret); } else { throw std::runtime_error("Vector::asElement unknown dimension for conversion"); } } template ::value == 1, int>> inline void Vector::asElement_impl(const T& arg, R& ret) const { this->asElement_impl_dofval(arg, ret); } template ::value == 2, int>> inline void Vector::asElement_impl(const T& arg, R& ret) const { this->asElement_impl_nodevec(arg, ret); } // asElement : implementation template inline void Vector::asElement_impl_dofval(const T& dofval, R& elemvec) const { static_assert( xt::get_rank::value == 3 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { elemvec(e, m, i) = dofval(m_dofs(m_conn(e, m), i)); } } } } template inline void Vector::asElement_impl_nodevec(const T& nodevec, R& elemvec) const { static_assert( xt::get_rank::value == 3 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec())); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { elemvec(e, m, i) = nodevec(m_conn(e, m), i); } } } } // assembleDofs template inline xt::xtensor Vector::AssembleDofs(const T& arg) const { xt::xtensor ret = xt::empty(this->shape_dofval()); this->assembleDofs_impl(arg, ret); return ret; } template void Vector::assembleDofs(const T& arg, R& ret) const { this->assembleDofs_impl(arg, ret); } // assembleDofs : implementation distribution template ::value, int>> inline void Vector::assembleDofs_impl(const T& arg, R& ret) const { if (arg.dimension() == 2) { this->assembleDofs_impl_nodevec(arg, ret); } else if (arg.dimension() == 3) { this->assembleDofs_impl_elemvec(arg, ret); } else { throw std::runtime_error("Vector::assembleDofs unknown dimension for conversion"); } } template ::value == 2, int>> inline void Vector::assembleDofs_impl(const T& arg, R& ret) const { this->assembleDofs_impl_nodevec(arg, ret); } template ::value == 3, int>> inline void Vector::assembleDofs_impl(const T& arg, R& ret) const { this->assembleDofs_impl_elemvec(arg, ret); } // asDofs : implementation template inline void Vector::assembleDofs_impl_nodevec(const T& nodevec, R& dofval) const { static_assert( xt::get_rank::value == 1 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec())); GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval())); dofval.fill(0.0); for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m, i)) += nodevec(m, i); } } } template inline void Vector::assembleDofs_impl_elemvec(const T& elemvec, R& dofval) const { static_assert( xt::get_rank::value == 1 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval())); dofval.fill(0.0); for (size_t e = 0; e < m_nelem; ++e) { for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m_conn(e, m), i)) += elemvec(e, m, i); } } } } // assembleNode template inline xt::xtensor Vector::AssembleNode(const T& arg) const { xt::xtensor ret = xt::empty(this->shape_nodevec()); this->assembleNode_impl(arg, ret); return ret; } template void Vector::assembleNode(const T& arg, R& ret) const { this->assembleNode_impl(arg, ret); } // assembleNode : implementation distribution template ::value, int>> inline void Vector::assembleNode_impl(const T& arg, R& ret) const { if (arg.dimension() == 3) { this->assembleNode_impl_elemvec(arg, ret); } else { throw std::runtime_error("Vector::assembleNode unknown dimension for conversion"); } } template ::value == 3, int>> inline void Vector::assembleNode_impl(const T& arg, R& ret) const { this->assembleNode_impl_elemvec(arg, ret); } // asNode : implementation template inline void Vector::assembleNode_impl_elemvec(const T& elemvec, R& nodevec) const { static_assert( xt::get_rank::value == 2 || !xt::has_fixed_rank_t::value, "Unknown rank 'ret'"); GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec())); GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec())); xt::xtensor dofval = this->AssembleDofs(elemvec); this->asNode(dofval, nodevec); } inline std::array Vector::shape_dofval() const { std::array shape; shape[0] = m_ndof; return shape; } inline std::array Vector::shape_nodevec() const { std::array shape; shape[0] = m_nnode; shape[1] = m_ndim; return shape; } inline std::array Vector::shape_elemvec() const { std::array shape; shape[0] = m_nelem; shape[1] = m_nne; shape[2] = m_ndim; return shape; } inline std::array Vector::shape_elemmat() const { std::array shape; shape[0] = m_nelem; shape[1] = m_nne * m_ndim; shape[2] = m_nne * m_ndim; return shape; } inline xt::xtensor Vector::allocate_dofval() const { xt::xtensor dofval = xt::empty(this->shape_dofval()); return dofval; } inline xt::xtensor Vector::allocate_nodevec() const { xt::xtensor nodevec = xt::empty(this->shape_nodevec()); return nodevec; } inline xt::xtensor Vector::allocate_elemvec() const { xt::xtensor elemvec = xt::empty(this->shape_elemvec()); return elemvec; } inline xt::xtensor Vector::allocate_elemmat() const { xt::xtensor elemmat = xt::empty(this->shape_elemmat()); return elemmat; } inline xt::xtensor Vector::allocate_dofval(double val) const { xt::xtensor dofval = xt::empty(this->shape_dofval()); dofval.fill(val); return dofval; } inline xt::xtensor Vector::allocate_nodevec(double val) const { xt::xtensor nodevec = xt::empty(this->shape_nodevec()); nodevec.fill(val); return nodevec; } inline xt::xtensor Vector::allocate_elemvec(double val) const { xt::xtensor elemvec = xt::empty(this->shape_elemvec()); elemvec.fill(val); return elemvec; } inline xt::xtensor Vector::allocate_elemmat(double val) const { xt::xtensor elemmat = xt::empty(this->shape_elemmat()); elemmat.fill(val); return elemmat; } } // namespace GooseFEM #endif diff --git a/python/Mesh.hpp b/python/Mesh.hpp index f473bdd..9ad500f 100644 --- a/python/Mesh.hpp +++ b/python/Mesh.hpp @@ -1,657 +1,683 @@ /** \file \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef PYGOOSEFEM_MESH_H #define PYGOOSEFEM_MESH_H #include #include #include #include #include namespace py = pybind11; template void register_Element_RegularBase(P& cls) { cls.def("nelem", &C::nelem); cls.def("nnode", &C::nnode); cls.def("nne", &C::nne); cls.def("ndim", &C::ndim); cls.def("nelx", &C::nelx); cls.def("nely", &C::nely); cls.def("h", &C::h); cls.def("getElementType", &C::getElementType); cls.def("coor", &C::coor); cls.def("conn", &C::conn); cls.def("dofs", &C::dofs); cls.def("dofsPeriodic", &C::dofsPeriodic); cls.def("nodesPeriodic", &C::nodesPeriodic); cls.def("nodesOrigin", &C::nodesOrigin); } template void register_Element_RegularBase2d(P& cls) { cls.def("nodesBottomEdge", &C::nodesBottomEdge); cls.def("nodesTopEdge", &C::nodesTopEdge); cls.def("nodesLeftEdge", &C::nodesLeftEdge); cls.def("nodesRightEdge", &C::nodesRightEdge); cls.def("nodesBottomOpenEdge", &C::nodesBottomOpenEdge); cls.def("nodesTopOpenEdge", &C::nodesTopOpenEdge); cls.def("nodesLeftOpenEdge", &C::nodesLeftOpenEdge); cls.def("nodesRightOpenEdge", &C::nodesRightOpenEdge); cls.def("nodesBottomLeftCorner", &C::nodesBottomLeftCorner); cls.def("nodesBottomRightCorner", &C::nodesBottomRightCorner); cls.def("nodesTopLeftCorner", &C::nodesTopLeftCorner); cls.def("nodesTopRightCorner", &C::nodesTopRightCorner); cls.def("nodesLeftBottomCorner", &C::nodesLeftBottomCorner); cls.def("nodesLeftTopCorner", &C::nodesLeftTopCorner); cls.def("nodesRightBottomCorner", &C::nodesRightBottomCorner); cls.def("nodesRightTopCorner", &C::nodesRightTopCorner); } template void register_Element_RegularBase3d(P& cls) { cls.def("nodesFront", &C::nodesFront); cls.def("nodesBack", &C::nodesBack); cls.def("nodesLeft", &C::nodesLeft); cls.def("nodesRight", &C::nodesRight); cls.def("nodesBottom", &C::nodesBottom); cls.def("nodesTop", &C::nodesTop); cls.def("nodesFrontFace", &C::nodesFrontFace); cls.def("nodesBackFace", &C::nodesBackFace); cls.def("nodesLeftFace", &C::nodesLeftFace); cls.def("nodesRightFace", &C::nodesRightFace); cls.def("nodesBottomFace", &C::nodesBottomFace); cls.def("nodesTopFace", &C::nodesTopFace); cls.def("nodesFrontBottomEdge", &C::nodesFrontBottomEdge); cls.def("nodesFrontTopEdge", &C::nodesFrontTopEdge); cls.def("nodesFrontLeftEdge", &C::nodesFrontLeftEdge); cls.def("nodesFrontRightEdge", &C::nodesFrontRightEdge); cls.def("nodesBackBottomEdge", &C::nodesBackBottomEdge); cls.def("nodesBackTopEdge", &C::nodesBackTopEdge); cls.def("nodesBackLeftEdge", &C::nodesBackLeftEdge); cls.def("nodesBackRightEdge", &C::nodesBackRightEdge); cls.def("nodesBottomLeftEdge", &C::nodesBottomLeftEdge); cls.def("nodesBottomRightEdge", &C::nodesBottomRightEdge); cls.def("nodesTopLeftEdge", &C::nodesTopLeftEdge); cls.def("nodesTopRightEdge", &C::nodesTopRightEdge); cls.def("nodesBottomFrontEdge", &C::nodesBottomFrontEdge); cls.def("nodesBottomBackEdge", &C::nodesBottomBackEdge); cls.def("nodesTopFrontEdge", &C::nodesTopFrontEdge); cls.def("nodesTopBackEdge", &C::nodesTopBackEdge); cls.def("nodesLeftBottomEdge", &C::nodesLeftBottomEdge); cls.def("nodesLeftFrontEdge", &C::nodesLeftFrontEdge); cls.def("nodesLeftBackEdge", &C::nodesLeftBackEdge); cls.def("nodesLeftTopEdge", &C::nodesLeftTopEdge); cls.def("nodesRightBottomEdge", &C::nodesRightBottomEdge); cls.def("nodesRightTopEdge", &C::nodesRightTopEdge); cls.def("nodesRightFrontEdge", &C::nodesRightFrontEdge); cls.def("nodesRightBackEdge", &C::nodesRightBackEdge); cls.def("nodesFrontBottomOpenEdge", &C::nodesFrontBottomOpenEdge); cls.def("nodesFrontTopOpenEdge", &C::nodesFrontTopOpenEdge); cls.def("nodesFrontLeftOpenEdge", &C::nodesFrontLeftOpenEdge); cls.def("nodesFrontRightOpenEdge", &C::nodesFrontRightOpenEdge); cls.def("nodesBackBottomOpenEdge", &C::nodesBackBottomOpenEdge); cls.def("nodesBackTopOpenEdge", &C::nodesBackTopOpenEdge); cls.def("nodesBackLeftOpenEdge", &C::nodesBackLeftOpenEdge); cls.def("nodesBackRightOpenEdge", &C::nodesBackRightOpenEdge); cls.def("nodesBottomLeftOpenEdge", &C::nodesBottomLeftOpenEdge); cls.def("nodesBottomRightOpenEdge", &C::nodesBottomRightOpenEdge); cls.def("nodesTopLeftOpenEdge", &C::nodesTopLeftOpenEdge); cls.def("nodesTopRightOpenEdge", &C::nodesTopRightOpenEdge); cls.def("nodesBottomFrontOpenEdge", &C::nodesBottomFrontOpenEdge); cls.def("nodesBottomBackOpenEdge", &C::nodesBottomBackOpenEdge); cls.def("nodesTopFrontOpenEdge", &C::nodesTopFrontOpenEdge); cls.def("nodesTopBackOpenEdge", &C::nodesTopBackOpenEdge); cls.def("nodesLeftBottomOpenEdge", &C::nodesLeftBottomOpenEdge); cls.def("nodesLeftFrontOpenEdge", &C::nodesLeftFrontOpenEdge); cls.def("nodesLeftBackOpenEdge", &C::nodesLeftBackOpenEdge); cls.def("nodesLeftTopOpenEdge", &C::nodesLeftTopOpenEdge); cls.def("nodesRightBottomOpenEdge", &C::nodesRightBottomOpenEdge); cls.def("nodesRightTopOpenEdge", &C::nodesRightTopOpenEdge); cls.def("nodesRightFrontOpenEdge", &C::nodesRightFrontOpenEdge); cls.def("nodesRightBackOpenEdge", &C::nodesRightBackOpenEdge); cls.def("nodesFrontBottomLeftCorner", &C::nodesFrontBottomLeftCorner); cls.def("nodesFrontBottomRightCorner", &C::nodesFrontBottomRightCorner); cls.def("nodesFrontTopLeftCorner", &C::nodesFrontTopLeftCorner); cls.def("nodesFrontTopRightCorner", &C::nodesFrontTopRightCorner); cls.def("nodesBackBottomLeftCorner", &C::nodesBackBottomLeftCorner); cls.def("nodesBackBottomRightCorner", &C::nodesBackBottomRightCorner); cls.def("nodesBackTopLeftCorner", &C::nodesBackTopLeftCorner); cls.def("nodesBackTopRightCorner", &C::nodesBackTopRightCorner); cls.def("nodesFrontLeftBottomCorner", &C::nodesFrontLeftBottomCorner); cls.def("nodesBottomFrontLeftCorner", &C::nodesBottomFrontLeftCorner); cls.def("nodesBottomLeftFrontCorner", &C::nodesBottomLeftFrontCorner); cls.def("nodesLeftFrontBottomCorner", &C::nodesLeftFrontBottomCorner); cls.def("nodesLeftBottomFrontCorner", &C::nodesLeftBottomFrontCorner); cls.def("nodesFrontRightBottomCorner", &C::nodesFrontRightBottomCorner); cls.def("nodesBottomFrontRightCorner", &C::nodesBottomFrontRightCorner); cls.def("nodesBottomRightFrontCorner", &C::nodesBottomRightFrontCorner); cls.def("nodesRightFrontBottomCorner", &C::nodesRightFrontBottomCorner); cls.def("nodesRightBottomFrontCorner", &C::nodesRightBottomFrontCorner); cls.def("nodesFrontLeftTopCorner", &C::nodesFrontLeftTopCorner); cls.def("nodesTopFrontLeftCorner", &C::nodesTopFrontLeftCorner); cls.def("nodesTopLeftFrontCorner", &C::nodesTopLeftFrontCorner); cls.def("nodesLeftFrontTopCorner", &C::nodesLeftFrontTopCorner); cls.def("nodesLeftTopFrontCorner", &C::nodesLeftTopFrontCorner); cls.def("nodesFrontRightTopCorner", &C::nodesFrontRightTopCorner); cls.def("nodesTopFrontRightCorner", &C::nodesTopFrontRightCorner); cls.def("nodesTopRightFrontCorner", &C::nodesTopRightFrontCorner); cls.def("nodesRightFrontTopCorner", &C::nodesRightFrontTopCorner); cls.def("nodesRightTopFrontCorner", &C::nodesRightTopFrontCorner); cls.def("nodesBackLeftBottomCorner", &C::nodesBackLeftBottomCorner); cls.def("nodesBottomBackLeftCorner", &C::nodesBottomBackLeftCorner); cls.def("nodesBottomLeftBackCorner", &C::nodesBottomLeftBackCorner); cls.def("nodesLeftBackBottomCorner", &C::nodesLeftBackBottomCorner); cls.def("nodesLeftBottomBackCorner", &C::nodesLeftBottomBackCorner); cls.def("nodesBackRightBottomCorner", &C::nodesBackRightBottomCorner); cls.def("nodesBottomBackRightCorner", &C::nodesBottomBackRightCorner); cls.def("nodesBottomRightBackCorner", &C::nodesBottomRightBackCorner); cls.def("nodesRightBackBottomCorner", &C::nodesRightBackBottomCorner); cls.def("nodesRightBottomBackCorner", &C::nodesRightBottomBackCorner); cls.def("nodesBackLeftTopCorner", &C::nodesBackLeftTopCorner); cls.def("nodesTopBackLeftCorner", &C::nodesTopBackLeftCorner); cls.def("nodesTopLeftBackCorner", &C::nodesTopLeftBackCorner); cls.def("nodesLeftBackTopCorner", &C::nodesLeftBackTopCorner); cls.def("nodesLeftTopBackCorner", &C::nodesLeftTopBackCorner); cls.def("nodesBackRightTopCorner", &C::nodesBackRightTopCorner); cls.def("nodesTopBackRightCorner", &C::nodesTopBackRightCorner); cls.def("nodesTopRightBackCorner", &C::nodesTopRightBackCorner); cls.def("nodesRightBackTopCorner", &C::nodesRightBackTopCorner); cls.def("nodesRightTopBackCorner", &C::nodesRightTopBackCorner); } void init_Mesh(py::module& mod) { py::enum_( mod, "ElementType", "ElementType." "See :cpp:enum:`GooseFEM::Mesh::ElementType`.") .value("Unknown", GooseFEM::Mesh::ElementType::Unknown) .value("Tri3", GooseFEM::Mesh::ElementType::Tri3) .value("Quad4", GooseFEM::Mesh::ElementType::Quad4) .value("Hex8", GooseFEM::Mesh::ElementType::Hex8) .export_values(); mod.def( "overlapping", &GooseFEM::Mesh::overlapping, xt::pytensor>, "Find overlapping nodes." "See :cpp:func:`GooseFEM::Mesh::overlapping`.", py::arg("coor_a"), py::arg("coor_b"), py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8); py::class_(mod, "ManualStitch") .def( py::init< const xt::pytensor&, const xt::pytensor&, const xt::pytensor&, const xt::pytensor&, const xt::pytensor&, const xt::pytensor&, bool, double, double>(), "Manually stitch meshes." "See :cpp:class:`GooseFEM::Mesh::ManualStitch`.", py::arg("coor_a"), py::arg("conn_a"), py::arg("overlapping_nodes_a"), py::arg("coor_b"), py::arg("conn_b"), py::arg("overlapping_nodes_b"), py::arg("check_position") = true, py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8) .def( "nmesh", &GooseFEM::Mesh::ManualStitch::nmesh, "Number of sub meshes." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nmesh`.") .def( "nnode", &GooseFEM::Mesh::ManualStitch::nnode, "Number of nodes of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nnode`.") .def( "nelem", &GooseFEM::Mesh::ManualStitch::nelem, "Number of elements of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nelem`.") .def( "ndim", &GooseFEM::Mesh::ManualStitch::ndim, "Number of dimensions (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::ndim`.") .def( "nne", &GooseFEM::Mesh::ManualStitch::nne, "Number of nodes per element (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nne`.") .def( "coor", &GooseFEM::Mesh::ManualStitch::coor, "Coordinates of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::coor`.") .def( "conn", &GooseFEM::Mesh::ManualStitch::conn, "Connectivity of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::conn`.") .def( "dofs", &GooseFEM::Mesh::ManualStitch::dofs, "DOF numbers per node." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::dofs`.") .def( "nodemap", py::overload_cast<>(&GooseFEM::Mesh::ManualStitch::nodemap, py::const_), "Node-map for givall sub-meshes." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.") .def( "elemmap", py::overload_cast<>(&GooseFEM::Mesh::ManualStitch::elemmap, py::const_), "Element-map for all sub-meshes." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemmap`.") .def( "nodemap", py::overload_cast(&GooseFEM::Mesh::ManualStitch::nodemap, py::const_), "Node-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.", py::arg("mesh_index")) .def( "elemmap", py::overload_cast(&GooseFEM::Mesh::ManualStitch::elemmap, py::const_), "Element-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemmap`.", py::arg("mesh_index")) .def( "nodeset", &GooseFEM::Mesh::ManualStitch::nodeset>, "Convert node-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodeset`.", py::arg("set"), py::arg("mesh_index")) .def( "elemset", &GooseFEM::Mesh::ManualStitch::elemset>, "Convert element-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemset`.", py::arg("set"), py::arg("mesh_index")) .def("__repr__", [](const GooseFEM::Mesh::ManualStitch&) { return ""; }); py::class_(mod, "Stitch") .def( py::init(), "Stitch meshes." "See :cpp:class:`GooseFEM::Mesh::Stitch`.", py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8) .def( "push_back", &GooseFEM::Mesh::Stitch::push_back, xt::pytensor>, "Add mesh to be stitched." "See :cpp:func:`GooseFEM::Mesh::Stitch::push_back`.", py::arg("coor"), py::arg("conn")) .def( "nmesh", &GooseFEM::Mesh::Stitch::nmesh, "Number of sub meshes." "See :cpp:func:`GooseFEM::Mesh::Stitch::nmesh`.") .def( "nnode", &GooseFEM::Mesh::Stitch::nnode, "Number of nodes of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nnode`.") .def( "nelem", &GooseFEM::Mesh::Stitch::nelem, "Number of elements of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nelem`.") .def( "ndim", &GooseFEM::Mesh::Stitch::ndim, "Number of dimensions (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::Stitch::ndim`.") .def( "nne", &GooseFEM::Mesh::Stitch::nne, "Number of nodes per element (of stitched mesh)." "See :cpp:func:`GooseFEM::Mesh::Stitch::nne`.") .def( "coor", &GooseFEM::Mesh::Stitch::coor, "Coordinates of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::coor`.") .def( "conn", &GooseFEM::Mesh::Stitch::conn, "Connectivity of stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::conn`.") .def( "dofs", &GooseFEM::Mesh::Stitch::dofs, "DOF numbers per node." "See :cpp:func:`GooseFEM::Mesh::Stitch::dofs`.") .def( "nodemap", py::overload_cast<>(&GooseFEM::Mesh::Stitch::nodemap, py::const_), "Node-map for givall sub-meshes." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodemap`.") .def( "elemmap", py::overload_cast<>(&GooseFEM::Mesh::Stitch::elemmap, py::const_), "Element-map for all sub-meshes." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemmap`.") .def( "nodemap", py::overload_cast(&GooseFEM::Mesh::Stitch::nodemap, py::const_), "Node-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodemap`.", py::arg("mesh_index")) .def( "elemmap", py::overload_cast(&GooseFEM::Mesh::Stitch::elemmap, py::const_), "Element-map for given sub-mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemmap`.", py::arg("mesh_index")) .def( "nodeset", static_cast (GooseFEM::Mesh::Stitch::*)( const xt::pytensor&, size_t) const>(&GooseFEM::Mesh::Stitch::nodeset), "Convert node-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.", py::arg("set"), py::arg("mesh_index")) .def( "elemset", static_cast (GooseFEM::Mesh::Stitch::*)( const xt::pytensor&, size_t) const>(&GooseFEM::Mesh::Stitch::elemset), "Convert element-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.", py::arg("set"), py::arg("mesh_index")) .def( "nodeset", static_cast (GooseFEM::Mesh::Stitch::*)( const std::vector>&) const>( &GooseFEM::Mesh::Stitch::nodeset), "Convert node-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.", py::arg("set")) .def( "elemset", static_cast (GooseFEM::Mesh::Stitch::*)( const std::vector>&) const>( &GooseFEM::Mesh::Stitch::elemset), "Convert element-set to the stitched mesh." "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.", py::arg("set")) .def("__repr__", [](const GooseFEM::Mesh::Stitch&) { return ""; }); py::class_(mod, "Vstack") .def( py::init(), "Vstack meshes." "See :cpp:class:`GooseFEM::Mesh::Vstack`.", py::arg("check_overlap") = true, py::arg("rtol") = 1e-5, py::arg("atol") = 1e-8) .def( "push_back", &GooseFEM::Mesh::Vstack::push_back< xt::pytensor, xt::pytensor, xt::pytensor>, "Add mesh to be stitched." "See :cpp:func:`GooseFEM::Mesh::Vstack::push_back`.", py::arg("coor"), py::arg("conn"), py::arg("nodes_bottom"), py::arg("nodes_top")) .def("__repr__", [](const GooseFEM::Mesh::Vstack&) { return ""; }); py::class_(mod, "Renumber") .def( py::init&>(), "Renumber to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Renumber`.", py::arg("dofs")) .def( "get", &GooseFEM::Mesh::Renumber::apply>, "Get renumbered DOFs." "See :cpp:func:`GooseFEM::Mesh::Renumber::get`.") .def( "apply", &GooseFEM::Mesh::Renumber::apply>, "Get renumbered list." "See :cpp:func:`GooseFEM::Mesh::Renumber::apply`.") .def( "index", &GooseFEM::Mesh::Renumber::index, "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``." "See :cpp:func:`GooseFEM::Mesh::Renumber::index`.") .def( "__repr__", [](const GooseFEM::Mesh::Renumber&) { return ""; }); py::class_(mod, "Reorder") .def( py::init([](xt::pytensor& a) { return new GooseFEM::Mesh::Reorder({a}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def( py::init([](xt::pytensor& a, xt::pytensor& b) { return new GooseFEM::Mesh::Reorder({a, b}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def( py::init([](xt::pytensor& a, xt::pytensor& b, xt::pytensor& c) { return new GooseFEM::Mesh::Reorder({a, b, c}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def( py::init([](xt::pytensor& a, xt::pytensor& b, xt::pytensor& c, xt::pytensor& d) { return new GooseFEM::Mesh::Reorder({a, b, c, d}); }), "Reorder to lowest possible index." "See :cpp:class:`GooseFEM::Mesh::Reorder`.") .def( "get", &GooseFEM::Mesh::Reorder::apply>, "Reorder matrix (e.g. ``dofs``)." "See :cpp:func:`GooseFEM::Mesh::Reorder::get`.", py::arg("dofs")) .def( "apply", &GooseFEM::Mesh::Reorder::apply>, "Get reordered list." "See :cpp:func:`GooseFEM::Mesh::Reorder::apply`.") .def( "index", &GooseFEM::Mesh::Reorder::index, "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``." "See :cpp:func:`GooseFEM::Mesh::Reorder::index`.") .def("__repr__", [](const GooseFEM::Mesh::Reorder&) { return ""; }); mod.def( "dofs", &GooseFEM::Mesh::dofs, "List with DOF-numbers in sequential order." "See :cpp:func:`GooseFEM::Mesh::dofs`.", py::arg("nnode"), py::arg("ndim")); mod.def( "renumber", &GooseFEM::Mesh::renumber>, "Renumber to lowest possible indices." "See :cpp:func:`GooseFEM::Mesh::renumber`.", py::arg("dofs")); mod.def( "coordination", &GooseFEM::Mesh::coordination>, "Coordination number of each node." "See :cpp:func:`GooseFEM::Mesh::coordination`.", py::arg("conn")); mod.def( "elem2node", &GooseFEM::Mesh::elem2node>, "Element-numbers connected to each node." "See :cpp:func:`GooseFEM::Mesh::elem2node`.", py::arg("conn"), py::arg("sorted") = true); mod.def( "edgesize", py::overload_cast&, const xt::pytensor&>( &GooseFEM::Mesh::edgesize, xt::pytensor>), "Get the edge size of all elements." "See :cpp:func:`GooseFEM::Mesh::edgesize`.", py::arg("coor"), py::arg("conn")); mod.def( "edgesize", py::overload_cast< const xt::pytensor&, const xt::pytensor&, GooseFEM::Mesh::ElementType>( &GooseFEM::Mesh::edgesize, xt::pytensor>), "Get the edge size of all elements." "See :cpp:func:`GooseFEM::Mesh::edgesize`.", py::arg("coor"), py::arg("conn"), py::arg("type")); mod.def( "centers", py::overload_cast&, const xt::pytensor&>( &GooseFEM::Mesh::centers, xt::pytensor>), "Coordinates of the center of each element." "See :cpp:func:`GooseFEM::Mesh::centers`.", py::arg("coor"), py::arg("conn")); mod.def( "centers", py::overload_cast< const xt::pytensor&, const xt::pytensor&, GooseFEM::Mesh::ElementType>( &GooseFEM::Mesh::centers, xt::pytensor>), "Coordinates of the center of each element." "See :cpp:func:`GooseFEM::Mesh::centers`.", py::arg("coor"), py::arg("conn"), py::arg("type")); mod.def( "elemmap2nodemap", py::overload_cast< const xt::pytensor&, const xt::pytensor&, const xt::pytensor&>(&GooseFEM::Mesh::elemmap2nodemap< xt::pytensor, xt::pytensor, xt::pytensor>), "Convert an element-map to a node-map." "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.", py::arg("elem_map"), py::arg("coor"), py::arg("conn")); mod.def( "elemmap2nodemap", py::overload_cast< const xt::pytensor&, const xt::pytensor&, const xt::pytensor&, GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::elemmap2nodemap< xt::pytensor, xt::pytensor, xt::pytensor>), "Convert an element-map to a node-map." "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.", py::arg("elem_map"), py::arg("coor"), py::arg("conn"), py::arg("type")); + + mod.def( + "center_of_gravity", + py::overload_cast< + const xt::pytensor&, + const xt::pytensor&>(&GooseFEM::Mesh::center_of_gravity< + xt::pytensor, + xt::pytensor>), + "Compute center of gravity of a mesh." + "See :cpp:func:`GooseFEM::Mesh::center_of_gravity`.", + py::arg("coor"), + py::arg("conn")); + + mod.def( + "center_of_gravity", + py::overload_cast< + const xt::pytensor&, + const xt::pytensor&, + GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::center_of_gravity< + xt::pytensor, + xt::pytensor>), + "Compute center of gravity of a mesh." + "See :cpp:func:`GooseFEM::Mesh::center_of_gravity`.", + py::arg("coor"), + py::arg("conn"), + py::arg("type")); } #endif diff --git a/tests/basic-python/Mesh.py b/tests/basic-python/Mesh.py new file mode 100644 index 0000000..ea955ab --- /dev/null +++ b/tests/basic-python/Mesh.py @@ -0,0 +1,26 @@ +import unittest + +import GooseFEM +import numpy as np + + +class Test_Mesh(unittest.TestCase): + """ + Simple mesh operations. + """ + + def test_center_of_gravity(self): + """ + Get the center of gravity of a mesh. + """ + + mesh = GooseFEM.Mesh.Quad4.Regular(3, 3) + coor = mesh.coor() + conn = mesh.conn() + + self.assertTrue(np.allclose([1.5, 1.5], GooseFEM.Mesh.center_of_gravity(coor, conn))) + + +if __name__ == "__main__": + + unittest.main() diff --git a/tests/basic-python/MeshQuad4.py b/tests/basic-python/MeshQuad4.py index 6705c6d..f63db2f 100644 --- a/tests/basic-python/MeshQuad4.py +++ b/tests/basic-python/MeshQuad4.py @@ -1,18 +1,25 @@ import unittest import GooseFEM as gf import numpy as np class Test_MeshQuad4(unittest.TestCase): + """ + Generation/manipulation of Quad4 meshes. + """ + def test_FineLayer_replica(self): + """ + Reconstruct Mesh.Quad4.FineLayer object from existing mesh. + """ mesh = gf.Mesh.Quad4.FineLayer(27, 27) replica = gf.Mesh.Quad4.FineLayer(mesh.coor(), mesh.conn()) self.assertTrue(np.allclose(mesh.coor(), replica.coor())) self.assertTrue(np.all(np.equal(mesh.conn(), replica.conn()))) if __name__ == "__main__": unittest.main()