diff --git a/README.md b/README.md index e1c4fbe..cd2533b 100644 --- a/README.md +++ b/README.md @@ -1,23 +1,23 @@ [![CI](https://github.com/tdegeus/GooseFEM/workflows/CI/badge.svg)](https://github.com/tdegeus/GooseFEM/actions) -[![Doxygen -> gh-pages](https://github.com/tdegeus/FrictionQPotFEM/workflows/gh-pages/badge.svg)](https://tdegeus.github.io/FrictionQPotFEM) +[![Doxygen -> gh-pages](https://github.com/tdegeus/GooseFEM/workflows/gh-pages/badge.svg)](https://tdegeus.github.io/GooseFEM) [![readthedocs](https://readthedocs.org/projects/goosefem/badge/?version=latest)](https://readthedocs.org/projects/goosefem/badge/?version=latest) ## GooseFEM Library to perform static or dynamic finite elements computations. The core of the implementation is a C++ library. For user convenience a Python interface is provided too. Please consult the documentation for more information: https://goosefem.readthedocs.io ## Credit / copyright (c) T.W.J. de Geus | [www.geus.me](http://www.geus.me) | [tom@geus.me](mailto:tom@geus.me) Tom de Geus was financially supported by: * [Swiss National Science Foundation (FNSF), Switzerland](http://www.snfs.ch) * [École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland](http://www.epfl.ch) * [Eindhoven University of Technology (TU/e), Eindhoven, The Netherlands](http://www.tue.nl) * [The Netherlands Research Council (NWO), The Netherlands](http://www.nwo.nl) * [Materials Innovation Institute (M2i), The Netherlands](http://www.m2i.nl) diff --git a/include/GooseFEM/ElementQuad4.h b/include/GooseFEM/ElementQuad4.h index ef9928b..0464cbd 100644 --- a/include/GooseFEM/ElementQuad4.h +++ b/include/GooseFEM/ElementQuad4.h @@ -1,139 +1,166 @@ /** Quadrature for 4-noded quadrilateral element (GooseFEM::Mesh::ElementType::Quad4). \file ElementQuad4.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENTQUAD4_H #define GOOSEFEM_ELEMENTQUAD4_H #include "config.h" #include "Element.h" namespace GooseFEM { namespace Element { namespace Quad4 { template inline double inv(const T& A, T& Ainv); namespace Gauss { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // namespace Gauss +/** +nodal quadrature: quadrature points coincide with the nodes. +The order is the same as in the connectivity:: + + 3 -- 2 + | | + 0 -- 1 +*/ namespace Nodal { -inline size_t nip(); // number of integration points -inline xt::xtensor xi(); // integration point coordinates (local coordinates) -inline xt::xtensor w(); // integration point weights + + /** + Number of integration points. + + \return unsigned int + */ + inline size_t nip(); // number of integration points + + /** + Integration point coordinates (local coordinates). + + \return Coordinates ``[nne, ndim]``. + */ + inline xt::xtensor xi(); + + /** + Integration point weights. + + \return Coordinates ``[nne]``. + */ + inline xt::xtensor w(); + } // namespace Nodal namespace MidPoint { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // namespace MidPoint class Quadrature : public GooseFEM::Element::QuadratureBase<4, 2, 2> { public: // Fixed dimensions: // ndim = 2 - number of dimensions // nne = 4 - number of nodes per element // // Naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - [nelem, nip] // Constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature() = default; Quadrature(const xt::xtensor& x); Quadrature( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w); // Update the nodal positions (shape of "x" should match the earlier definition) void update_x(const xt::xtensor& x); // Return shape function gradients xt::xtensor GradN() const; // Return integration volume xt::xtensor dV() const; /** Interpolate element vector. \param elemvec Nodal vector stored per element (shape: ``[nelem, nne, ndim]``). \param qvector Output: integration point vector (shape: ``[nelem, nip, ndim]``). */ template void interp_N_vector(const xt::xtensor& elemvec, xt::xtensor& qvector) const; /** Same as interp_N_vector(), but returns auto-allocated data. \param elemvec Nodal vector stored per element (shape: ``[nelem, nne, ndim]``). \return integration point vector (shape: ``[nelem, nip, ndim]``). */ template xt::xtensor Interp_N_vector(const xt::xtensor& elemvec) const; // Dyadic product (and its transpose and symmetric part) // qtensor(i,j) += dNdx(m,i) * elemvec(m,j) void gradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void gradN_vector_T(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void symGradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; // Integral of the scalar product // elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV void int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const; // Integral of the dot product // elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV void int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const; // Integral of the dot product // elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV void int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const; // Auto-allocation of the functions above xt::xtensor GradN_vector(const xt::xtensor& elemvec) const; xt::xtensor GradN_vector_T(const xt::xtensor& elemvec) const; xt::xtensor SymGradN_vector(const xt::xtensor& elemvec) const; xt::xtensor Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const; xt::xtensor Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const; xt::xtensor Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const; private: // Compute "vol" and "dNdx" based on current "x" void compute_dN(); private: xt::xtensor m_x; // nodal positions stored per element [nelem, nne, ndim] xt::xtensor m_w; // weight of each integration point [nip] xt::xtensor m_xi; // local coordinate of each integration point [nip, ndim] xt::xtensor m_N; // shape functions [nip, nne] xt::xtensor m_dNxi; // shape function grad. wrt local coor. [nip, nne, ndim] xt::xtensor m_dNx; // shape function grad. wrt global coor. [nelem, nip, nne, ndim] xt::xtensor m_vol; // integration point volume [nelem, nip] }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #include "ElementQuad4.hpp" #endif diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h index e5298af..3822777 100644 --- a/include/GooseFEM/MeshQuad4.h +++ b/include/GooseFEM/MeshQuad4.h @@ -1,325 +1,382 @@ /** Generate mesh with 4-noded quadrilateral elements. \file MeshQuad4.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MESHQUAD4_H #define GOOSEFEM_MESHQUAD4_H #include "config.h" namespace GooseFEM { namespace Mesh { namespace Quad4 { // pre-allocation namespace Map { class FineLayer2Regular; } -// Regular mesh: equi-sized elements - +/** +Regular mesh: equi-sized elements. +*/ class Regular { public: + Regular() = default; + + /** + Constructor. + + \param nelx Number of elements in horizontal (x) direction. + \param nely Number of elements in vertical (y) direction. + \param h Edge size (width == height). + */ Regular(size_t nelx, size_t nely, double h = 1.0); - // size - size_t nelem() const; // number of elements - size_t nnode() const; // number of nodes - size_t nne() const; // number of nodes-per-element - size_t ndim() const; // number of dimensions - size_t nelx() const; // number of elements in x-direction - size_t nely() const; // number of elements in y-direction - double h() const; // edge size + /** + Number of elements. + + \return unsigned int. + */ + size_t nelem() const; + + /** + Number of nodes. + + \return unsigned int. + */ + size_t nnode() const; + + /** + Number of nodes-per-element. + + \return unsigned int. + */ + size_t nne() const; + + /** + Number of dimensions. + + \return unsigned int. + */ + size_t ndim() const; + + /** + Number of elements in x-direction. + + \return unsigned int. + */ + size_t nelx() const; + + /** + Number of elements in y-direction. + + \return unsigned int. + */ + size_t nely() const; + + /** + Edge size. + + \return double. + */ + double h() const; /** Get the element type. \return GooseFEM::Mesh::ElementType(). */ ElementType getElementType() const; // mesh xt::xtensor coor() const; // nodal positions [nnode, ndim] xt::xtensor conn() const; // connectivity [nelem, nne] // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor nodesRightOpenEdge() const; // boundary nodes: corners (including aliases) size_t nodesBottomLeftCorner() const; size_t nodesBottomRightCorner() const; size_t nodesTopLeftCorner() const; size_t nodesTopRightCorner() const; size_t nodesLeftBottomCorner() const; size_t nodesLeftTopCorner() const; size_t nodesRightBottomCorner() const; size_t nodesRightTopCorner() const; // DOF-numbers for each component of each node (sequential) xt::xtensor dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; // element numbers as matrix xt::xtensor elementgrid() const; private: double m_h; // elementary element edge-size (in all directions) size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") size_t m_nelem; // number of elements size_t m_nnode; // number of nodes static const size_t m_nne = 4; // number of nodes-per-element static const size_t m_ndim = 2; // number of dimensions }; // Mesh with fine middle layer, and coarser elements towards the top and bottom class FineLayer { public: FineLayer() = default; FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1); // Reconstruct class for given coordinates / connectivity FineLayer(const xt::xtensor& coor, const xt::xtensor& conn); // size size_t nelem() const; // number of elements size_t nnode() const; // number of nodes size_t nne() const; // number of nodes-per-element size_t ndim() const; // number of dimensions size_t nelx() const; // number of elements in x-direction size_t nely() const; // number of elements in y-direction double h() const; // edge size // edge size, per row of elements (in units of "h") xt::xtensor elemrow_nhx() const; xt::xtensor elemrow_nhy() const; xt::xtensor elemrow_nelem() const; // type ElementType getElementType() const; // mesh xt::xtensor coor() const; // nodal positions [nnode, ndim] xt::xtensor conn() const; // connectivity [nelem, nne] // elements in the middle (fine) layer xt::xtensor elementsMiddleLayer() const; // extract elements along a layer xt::xtensor elementsLayer(size_t layer) const; // select region of elements from 'matrix' of element numbers xt::xtensor elementgrid_ravel( std::vector rows_start_stop, std::vector cols_start_stop) const; /** Select region of elements from 'matrix' of element numbers around an element: square box with edge-size ``(2 * size + 1) * h``, around ``element``. \param element The element around which to select elements. \param size Edge size of the square box encapsulating the selected element. \param periodic Assume the mesh periodic. \returns List of elements. */ xt::xtensor elementgrid_around_ravel( size_t element, size_t size, bool periodic = true); /** Select region of elements from 'matrix' of element numbers around an element: left/right from ``element`` (on the same layer). \param element The element around which to select elements. \param left Number of elements to select to the left. \param right Number of elements to select to the right. \param periodic Assume the mesh periodic. \returns List of elements. */ // - xt::xtensor elementgrid_leftright( size_t element, size_t left, size_t right, bool periodic = true); // boundary nodes: edges xt::xtensor nodesBottomEdge() const; xt::xtensor nodesTopEdge() const; xt::xtensor nodesLeftEdge() const; xt::xtensor nodesRightEdge() const; // boundary nodes: edges, without corners xt::xtensor nodesBottomOpenEdge() const; xt::xtensor nodesTopOpenEdge() const; xt::xtensor nodesLeftOpenEdge() const; xt::xtensor nodesRightOpenEdge() const; // boundary nodes: corners (including aliases) size_t nodesBottomLeftCorner() const; size_t nodesBottomRightCorner() const; size_t nodesTopLeftCorner() const; size_t nodesTopRightCorner() const; size_t nodesLeftBottomCorner() const; size_t nodesLeftTopCorner() const; size_t nodesRightBottomCorner() const; size_t nodesRightTopCorner() const; // DOF-numbers for each component of each node (sequential) xt::xtensor dofs() const; // DOF-numbers for the case that the periodicity if fully eliminated xt::xtensor dofsPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) xt::xtensor nodesPeriodic() const; // front-bottom-left node, used as reference for periodicity size_t nodesOrigin() const; // mapping to 'roll' periodically in the x-direction, // returns element mapping, such that: new_elemvar = elemvar[elem_map] xt::xtensor roll(size_t n); private: double m_h; // elementary element edge-size (in all directions) double m_Lx; // mesh size in "x" size_t m_nelem; // number of elements size_t m_nnode; // number of nodes static const size_t m_nne = 4; // number of nodes-per-element static const size_t m_ndim = 2; // number of dimensions xt::xtensor m_nelx; // number of elements in "x" (*) xt::xtensor m_nnd; // total number of nodes in the main node layer (**) xt::xtensor m_nhx; // element size in x-direction (*) xt::xtensor m_nhy; // element size in y-direction (*) xt::xtensor m_refine; // refine direction (-1:no refine, 0:"x" (*) xt::xtensor m_startElem; // start element (*) xt::xtensor m_startNode; // start node (**) // (*) per element layer in "y" // (**) per node layer in "y" void init(size_t nelx, size_t nely, double h, size_t nfine = 1); void map(const xt::xtensor& coor, const xt::xtensor& conn); friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular; }; // Mesh mappings namespace Map { // Return "FineLayer"-class responsible for generating a connectivity // Throws if conversion is not possible GooseFEM::Mesh::Quad4::FineLayer FineLayer( const xt::xtensor& coor, const xt::xtensor& conn); // Refine a regular mesh: sub-divide elements in several smaller elements class RefineRegular { public: // Constructors RefineRegular() = default; RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny); // return the coarse or the fine mesh objects GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const; GooseFEM::Mesh::Quad4::Regular getFineMesh() const; // elements of the fine mesh per element of the coarse mesh xt::xtensor getMap() const; // map field xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToCoarse(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToCoarse(const xt::xtensor& data) const; // tensor per intpnt // map field xt::xtensor mapToFine(const xt::xtensor& data) const; // scalar per el xt::xtensor mapToFine(const xt::xtensor& data) const; // scalar per intpnt xt::xtensor mapToFine(const xt::xtensor& data) const; // tensor per intpnt private: // the meshes GooseFEM::Mesh::Quad4::Regular m_coarse; GooseFEM::Mesh::Quad4::Regular m_fine; // mapping xt::xtensor m_fine2coarse; xt::xtensor m_fine2coarse_index; xt::xtensor m_coarse2fine; }; class FineLayer2Regular { public: // constructor FineLayer2Regular() = default; FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh); // return either of the meshes GooseFEM::Mesh::Quad4::Regular getRegularMesh() const; GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const; // elements of the Regular mesh per element of the FineLayer mesh // and the fraction by which the overlap is std::vector> getMap() const; std::vector> getMapFraction() const; /** - Map integration point quantities to Regular. + Map element quantities to Regular. + The mapping is a bit simplistic: no interpolation is involved, the function just + accounts the fraction of overlap between the FineLayer element and the Regular element. + The mapping is such that:: + + ret[e, :, :] <- arg[e, :, :] + + (for arbitrary rank). \tparam T The type of the data (e.g. ``double``). \tparam rank Rank of the data. \param arg The data. \return The mapped data. */ template xt::xtensor mapToRegular(const xt::xtensor& arg) const; private: // the "FineLayer" mesh to map GooseFEM::Mesh::Quad4::FineLayer m_finelayer; // the new "Regular" mesh to which to map GooseFEM::Mesh::Quad4::Regular m_regular; // mapping std::vector> m_elem_regular; std::vector> m_frac_regular; }; } // namespace Map } // namespace Quad4 } // namespace Mesh } // namespace GooseFEM #include "MeshQuad4.hpp" #endif diff --git a/include/GooseFEM/MeshQuad4.hpp b/include/GooseFEM/MeshQuad4.hpp index ac7fc7a..5f44ee8 100644 --- a/include/GooseFEM/MeshQuad4.hpp +++ b/include/GooseFEM/MeshQuad4.hpp @@ -1,1618 +1,1617 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_MESHQUAD4_HPP #define GOOSEFEM_MESHQUAD4_HPP #include "MeshQuad4.h" namespace GooseFEM { namespace Mesh { namespace Quad4 { inline Regular::Regular(size_t nelx, size_t nely, double h) : m_h(h), m_nelx(nelx), m_nely(nely) { GOOSEFEM_ASSERT(m_nelx >= 1ul); GOOSEFEM_ASSERT(m_nely >= 1ul); m_nnode = (m_nelx + 1) * (m_nely + 1); m_nelem = m_nelx * m_nely; } inline size_t Regular::nelem() const { return m_nelem; } inline size_t Regular::nnode() const { return m_nnode; } inline size_t Regular::nne() const { return m_nne; } inline size_t Regular::ndim() const { return m_ndim; } inline size_t Regular::nelx() const { return m_nelx; } inline size_t Regular::nely() const { return m_nely; } inline double Regular::h() const { return m_h; } inline ElementType Regular::getElementType() const { return ElementType::Quad4; } inline xt::xtensor Regular::coor() const { xt::xtensor ret = xt::empty({m_nnode, m_ndim}); xt::xtensor x = xt::linspace(0.0, m_h * static_cast(m_nelx), m_nelx + 1); xt::xtensor y = xt::linspace(0.0, m_h * static_cast(m_nely), m_nely + 1); size_t inode = 0; for (size_t iy = 0; iy < m_nely + 1; ++iy) { for (size_t ix = 0; ix < m_nelx + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ++inode; } } return ret; } inline xt::xtensor Regular::conn() const { xt::xtensor ret = xt::empty({m_nelem, m_nne}); size_t ielem = 0; for (size_t iy = 0; iy < m_nely; ++iy) { for (size_t ix = 0; ix < m_nelx; ++ix) { ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix); ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1); ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix); ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1); ++ielem; } } return ret; } inline xt::xtensor Regular::nodesBottomEdge() const { return xt::arange(m_nelx + 1); } inline xt::xtensor Regular::nodesTopEdge() const { return xt::arange(m_nelx + 1) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightEdge() const { return xt::arange(m_nely + 1) * (m_nelx + 1) + m_nelx; } inline xt::xtensor Regular::nodesBottomOpenEdge() const { return xt::arange(1, m_nelx); } inline xt::xtensor Regular::nodesTopOpenEdge() const { return xt::arange(1, m_nelx) + m_nely * (m_nelx + 1); } inline xt::xtensor Regular::nodesLeftOpenEdge() const { return xt::arange(1, m_nely) * (m_nelx + 1); } inline xt::xtensor Regular::nodesRightOpenEdge() const { return xt::arange(1, m_nely) * (m_nelx + 1) + m_nelx; } inline size_t Regular::nodesBottomLeftCorner() const { return 0; } inline size_t Regular::nodesBottomRightCorner() const { return m_nelx; } inline size_t Regular::nodesTopLeftCorner() const { return m_nely * (m_nelx + 1); } inline size_t Regular::nodesTopRightCorner() const { return m_nely * (m_nelx + 1) + m_nelx; } inline size_t Regular::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } inline size_t Regular::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } inline size_t Regular::nodesRightBottomCorner() const { return nodesBottomRightCorner(); } inline size_t Regular::nodesRightTopCorner() const { return nodesTopRightCorner(); } inline xt::xtensor Regular::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(shape); ret(0, 0) = nodesBottomLeftCorner(); ret(0, 1) = nodesBottomRightCorner(); ret(1, 0) = nodesBottomLeftCorner(); ret(1, 1) = nodesTopRightCorner(); ret(2, 0) = nodesBottomLeftCorner(); ret(2, 1) = nodesTopLeftCorner(); size_t i = 3; xt::view(ret, xt::range(i, i + bot.size()), 0) = bot; xt::view(ret, xt::range(i, i + bot.size()), 1) = top; i += bot.size(); xt::view(ret, xt::range(i, i + lft.size()), 0) = lft; xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt; return ret; } inline size_t Regular::nodesOrigin() const { return nodesBottomLeftCorner(); } inline xt::xtensor Regular::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor Regular::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor dependent = xt::view(nodePer, xt::all(), 1); for (size_t j = 0; j < m_ndim; ++j) { xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j); } return GooseFEM::Mesh::renumber(ret); } inline xt::xtensor Regular::elementgrid() const { return xt::arange(m_nelem).reshape({m_nely, m_nelx}); } inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine) { this->init(nelx, nely, h, nfine); } inline FineLayer::FineLayer(const xt::xtensor& coor, const xt::xtensor& conn) { this->map(coor, conn); } inline void FineLayer::init(size_t nelx, size_t nely, double h, size_t nfine) { GOOSEFEM_ASSERT(nelx >= 1ul); GOOSEFEM_ASSERT(nely >= 1ul); m_h = h; m_Lx = m_h * static_cast(nelx); // compute element size in y-direction (use symmetry, compute upper half) // temporary variables size_t nmin, ntot; xt::xtensor nhx = xt::ones({nely}); xt::xtensor nhy = xt::ones({nely}); xt::xtensor refine = -1 * xt::ones({nely}); // minimum height in y-direction (half of the height because of symmetry) if (nely % 2 == 0) { nmin = nely / 2; } else { nmin = (nely + 1) / 2; } // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half) if (nfine % 2 == 0) { nfine = nfine / 2 + 1; } else { nfine = (nfine + 1) / 2; } if (nfine < 1) { nfine = 1; } if (nfine > nmin) { nfine = nmin; } // loop over element layers in y-direction, try to coarsen using these rules: // (1) element size in y-direction <= distance to origin in y-direction // (2) element size in x-direction should fit the total number of elements in x-direction // (3) a certain number of layers have the minimum size "1" (are fine) for (size_t iy = nfine;;) { // initialize current size in y-direction if (iy == nfine) { ntot = nfine; } // check to stop if (iy >= nely || ntot >= nmin) { nely = iy; break; } // rules (1,2) satisfied: coarsen in x-direction if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) { refine(iy) = 0; nhy(iy) *= 2; auto vnhy = xt::view(nhy, xt::range(iy + 1, _)); auto vnhx = xt::view(nhx, xt::range(iy, _)); vnhy *= 3; vnhx *= 3; } // update the number of elements in y-direction ntot += nhy(iy); // proceed to next element layer in y-direction ++iy; // check to stop if (iy >= nely || ntot >= nmin) { nely = iy; break; } } // symmetrize, compute full information // allocate mesh constructor parameters m_nhx = xt::empty({nely * 2 - 1}); m_nhy = xt::empty({nely * 2 - 1}); m_refine = xt::empty({nely * 2 - 1}); m_nelx = xt::empty({nely * 2 - 1}); m_nnd = xt::empty({nely * 2}); m_startElem = xt::empty({nely * 2 - 1}); m_startNode = xt::empty({nely * 2}); // fill // - lower half for (size_t iy = 0; iy < nely; ++iy) { m_nhx(iy) = nhx(nely - iy - 1); m_nhy(iy) = nhy(nely - iy - 1); m_refine(iy) = refine(nely - iy - 1); } // - upper half for (size_t iy = 0; iy < nely - 1; ++iy) { m_nhx(iy + nely) = nhx(iy + 1); m_nhy(iy + nely) = nhy(iy + 1); m_refine(iy + nely) = refine(iy + 1); } // update size nely = m_nhx.size(); // compute the number of elements per element layer in y-direction for (size_t iy = 0; iy < nely; ++iy) { m_nelx(iy) = nelx / m_nhx(iy); } // compute the number of nodes per node layer in y-direction for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) { m_nnd(iy) = m_nelx(iy) + 1; } for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { m_nnd(iy + 1) = m_nelx(iy) + 1; } // compute mesh dimensions // initialize m_nnode = 0; m_nelem = 0; m_startNode(0) = 0; // loop over element layers (bottom -> middle, elements become finer) for (size_t i = 0; i < (nely - 1) / 2; ++i) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if (m_refine(i) == 0) { m_nnode += (3 * m_nelx(i) + 1); } else { m_nnode += (m_nelx(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_nelx(i)); } else { m_nelem += (m_nelx(i)); } // - store the starting node of the next layer m_startNode(i + 1) = m_nnode; } // loop over element layers (middle -> top, elements become coarser) for (size_t i = (nely - 1) / 2; i < nely; ++i) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if (m_refine(i) == 0) { m_nnode += (5 * m_nelx(i) + 1); } else { m_nnode += (m_nelx(i) + 1); } // - add the elements of this layer if (m_refine(i) == 0) { m_nelem += (4 * m_nelx(i)); } else { m_nelem += (m_nelx(i)); } // - store the starting node of the next layer m_startNode(i + 1) = m_nnode; } // - add the top row of nodes m_nnode += m_nelx(nely - 1) + 1; } inline size_t FineLayer::nelem() const { return m_nelem; } inline size_t FineLayer::nnode() const { return m_nnode; } inline size_t FineLayer::nne() const { return m_nne; } inline size_t FineLayer::ndim() const { return m_ndim; } inline size_t FineLayer::nelx() const { return xt::amax(m_nelx)(); } inline size_t FineLayer::nely() const { return xt::sum(m_nhy)(); } inline double FineLayer::h() const { return m_h; } inline xt::xtensor FineLayer::elemrow_nhx() const { return m_nhx; } inline xt::xtensor FineLayer::elemrow_nhy() const { return m_nhy; } inline xt::xtensor FineLayer::elemrow_nelem() const { return m_nelx; } inline ElementType FineLayer::getElementType() const { return ElementType::Quad4; } inline xt::xtensor FineLayer::coor() const { // allocate output xt::xtensor ret = xt::empty({m_nnode, m_ndim}); // current node, number of element layers size_t inode = 0; size_t nely = static_cast(m_nhy.size()); // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) // - allocate xt::xtensor y = xt::empty({nely + 1}); // - initialize y(0) = 0.0; // - compute for (size_t iy = 1; iy < nely + 1; ++iy) { y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h; } // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes for (size_t iy = 0;; ++iy) { // get positions along the x- and z-axis xt::xtensor x = xt::linspace(0.0, m_Lx, m_nelx(iy) + 1); // add nodes of the bottom layer of this element for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy); ++inode; } // stop at middle layer if (iy == (nely - 1) / 2) { break; } // add extra nodes of the intermediate layer, for refinement in x-direction if (m_refine(iy) == 0) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t ix = 0; ix < m_nelx(iy); ++ix) { for (size_t j = 0; j < 2; ++j) { ret(inode, 0) = x(ix) + dx * static_cast(j + 1); ret(inode, 1) = y(iy) + dy; ++inode; } } } } // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) { // get positions along the x- and z-axis xt::xtensor x = xt::linspace(0.0, m_Lx, m_nelx(iy) + 1); // add extra nodes of the intermediate layer, for refinement in x-direction if (m_refine(iy) == 0) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy) / 3); double dy = m_h * static_cast(m_nhy(iy) / 2); // - add nodes of the intermediate layer for (size_t ix = 0; ix < m_nelx(iy); ++ix) { for (size_t j = 0; j < 2; ++j) { ret(inode, 0) = x(ix) + dx * static_cast(j + 1); ret(inode, 1) = y(iy) + dy; ++inode; } } } // add nodes of the top layer of this element for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) { ret(inode, 0) = x(ix); ret(inode, 1) = y(iy + 1); ++inode; } } return ret; } inline xt::xtensor FineLayer::conn() const { // allocate output xt::xtensor ret = xt::empty({m_nelem, m_nne}); // current element, number of element layers, starting nodes of each node layer size_t ielem = 0; size_t nely = static_cast(m_nhy.size()); size_t bot, mid, top; // loop over all element layers for (size_t iy = 0; iy < nely; ++iy) { // - get: starting nodes of bottom(, middle) and top layer bot = m_startNode(iy); mid = m_startNode(iy) + m_nnd(iy); top = m_startNode(iy + 1); // - define connectivity: no coarsening/refinement if (m_refine(iy) == -1) { for (size_t ix = 0; ix < m_nelx(iy); ++ix) { ret(ielem, 0) = bot + (ix); ret(ielem, 1) = bot + (ix + 1); ret(ielem, 2) = top + (ix + 1); ret(ielem, 3) = top + (ix); ielem++; } } // - define connectivity: refinement along the x-direction (below the middle layer) else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) { for (size_t ix = 0; ix < m_nelx(iy); ++ix) { // -- bottom element ret(ielem, 0) = bot + (ix); ret(ielem, 1) = bot + (ix + 1); ret(ielem, 2) = mid + (2 * ix + 1); ret(ielem, 3) = mid + (2 * ix); ielem++; // -- top-right element ret(ielem, 0) = bot + (ix + 1); ret(ielem, 1) = top + (3 * ix + 3); ret(ielem, 2) = top + (3 * ix + 2); ret(ielem, 3) = mid + (2 * ix + 1); ielem++; // -- top-center element ret(ielem, 0) = mid + (2 * ix); ret(ielem, 1) = mid + (2 * ix + 1); ret(ielem, 2) = top + (3 * ix + 2); ret(ielem, 3) = top + (3 * ix + 1); ielem++; // -- top-left element ret(ielem, 0) = bot + (ix); ret(ielem, 1) = mid + (2 * ix); ret(ielem, 2) = top + (3 * ix + 1); ret(ielem, 3) = top + (3 * ix); ielem++; } } // - define connectivity: coarsening along the x-direction (above the middle layer) else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) { for (size_t ix = 0; ix < m_nelx(iy); ++ix) { // -- lower-left element ret(ielem, 0) = bot + (3 * ix); ret(ielem, 1) = bot + (3 * ix + 1); ret(ielem, 2) = mid + (2 * ix); ret(ielem, 3) = top + (ix); ielem++; // -- lower-center element ret(ielem, 0) = bot + (3 * ix + 1); ret(ielem, 1) = bot + (3 * ix + 2); ret(ielem, 2) = mid + (2 * ix + 1); ret(ielem, 3) = mid + (2 * ix); ielem++; // -- lower-right element ret(ielem, 0) = bot + (3 * ix + 2); ret(ielem, 1) = bot + (3 * ix + 3); ret(ielem, 2) = top + (ix + 1); ret(ielem, 3) = mid + (2 * ix + 1); ielem++; // -- upper element ret(ielem, 0) = mid + (2 * ix); ret(ielem, 1) = mid + (2 * ix + 1); ret(ielem, 2) = top + (ix + 1); ret(ielem, 3) = top + (ix); ielem++; } } } return ret; } inline xt::xtensor FineLayer::elementsMiddleLayer() const { size_t nely = m_nhy.size(); size_t iy = (nely - 1) / 2; return m_startElem(iy) + xt::arange(m_nelx(iy)); } inline xt::xtensor FineLayer::elementsLayer(size_t iy) const { GOOSEFEM_ASSERT(iy < m_nelx.size()); size_t n = m_nelx(iy); if (m_refine(iy) != -1) { n *= 4; } return m_startElem(iy) + xt::arange(n); } inline xt::xtensor FineLayer::elementgrid_ravel( std::vector start_stop_rows, std::vector start_stop_cols) const { GOOSEFEM_ASSERT(start_stop_rows.size() == 0 || start_stop_rows.size() == 2); GOOSEFEM_ASSERT(start_stop_cols.size() == 0 || start_stop_cols.size() == 2); std::array rows; std::array cols; if (start_stop_rows.size() == 2) { std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin()); GOOSEFEM_ASSERT(rows[0] <= this->nely()); GOOSEFEM_ASSERT(rows[1] <= this->nely()); } else { rows[0] = 0; rows[1] = this->nely(); } if (start_stop_cols.size() == 2) { std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin()); GOOSEFEM_ASSERT(cols[0] <= this->nelx()); GOOSEFEM_ASSERT(cols[1] <= this->nelx()); } else { cols[0] = 0; cols[1] = this->nelx(); } if (rows[0] == rows[1] || cols[0] == cols[1]) { xt::xtensor ret = xt::empty({0}); return ret; } // Compute dimensions auto H = xt::cumsum(m_nhy); size_t yl = 0; if (rows[0] > 0) { yl = xt::argmax(H > rows[0])(); } size_t yu = xt::argmax(H >= rows[1])(); size_t hx = std::max(m_nhx(yl), m_nhx(yu)); size_t xl = (cols[0] - cols[0] % hx) / hx; size_t xu = (cols[1] - cols[1] % hx) / hx; // Allocate output size_t N = 0; for (size_t i = yl; i <= yu; ++i) { // no refinement size_t n = (xu - xl) * hx / m_nhx(i); // refinement if (m_refine(i) != -1) { n *= 4; } N += n; } xt::xtensor ret = xt::empty({N}); // Write output N = 0; for (size_t i = yl; i <= yu; ++i) { // no refinement size_t n = (xu - xl) * hx / m_nhx(i); size_t h = hx; // refinement if (m_refine(i) != -1) { n *= 4; h *= 4; } xt::xtensor e = m_startElem(i) + xl * h / m_nhx(i) + xt::arange(n); xt::view(ret, xt::range(N, N + n)) = e; N += n; } return ret; } inline xt::xtensor FineLayer::elementgrid_around_ravel( size_t e, size_t size, bool periodic) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_nelx(iy); GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2); auto H = xt::cumsum(m_nhy); if (2 * size >= H(H.size() - 1)) { return xt::arange(this->nelem()); } size_t hy = H(iy); size_t l = xt::argmax(H > (hy - size - 1))(); size_t u = xt::argmax(H >= (hy + size))(); size_t lh = 0; if (l > 0) { lh = H(l - 1); } size_t uh = H(u); size_t step = xt::amax(m_nhx)(); size_t relx = (e - m_startElem(iy)) % step; size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx; size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step; size_t dx = m_nhx(u); size_t xl = 0; size_t xu = nel; if (mid >= size) { xl = mid - size; } if (mid + size < nel) { xu = mid + size + 1; } xl = xl - xl % dx; xu = xu - xu % dx; if (mid - xl < size) { if (xl < dx) { xl = 0; } else { xl -= dx; } } if (xu - mid < size) { if (xu > nel - dx) { xu = nel; } else { xu += dx; } } auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu}); auto map = this->roll(nroll); return xt::view(map, xt::keep(ret)); } inline xt::xtensor FineLayer::elementgrid_leftright( size_t e, size_t left, size_t right, bool periodic) { GOOSEFEM_WIP_ASSERT(periodic == true); size_t iy = xt::argmin(m_startElem <= e)() - 1; size_t nel = m_nelx(iy); GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2); size_t step = xt::amax(m_nhx)(); size_t relx = (e - m_startElem(iy)) % step; size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx; size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step; size_t dx = m_nhx(iy); size_t xl = 0; size_t xu = nel; if (mid >= left) { xl = mid - left; } if (mid + right < nel) { xu = mid + right + 1; } xl = xl - xl % dx; xu = xu - xu % dx; if (mid - xl < left) { if (xl < dx) { xl = 0; } else { xl -= dx; } } if (xu - mid < right) { if (xu > nel - dx) { xu = nel; } else { xu += dx; } } auto H = xt::cumsum(m_nhy); size_t yl = 0; if (iy > 0) { yl = H(iy - 1); } auto ret = this->elementgrid_ravel({yl, H(iy)}, {xl, xu}); auto map = this->roll(nroll); return xt::view(map, xt::keep(ret)); } inline xt::xtensor FineLayer::nodesBottomEdge() const { return m_startNode(0) + xt::arange(m_nelx(0) + 1); } inline xt::xtensor FineLayer::nodesTopEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(m_nelx(nely - 1) + 1); } inline xt::xtensor FineLayer::nodesLeftEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely + 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j)); xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1)); return ret; } inline xt::xtensor FineLayer::nodesRightEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely + 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j)) + xt::view(m_nelx, xt::range(i, j)); xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1)) + xt::view(m_nelx, xt::range(k, l)); return ret; } inline xt::xtensor FineLayer::nodesBottomOpenEdge() const { return m_startNode(0) + xt::arange(1, m_nelx(0)); } inline xt::xtensor FineLayer::nodesTopOpenEdge() const { size_t nely = m_nhy.size(); return m_startNode(nely) + xt::arange(1, m_nelx(nely - 1)); } inline xt::xtensor FineLayer::nodesLeftOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely - 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j)); xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l)); return ret; } inline xt::xtensor FineLayer::nodesRightOpenEdge() const { size_t nely = m_nhy.size(); xt::xtensor ret = xt::empty({nely - 1}); size_t i = 0; size_t j = (nely + 1) / 2; size_t k = (nely - 1) / 2; size_t l = nely; xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j)) + xt::view(m_nelx, xt::range(i + 1, j)); xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l)) + xt::view(m_nelx, xt::range(k, l - 1)); return ret; } inline size_t FineLayer::nodesBottomLeftCorner() const { return m_startNode(0); } inline size_t FineLayer::nodesBottomRightCorner() const { return m_startNode(0) + m_nelx(0); } inline size_t FineLayer::nodesTopLeftCorner() const { size_t nely = m_nhy.size(); return m_startNode(nely); } inline size_t FineLayer::nodesTopRightCorner() const { size_t nely = m_nhy.size(); return m_startNode(nely) + m_nelx(nely - 1); } inline size_t FineLayer::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } inline size_t FineLayer::nodesRightBottomCorner() const { return nodesBottomRightCorner(); } inline size_t FineLayer::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } inline size_t FineLayer::nodesRightTopCorner() const { return nodesTopRightCorner(); } inline xt::xtensor FineLayer::nodesPeriodic() const { xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); std::array shape = {bot.size() + lft.size() + 3ul, 2ul}; xt::xtensor ret = xt::empty(shape); ret(0, 0) = nodesBottomLeftCorner(); ret(0, 1) = nodesBottomRightCorner(); ret(1, 0) = nodesBottomLeftCorner(); ret(1, 1) = nodesTopRightCorner(); ret(2, 0) = nodesBottomLeftCorner(); ret(2, 1) = nodesTopLeftCorner(); size_t i = 3; xt::view(ret, xt::range(i, i + bot.size()), 0) = bot; xt::view(ret, xt::range(i, i + bot.size()), 1) = top; i += bot.size(); xt::view(ret, xt::range(i, i + lft.size()), 0) = lft; xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt; return ret; } inline size_t FineLayer::nodesOrigin() const { return nodesBottomLeftCorner(); } inline xt::xtensor FineLayer::dofs() const { return GooseFEM::Mesh::dofs(m_nnode, m_ndim); } inline xt::xtensor FineLayer::dofsPeriodic() const { xt::xtensor ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim); xt::xtensor nodePer = nodesPeriodic(); xt::xtensor independent = xt::view(nodePer, xt::all(), 0); xt::xtensor dependent = xt::view(nodePer, xt::all(), 1); for (size_t j = 0; j < m_ndim; ++j) { xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j); } return GooseFEM::Mesh::renumber(ret); } inline xt::xtensor FineLayer::roll(size_t n) { auto conn = this->conn(); size_t nely = static_cast(m_nhy.size()); xt::xtensor ret = xt::empty({m_nelem}); // loop over all element layers for (size_t iy = 0; iy < nely; ++iy) { // no refinement size_t shift = n * (m_nelx(iy) / m_nelx(0)); size_t nel = m_nelx(iy); // refinement if (m_refine(iy) != -1) { shift = n * (m_nelx(iy) / m_nelx(0)) * 4; nel = m_nelx(iy) * 4; } // element numbers of the layer, and roll them auto e = m_startElem(iy) + xt::arange(nel); xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift); } return ret; } inline void FineLayer::map(const xt::xtensor& coor, const xt::xtensor& conn) { GOOSEFEM_ASSERT(coor.shape(1) == 2); GOOSEFEM_ASSERT(conn.shape(1) == 4); GOOSEFEM_ASSERT(conn.shape(0) > 0); GOOSEFEM_ASSERT(coor.shape(0) >= 4); GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0)); if (conn.shape(0) == 1) { this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0)); GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn))); GOOSEFEM_CHECK(xt::allclose(this->coor(), coor)); return; } // Identify the middle layer size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2; size_t eleft = emid; size_t eright = emid; while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) { eleft--; } while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) { eright++; } GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0)); // Get element sizes along the middle layer auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0); auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1); auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 2); auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0); auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1); auto hx = xt::amin(dx)(); auto hy = xt::amin(dy)(); GOOSEFEM_CHECK(xt::allclose(hx, hy)); GOOSEFEM_CHECK(xt::allclose(dx, hx)); GOOSEFEM_CHECK(xt::allclose(dy, hy)); // Extract shape and initialise size_t nelx = eright - eleft + 1; size_t nely = static_cast((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx); this->init(nelx, nely, hx); GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn))); GOOSEFEM_CHECK(xt::allclose(this->coor(), coor)); GOOSEFEM_CHECK(xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange(nelx)))); } namespace Map { inline RefineRegular::RefineRegular( const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny) : m_coarse(mesh) { m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h()); xt::xtensor elmat_coarse = m_coarse.elementgrid(); xt::xtensor elmat_fine = m_fine.elementgrid(); m_coarse2fine = xt::empty({m_coarse.nelem(), nx * ny}); for (size_t i = 0; i < elmat_coarse.shape(0); ++i) { for (size_t j = 0; j < elmat_coarse.shape(1); ++j) { xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view( elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx))); } } } inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getCoarseMesh() const { return m_coarse; } inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getFineMesh() const { return m_fine; } inline xt::xtensor RefineRegular::getMap() const { return m_coarse2fine; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); size_t m = m_coarse2fine.shape(0); size_t n = m_coarse2fine.shape(1); xt::xtensor ret = xt::empty({m, n}); for (size_t i = 0; i < m; ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); xt::view(ret, i) = xt::view(data, xt::keep(e)); } return ret; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); size_t m = m_coarse2fine.shape(0); size_t n = m_coarse2fine.shape(1); size_t N = data.shape(1); xt::xtensor ret = xt::empty({m, n * N}); for (size_t i = 0; i < m; ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); for (size_t q = 0; q < data.shape(1); ++q) { xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q); } } return ret; } inline xt::xtensor RefineRegular::mapToCoarse(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size()); size_t m = m_coarse2fine.shape(0); size_t n = m_coarse2fine.shape(1); size_t N = data.shape(1); xt::xtensor ret = xt::empty({m, n * N, data.shape(2), data.shape(3)}); for (size_t i = 0; i < m; ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); for (size_t q = 0; q < data.shape(1); ++q) { xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q); } } return ret; } inline xt::xtensor RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({m_coarse2fine.size()}); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); xt::view(ret, xt::keep(e)) = data(i); } return ret; } inline xt::xtensor RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({m_coarse2fine.size(), data.shape(1)}); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); xt::view(ret, xt::keep(e)) = xt::view(data, i); } return ret; } inline xt::xtensor RefineRegular::mapToFine(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0)); xt::xtensor ret = xt::empty({m_coarse2fine.size(), data.shape(1), data.shape(2), data.shape(3)}); for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) { auto e = xt::view(m_coarse2fine, i, xt::all()); xt::view(ret, xt::keep(e)) = xt::view(data, i); } return ret; } inline FineLayer2Regular::FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh) : m_finelayer(mesh) { // ------------ // Regular-mesh // ------------ m_regular = GooseFEM::Mesh::Quad4::Regular( xt::amax(m_finelayer.m_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h); // ------- // mapping // ------- // allocate mapping m_elem_regular.resize(m_finelayer.m_nelem); m_frac_regular.resize(m_finelayer.m_nelem); // alias xt::xtensor nhx = m_finelayer.m_nhx; xt::xtensor nhy = m_finelayer.m_nhy; xt::xtensor nelx = m_finelayer.m_nelx; xt::xtensor start = m_finelayer.m_startElem; // 'matrix' of element numbers of the Regular-mesh xt::xtensor elementgrid = m_regular.elementgrid(); // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh xt::xtensor cum_nhy = xt::concatenate(xt::xtuple(xt::zeros({1}), xt::cumsum(nhy))); // number of element layers in y-direction of the FineLayer-mesh size_t nely = nhy.size(); // loop over layers of the FineLayer-mesh for (size_t iy = 0; iy < nely; ++iy) { // element numbers of the Regular-mesh along this layer of the FineLayer-mesh auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all()); // no coarsening/refinement // ------------------------ if (m_finelayer.m_refine(iy) == -1) { // element numbers of the FineLayer-mesh along this layer xt::xtensor el_old = start(iy) + xt::arange(nelx(iy)); // loop along this layer of the FineLayer-mesh for (size_t ix = 0; ix < nelx(iy); ++ix) { // get the element numbers of the Regular-mesh for this element of the // FineLayer-mesh auto block = xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy))); // write to mapping for (auto& i : block) { m_elem_regular[el_old(ix)].push_back(i); m_frac_regular[el_old(ix)].push_back(1.0); } } } // refinement along the x-direction (below the middle layer) // --------------------------------------------------------- else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) { // element numbers of the FineLayer-mesh along this layer // rows: coarse block, columns element numbers per block xt::xtensor el_old = start(iy) + xt::arange(nelx(iy) * 4ul).reshape({-1, 4}); // loop along this layer of the FineLayer-mesh for (size_t ix = 0; ix < nelx(iy); ++ix) { // get the element numbers of the Regular-mesh for this block of the FineLayer-mesh auto block = xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy))); // bottom: wide-to-narrow { for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, j, xt::range(j, nhx(iy) - j)); m_elem_regular[el_old(ix, 0)].push_back(e(0)); m_frac_regular[el_old(ix, 0)].push_back(0.5); for (size_t k = 1; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 0)].push_back(e(k)); m_frac_regular[el_old(ix, 0)].push_back(1.0); } m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 0)].push_back(0.5); } } // top: regular small element { auto e = xt::view( block, xt::range(nhy(iy) / 2, nhy(iy)), xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 2)].push_back(i); m_frac_regular[el_old(ix, 2)].push_back(1.0); } } // left { // left-bottom: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, j, xt::range(0, j + 1)); for (size_t k = 0; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 3)].push_back(e(k)); m_frac_regular[el_old(ix, 3)].push_back(1.0); } m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 3)].push_back(0.5); } // left-top: regular { auto e = xt::view( block, xt::range(nhy(iy) / 2, nhy(iy)), xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 3)].push_back(i); m_frac_regular[el_old(ix, 3)].push_back(1.0); } } } // right { // right-bottom: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy))); m_elem_regular[el_old(ix, 1)].push_back(e(0)); m_frac_regular[el_old(ix, 1)].push_back(0.5); for (size_t k = 1; k < e.size(); ++k) { m_elem_regular[el_old(ix, 1)].push_back(e(k)); m_frac_regular[el_old(ix, 1)].push_back(1.0); } } // right-top: regular { auto e = xt::view( block, xt::range(nhy(iy) / 2, nhy(iy)), xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 1)].push_back(i); m_frac_regular[el_old(ix, 1)].push_back(1.0); } } } } } // coarsening along the x-direction (above the middle layer) else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) { // element numbers of the FineLayer-mesh along this layer // rows: coarse block, columns element numbers per block xt::xtensor el_old = start(iy) + xt::arange(nelx(iy) * 4ul).reshape({-1, 4}); // loop along this layer of the FineLayer-mesh for (size_t ix = 0; ix < nelx(iy); ++ix) { // get the element numbers of the Regular-mesh for this block of the FineLayer-mesh auto block = xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy))); // top: narrow-to-wide { for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view( block, nhy(iy) / 2 + j, xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1)); m_elem_regular[el_old(ix, 3)].push_back(e(0)); m_frac_regular[el_old(ix, 3)].push_back(0.5); for (size_t k = 1; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 3)].push_back(e(k)); m_frac_regular[el_old(ix, 3)].push_back(1.0); } m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 3)].push_back(0.5); } } // bottom: regular small element { auto e = xt::view( block, xt::range(0, nhy(iy) / 2), xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 1)].push_back(i); m_frac_regular[el_old(ix, 1)].push_back(1.0); } } // left { // left-bottom: regular { auto e = xt::view( block, xt::range(0, nhy(iy) / 2), xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 0)].push_back(i); m_frac_regular[el_old(ix, 0)].push_back(1.0); } } // left-top: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j)); for (size_t k = 0; k < e.size() - 1; ++k) { m_elem_regular[el_old(ix, 0)].push_back(e(k)); m_frac_regular[el_old(ix, 0)].push_back(1.0); } m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1)); m_frac_regular[el_old(ix, 0)].push_back(0.5); } } // right { // right-bottom: regular { auto e = xt::view( block, xt::range(0, nhy(iy) / 2), xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3)); for (auto& i : e) { m_elem_regular[el_old(ix, 2)].push_back(i); m_frac_regular[el_old(ix, 2)].push_back(1.0); } } // right-top: narrow-to-wide for (size_t j = 0; j < nhy(iy) / 2; ++j) { auto e = xt::view( block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy))); m_elem_regular[el_old(ix, 2)].push_back(e(0)); m_frac_regular[el_old(ix, 2)].push_back(0.5); for (size_t k = 1; k < e.size(); ++k) { m_elem_regular[el_old(ix, 2)].push_back(e(k)); m_frac_regular[el_old(ix, 2)].push_back(1.0); } } } } } } } inline GooseFEM::Mesh::Quad4::Regular FineLayer2Regular::getRegularMesh() const { return m_regular; } inline GooseFEM::Mesh::Quad4::FineLayer FineLayer2Regular::getFineLayerMesh() const { return m_finelayer; } inline std::vector> FineLayer2Regular::getMap() const { return m_elem_regular; } inline std::vector> FineLayer2Regular::getMapFraction() const { return m_frac_regular; } template -inline xt::xtensor -FineLayer2Regular::mapToRegular(const xt::xtensor& data) const +inline xt::xtensor FineLayer2Regular::mapToRegular(const xt::xtensor& data) const { GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem()); std::array shape; std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]); shape[0] = m_regular.nelem(); xt::xtensor ret = xt::zeros(shape); for (size_t e = 0; e < m_finelayer.nelem(); ++e) { for (size_t i = 0; i < m_elem_regular[e].size(); ++i) { xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e); } } return ret; } } // namespace Map } // namespace Quad4 } // namespace Mesh } // namespace GooseFEM #endif diff --git a/include/GooseFEM/Vector.h b/include/GooseFEM/Vector.h index 5f4d81c..6ad79cb 100644 --- a/include/GooseFEM/Vector.h +++ b/include/GooseFEM/Vector.h @@ -1,110 +1,111 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTOR_H #define GOOSEFEM_VECTOR_H #include "config.h" namespace GooseFEM { -/* - "nodevec" - nodal vectors - [nnode, ndim] - "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] - "dofval" - DOF values - [ndof] -*/ +/** +Class to switch between: +- "nodevec": nodal vectors, shape ``[nnode, ndim]``. +- "elemvec": nodal vectors stored per element, shape: ``[nelem, nne, ndim]``. +- "dofval": DOF values, shape: ``[ndof]``. +*/ class Vector { public: Vector() = default; /** Constructor. \param conn Connectivity, shape ``[nelem, nne]``. \param dofs DOFs per node, shape ``[nnode, ndim]``. */ Vector(const xt::xtensor& conn, const xt::xtensor& dofs); // Dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs // DOF lists xt::xtensor dofs() const; // DOFs // Copy nodevec to another nodevec void copy(const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const; // Convert to "dofval" (overwrite entries that occur more than once) -- (auto allocation below) void asDofs(const xt::xtensor& nodevec, xt::xtensor& dofval) const; void asDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const; // Convert to "nodevec" (overwrite entries that occur more than once) -- (auto allocation below) void asNode(const xt::xtensor& dofval, xt::xtensor& nodevec) const; void asNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const; // Convert to "elemvec" (overwrite entries that occur more than once) -- (auto allocation below) void asElement(const xt::xtensor& dofval, xt::xtensor& elemvec) const; void asElement(const xt::xtensor& nodevec, xt::xtensor& elemvec) const; // Assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) void assembleDofs(const xt::xtensor& nodevec, xt::xtensor& dofval) const; void assembleDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const; // Assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) void assembleNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const; // Auto-allocation of the functions above xt::xtensor AsDofs(const xt::xtensor& nodevec) const; xt::xtensor AsDofs(const xt::xtensor& elemvec) const; xt::xtensor AsNode(const xt::xtensor& dofval) const; xt::xtensor AsNode(const xt::xtensor& elemvec) const; xt::xtensor AsElement(const xt::xtensor& dofval) const; xt::xtensor AsElement(const xt::xtensor& nodevec) const; xt::xtensor AssembleDofs(const xt::xtensor& nodevec) const; xt::xtensor AssembleDofs(const xt::xtensor& elemvec) const; xt::xtensor AssembleNode(const xt::xtensor& elemvec) const; xt::xtensor Copy( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const; // Get shape of dofval, nodevec, elemvec std::array ShapeDofval() const; std::array ShapeNodevec() const; std::array ShapeElemvec() const; std::array ShapeElemmat() const; // Get zero-allocated dofval, nodevec, elemvec xt::xtensor AllocateDofval() const; xt::xtensor AllocateNodevec() const; xt::xtensor AllocateElemvec() const; xt::xtensor AllocateElemmat() const; xt::xtensor AllocateDofval(double val) const; xt::xtensor AllocateNodevec(double val) const; xt::xtensor AllocateElemvec(double val) const; xt::xtensor AllocateElemmat(double val) const; protected: xt::xtensor m_conn; ///< Connectivity ``[nelem, nne]`` xt::xtensor m_dofs; ///< DOF-numbers per node ``[nnode, ndim]`` size_t m_nelem; ///< Number of elements size_t m_nne; ///< Number of nodes per element size_t m_nnode; ///< Number of nodes size_t m_ndim; ///< Number of dimensions size_t m_ndof; ///< Number of DOFs }; } // namespace GooseFEM #include "Vector.hpp" #endif