diff --git a/include/GooseFEM/Element.h b/include/GooseFEM/Element.h index 26c6760..38c4906 100644 --- a/include/GooseFEM/Element.h +++ b/include/GooseFEM/Element.h @@ -1,630 +1,630 @@ /** Convenience methods for integration point data. \file Element.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENT_H #define GOOSEFEM_ELEMENT_H #include "config.h" #include "Allocate.h" namespace GooseFEM { /** Element quadrature and interpolation. */ namespace Element { /** Convert nodal vector with ("nodevec", shape:``[nnode, ndim]``) to nodal vector stored per element ("elemvec", shape: ``[nelem, nne, ndim]``). \param conn Connectivity. \param nodevec "nodevec". \return "elemvec". */ inline xt::xtensor asElementVector( const xt::xtensor& conn, const xt::xtensor& nodevec); /** Assemble nodal vector stored per element ("elemvec", shape ``[nelem, nne, ndim]``) to nodal vector ("nodevec", shape ``[nnode, ndim]``). \param conn Connectivity. \param elemvec "elemvec". \return "nodevec". */ inline xt::xtensor assembleNodeVector( const xt::xtensor& conn, const xt::xtensor& elemvec); /** Check that DOFs leave no holes. \param dofs DOFs ("nodevec") \return ``true`` if there are no holds. */ template inline bool isSequential(const E& dofs); /** Check that all of the matrices stored per elemmat (shape: ``[nelem, nne * ndim, nne * ndim]``) are diagonal. \param elemmat Element-vectors ("elemmat") \return ``true`` if all element matrices are diagonal. */ bool isDiagonal(const xt::xtensor& elemmat); /** Base quadrature-class. This class does not have a specific element-type in mind, it is used mostly internally to derive from such that common methods do not have to be reimplementation. \tparam ne Number of nodes per element. \tparam nd Number of dimensions for node vectors. \tparam td Number of dimensions for integration point tensors. */ template class QuadratureBase { public: QuadratureBase() = default; /** Constructor */ QuadratureBase(size_t nelem, size_t nip); /** Number of elements. \return Scalar. */ size_t nelem() const; /** Number of nodes per element. \return Scalar. */ size_t nne() const; /** Number of dimensions for node vectors. \return Scalar. */ size_t ndim() const; /** Number of dimensions for integration point tensors. \return Scalar. */ size_t tdim() const; /** Number of integration points. \return Scalar. */ size_t nip() const; /** Convert "qscalar" to "qtensor" of certain rank. Fully allocated output passed as reference, use AsTensor to allocate and return data. \param qscalar A "qscalar". \param qtensor A "qtensor". */ template void asTensor(const xt::xtensor& qscalar, xt::xtensor& qtensor) const; /** Convert "qscalar" to "qtensor" of certain rank. \param "qscalar". \return "qtensor". */ template xt::xtensor AsTensor(const xt::xtensor& qscalar) const; /** Convert "qscalar" to "qtensor" of certain rank. \param rank Tensor rank. \param qscalar A "qscalar". \return "qtensor". */ template xt::xarray AsTensor(size_t rank, const xt::xtensor& qscalar) const; /** Get the shape of an "elemvec". \returns [#nelem, #nne, #ndim]. */ std::array shape_elemvec() const; /** Get the shape of an "elemmat". \returns [#nelem, #nne * #ndim, #nne * #ndim]. */ std::array shape_elemmat() const; /** Get the shape of a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). Default: rank = 0, a.k.a. scalar. \returns [#nelem, #nip, #tdim, ...]. */ template std::array shape_qtensor() const; /** Get the shape of a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). \param rank The tensor rank. \returns [#nelem, #nip, #tdim, ...]. */ std::vector shape_qtensor(size_t rank) const; /** Get the shape of a "qscalar" (a "qtensor" of rank 0) \returns [#nelem, #nip]. */ std::vector shape_qscalar() const; /** Get an allocated `xt::xtensor` to store a "elemvec". Note: the container is not (zero-)initialised. \returns `xt::xarray` container of the correct shape. */ template xt::xtensor allocate_elemvec() const; /** Get an allocated and initialised `xt::xarray` to store a "elemvec". \param val The value to which to initialise all items. \returns `xt::xtensor` container of the correct shape. */ template xt::xtensor allocate_elemvec(T val) const; /** Get an allocated `xt::xtensor` to store a "elemmat". Note: the container is not (zero-)initialised. \returns `xt::xarray` container of the correct shape. */ template xt::xtensor allocate_elemmat() const; /** Get an allocated and initialised `xt::xarray` to store a "elemmat". \param val The value to which to initialise all items. \returns `xt::xtensor` container of the correct shape. */ template xt::xtensor allocate_elemmat(T val) const; /** Get an allocated `xt::xtensor` to store a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). Default: rank = 0, a.k.a. scalar. Note: the container is not (zero-)initialised. \returns [#nelem, #nip]. */ template xt::xtensor allocate_qtensor() const; /** Get an allocated and initialised `xt::xtensor` to store a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). Default: rank = 0, a.k.a. scalar. \param val The value to which to initialise all items. \returns `xt::xtensor` container of the correct shape (and rank). */ template xt::xtensor allocate_qtensor(T val) const; /** Get an allocated `xt::xarray` to store a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). Note: the container is not (zero-)initialised. \param rank The tensor rank. \returns `xt::xarray` container of the correct shape. */ template xt::xarray allocate_qtensor(size_t rank) const; /** Get an allocated and initialised `xt::xarray` to store a "qtensor" of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.). \param rank The tensor rank. \param val The value to which to initialise all items. \returns `xt::xtensor` container of the correct shape (and rank). */ template xt::xarray allocate_qtensor(size_t rank, T val) const; /** Get an allocated `xt::xtensor` to store a "qscalar" (a "qtensor" of rank 0). Note: the container is not (zero-)initialised. \returns `xt::xarray` container of the correct shape. */ template xt::xtensor allocate_qscalar() const; /** Get an allocated and initialised `xt::xarray` to store a "qscalar" (a "qtensor" of rank 0). \param val The value to which to initialise all items. \returns `xt::xtensor` container of the correct shape (and rank). */ template xt::xtensor allocate_qscalar(T val) const; /** \cond */ template [[ deprecated ]] std::array ShapeQtensor() const; [[ deprecated ]] std::vector ShapeQtensor(size_t rank) const; template [[ deprecated ]] xt::xtensor AllocateQtensor() const; [[ deprecated ]] std::vector ShapeQscalar() const; template [[ deprecated ]] xt::xtensor AllocateQtensor(T val) const; template [[ deprecated ]] xt::xarray AllocateQtensor(size_t rank) const; template [[ deprecated ]] xt::xarray AllocateQtensor(size_t rank, T val) const; template [[ deprecated ]] xt::xtensor AllocateQscalar() const; template [[ deprecated ]] xt::xtensor AllocateQscalar(T val) const; /** \endcond */ protected: /** Wrapper of constructor, for derived classes. */ void initQuadratureBase(size_t nelem, size_t nip); protected: size_t m_nelem; ///< Number of elements. size_t m_nip; ///< Number of integration points per element. constexpr static size_t m_nne = ne; ///< Number of nodes per element. constexpr static size_t m_ndim = nd; ///< Number of dimensions for nodal vectors. constexpr static size_t m_tdim = td; ///< Number of dimensions for integration point tensors. }; /** Interpolation and quadrature for a generic element in Cartesian coordinates. 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, #tdim, #tdim] - ``qscalar``: integration point scalar, [#nelem, #nip] */ template class QuadratureBaseCartesian : public QuadratureBase { public: QuadratureBaseCartesian() = default; virtual ~QuadratureBaseCartesian(){}; /** Constructor with custom integration. During construction the values of the shape functions and the shape function gradients (in local and global) coordinates are computed. They can be reused without any cost. They only have to be recomputed when the nodal position changes. In that case use update_x() to update the nodal positions and to recompute the shape functions and their gradients. \param x nodal coordinates (``elemvec``). \param xi Integration point coordinates (local coordinates) [#nip]. \param w Integration point weights [#nip]. \param N Shape functions in local coordinates [#nip, #nne]. \param dNdxi Shape function gradient w.r.t. local coordinates [#nip, #nne, #ndim]. */ QuadratureBaseCartesian( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w, const xt::xtensor& N, const xt::xtensor& dNdxi); /** Update the nodal positions. This recomputes the values of the shape functions and the shape function gradients (in local and global) coordinates. Under the small deformations assumption this should not be called. \param x nodal coordinates (``elemvec``). Shape should match the earlier definition. */ void update_x(const xt::xtensor& x); /** Get the shape function gradients (in global coordinates). - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \return ``gradN`` stored per element, per integration point [#nelem, #nip, #nne, #ndim]. */ - xt::xtensor GradN() const; + virtual xt::xtensor GradN() const; /** Get the integration volume. - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \return volume stored per element, per integration point [#nelem, #nip]. */ xt::xtensor dV() const; /** Interpolate element vector. - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \param elemvec nodal vector stored per element (shape: [#nelem, #nne, #ndim]). \return integration point vector (shape: [#nelem, #nip, #ndim]). */ xt::xtensor Interp_N_vector(const xt::xtensor& elemvec) const; /** Same as Interp_N_vector(), but writing to preallocated return. \param elemvec nodal vector stored per element (shape: [#nelem, #nne, #ndim]). \param qvector integration point vector, overwritten (shape: [#nelem, #nip, #ndim]). */ virtual void interp_N_vector( const xt::xtensor& elemvec, xt::xtensor& qvector) const; /** Element-by-element: dyadic product of the shape function gradients and a nodal vector. Typical input: nodal displacements. Typical output: quadrature point strains. Within one element:: for e in range(nelem): for q in range(nip): for m in range(nne): qtensor(e, q, i, j) += dNdx(e, q, m, i) * elemvec(e, m, j) - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \param elemvec [#nelem, #nne, #ndim] \return qtensor [#nelem, #nip, #tdim, #tdim] */ xt::xtensor GradN_vector(const xt::xtensor& elemvec) const; /** Same as GradN_vector(), but writing to preallocated return. \param elemvec [#nelem, #nne, #ndim] \param qtensor overwritten [#nelem, #nip, #tdim, #tdim] */ virtual void gradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const; /** The transposed output of GradN_vector(). Within one element:: for e in range(nelem): for q in range(nip): for m in range(nne): qtensor(e, q, j, i) += dNdx(e, q, m, i) * elemvec(e, m, j) \param elemvec [#nelem, #nne, #ndim] \return qtensor [#nelem, #nip, #tdim, #tdim] */ xt::xtensor GradN_vector_T(const xt::xtensor& elemvec) const; /** Same as GradN_vector_T(), but writing to preallocated return. \param elemvec [#nelem, #nne, #ndim] \param qtensor overwritten [#nelem, #nip, #tdim, #tdim] */ virtual void gradN_vector_T( const xt::xtensor& elemvec, xt::xtensor& qtensor) const; /** The symmetric output of GradN_vector(). Without one element:: for e in range(nelem): for q in range(nip): for m in range(nne): qtensor(e, q, i, j) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j) qtensor(e, q, j, i) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j) \param elemvec [#nelem, #nne, #ndim] \return qtensor [#nelem, #nip, #tdim, #tdim] */ xt::xtensor SymGradN_vector(const xt::xtensor& elemvec) const; /** Same as SymGradN_vector(), but writing to preallocated return. \param elemvec [#nelem, #nne, #ndim] \param qtensor overwritten [#nelem, #nip, #tdim, #tdim] */ virtual void symGradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const; /** Element-by-element: integral of the scalar product of the shape function with a scalar. Within one one element:: for e in range(nelem): for q in range(nip): for m in range(nne): for n in range(nne): elemmat(e, m * ndim + i, n * ndim + i) += N(e, q, m) * qscalar(e, q) * N(e, q, n) * dV(e, q) with ``i`` a tensor dimension. - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \param qscalar [#nelem, #nip] \return elemmat [#nelem, #nne * #ndim, #nne * #ndim] */ xt::xtensor Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const; /** Same as Int_N_scalar_NT_dV(), but writing to preallocated return. \param qscalar [#nelem, #nip] \param elemmat overwritten [#nelem, #nne * #ndim, #nne * #ndim] */ virtual void int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const; /** Element-by-element: integral of the dot product of the shape function gradients with a second order tensor. Typical input: stress. Typical output: nodal force. Within one one element:: for e in range(nelem): for q in range(nip): for m in range(nne): elemvec(e, m, j) += dNdx(e, q, m, i) * qtensor(e, q, i, j) * dV(e, q) with ``i`` and ``j`` tensor dimensions. - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \param qtensor [#nelem, #nip, #ndim, #ndim] \return elemvec [#nelem, #nne. #ndim] */ xt::xtensor Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const; /** Same as Int_gradN_dot_tensor2_dV(), but writing to preallocated return. \param qtensor [#nelem, #nip, #ndim, #ndim] \param elemvec overwritten [#nelem, #nne. #ndim] */ virtual 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 /** Element-by-element: integral of the dot products of the shape function gradients with a fourth order tensor. Typical input: stiffness tensor. Typical output: stiffness matrix. Within one one element:: for e in range(nelem): for q in range(nip): for m in range(nne): for n in range(nne): elemmat(e, m * ndim + j, n * ndim + k) += dNdx(e, q, m, i) * qtensor(e, q, i, j, k, l) * dNdx(e, q, n, l) * dV(e, q) with ``i``, ``j``, ``k``, and ``l`` tensor dimensions. - Note that the function and their gradient are precomputed upon construction, + Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x(). \param qtensor [#nelem, #nip, #ndim, #ndim, #ndim, #ndim] \return elemmat [#nelem, #nne * #ndim, #nne * #ndim] */ xt::xtensor Int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor) const; /** Same as Int_gradN_dot_tensor4_dot_gradNT_dV(), but writing to preallocated return. \param qtensor [#nelem, #nip, #ndim, #ndim, #ndim, #ndim] \param elemmat overwritten [#nelem, #nne * #ndim, #nne * #ndim] */ virtual void int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const; protected: /** Constructor alias. \param x nodal coordinates (``elemvec``). \param xi Integration point coordinates (local coordinates) [#nip]. \param w Integration point weights [#nip]. \param N Shape functions in the integration points [#nip, #nne]. \param dNdxi Shape function gradient w.r.t. local coordinates [#nip, #nne, #ndim]. */ void initQuadratureBaseCartesian( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w, const xt::xtensor& N, const xt::xtensor& dNdxi); /** Compute m_vol and m_dNdx based on current m_x. */ virtual void compute_dN(); protected: using QuadratureBase::m_nelem; using QuadratureBase::m_nip; using QuadratureBase::m_nne; using QuadratureBase::m_ndim; using QuadratureBase::m_tdim; 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 Element } // namespace GooseFEM #include "Element.hpp" #endif diff --git a/include/GooseFEM/ElementQuad4Axisymmetric.h b/include/GooseFEM/ElementQuad4Axisymmetric.h index af44e4c..a601b21 100644 --- a/include/GooseFEM/ElementQuad4Axisymmetric.h +++ b/include/GooseFEM/ElementQuad4Axisymmetric.h @@ -1,99 +1,93 @@ /** Quadrature for 4-noded quadrilateral element in 2d (GooseFEM::Mesh::ElementType::Quad4), in an axisymmetric coordinated system. \file ElementQuad4Axisymmetric.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENTQUAD4AXISYMMETRIC_H #define GOOSEFEM_ELEMENTQUAD4AXISYMMETRIC_H #include "config.h" namespace GooseFEM { namespace Element { namespace Quad4 { -class QuadratureAxisymmetric : public GooseFEM::Element::QuadratureBase<4, 2, 3> { +/** +Interpolation and quadrature. + +Fixed dimensions: +- ``ndim = 2``: number of dimensions. +- ``tdim = 3``: number of dimensions or tensor. +- ``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, #tdim, #tdim] +- ``qscalar``: integration point scalar, [#nelem, #nip] +*/ +class QuadratureAxisymmetric : public GooseFEM::Element::QuadratureBaseCartesian<4, 2, 3> { public: - // Fixed dimensions: - // ndim = 2 - number of dimensions - // nne = 4 - number of nodes per element - // tdim = 3 - number of dimensions of tensors - // - // 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, tdim, tdim] - // "qscalar" - integration point scalar - [nelem, nip] - - // Constructor: integration point coordinates and weights are optional (default: Gauss) + QuadratureAxisymmetric() = default; QuadratureAxisymmetric(const xt::xtensor& x); QuadratureAxisymmetric( 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); + /** + Get the B-matrix (shape function gradients) (in global coordinates). + Note that the functions and their gradients are precomputed upon construction, + or updated when calling update_x(). + + \return ``B`` matrix stored per element, per integration point [#nelem, #nne, #tdim, #tdim, #tdim] + */ + xt::xtensor B() const; - xt::xtensor GradN() const; + // qtensor(e, q, i, j) += B(e, q, m, i, j, k) * elemvec(e, q, m, k) + void gradN_vector( + const xt::xtensor& elemvec, xt::xtensor& qtensor) const override; - // Return integration volume - xt::xtensor dV() const; + void gradN_vector_T( + const xt::xtensor& elemvec, xt::xtensor& qtensor) const override; - // Dyadic product (and its transpose and symmetric part) - // qtensor(i,j) += B(m,i,j,k) * elemvec(m,k) - 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; + void symGradN_vector( + const xt::xtensor& elemvec, xt::xtensor& qtensor) const override; - // Integral of the scalar product - // elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV + // elemmat(e, q, m * ndim + i, n * ndim + i) += + // N(e, q, m) * qscalar(e, q) * N(e, q, n) * dV(e, q) void int_N_scalar_NT_dV( - const xt::xtensor& qscalar, xt::xtensor& elemmat) const; + const xt::xtensor& qscalar, xt::xtensor& elemmat) const override; - // Integral of the assembled product // fm = ( Bm^T : qtensor ) dV void int_gradN_dot_tensor2_dV( - const xt::xtensor& qtensor, xt::xtensor& elemvec) const; + const xt::xtensor& qtensor, xt::xtensor& elemvec) const override; - // Integral of the assembled product // Kmn = ( Bm^T : qtensor : Bn ) dV void int_gradN_dot_tensor4_dot_gradNT_dV( - const xt::xtensor& qtensor, xt::xtensor& elemmat) const; + const xt::xtensor& qtensor, xt::xtensor& elemmat) const override; - // 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; +protected: + void compute_dN() override; private: - // Compute "vol" and "B" based on current "x" - void compute_dN(); + xt::xtensor GradN() const override; 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_B; // B-matrix [nelem, nne, tdim, tdim, tdim] - xt::xtensor m_vol; // integration point volume [nelem, nip] + xt::xtensor m_B; ///< B-matrix [#nelem, #nne, #tdim, #tdim, #tdim] }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #include "ElementQuad4Axisymmetric.hpp" #endif diff --git a/include/GooseFEM/ElementQuad4Axisymmetric.hpp b/include/GooseFEM/ElementQuad4Axisymmetric.hpp index 9f3a658..31621f3 100644 --- a/include/GooseFEM/ElementQuad4Axisymmetric.hpp +++ b/include/GooseFEM/ElementQuad4Axisymmetric.hpp @@ -1,432 +1,385 @@ /** Implementation of ElementQuad4Axisymmetric.h \file ElementQuad4Axisymmetric.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ELEMENTQUADAXISYMMETRIC_HPP #define GOOSEFEM_ELEMENTQUADAXISYMMETRIC_HPP #include "ElementQuad4Axisymmetric.h" #include "detail.hpp" namespace GooseFEM { namespace Element { namespace Quad4 { inline QuadratureAxisymmetric::QuadratureAxisymmetric(const xt::xtensor& x) : QuadratureAxisymmetric(x, Gauss::xi(), Gauss::w()) { } inline QuadratureAxisymmetric::QuadratureAxisymmetric( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w) - : m_x(x), m_w(w), m_xi(xi) { + m_x = x; + m_w = w; + m_xi = xi; + this->initQuadratureBase(m_x.shape(0), m_w.size()); + GOOSEFEM_ASSERT(m_x.shape(1) == m_nne); GOOSEFEM_ASSERT(m_x.shape(2) == m_ndim); - GOOSEFEM_ASSERT(m_xi.shape(0) == m_nip); - GOOSEFEM_ASSERT(m_xi.shape(1) == m_ndim); - GOOSEFEM_ASSERT(m_w.size() == m_nip); + GOOSEFEM_ASSERT(xt::has_shape(m_xi, {m_nip, m_ndim})); + GOOSEFEM_ASSERT(xt::has_shape(m_w, {m_nip})); m_N = xt::empty({m_nip, m_nne}); m_dNxi = xt::empty({m_nip, m_nne, m_ndim}); m_B = xt::empty({m_nelem, m_nip, m_nne, m_tdim, m_tdim, m_tdim}); - m_vol = xt::empty({m_nelem, m_nip}); + m_vol = xt::empty(this->shape_qscalar()); for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 1) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 2) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_N(q, 3) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 1, 0) = +0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 2, 0) = +0.25 * (1.0 + m_xi(q, 1)); m_dNxi(q, 3, 0) = -0.25 * (1.0 + m_xi(q, 1)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.25 * (1.0 - m_xi(q, 0)); m_dNxi(q, 1, 1) = -0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 2, 1) = +0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 3, 1) = +0.25 * (1.0 - m_xi(q, 0)); } - compute_dN(); -} - -inline xt::xtensor QuadratureAxisymmetric::dV() const -{ - return m_vol; -} - -inline void QuadratureAxisymmetric::update_x(const xt::xtensor& x) -{ - GOOSEFEM_ASSERT(x.shape() == m_x.shape()); - xt::noalias(m_x) = x; - compute_dN(); + this->compute_dN(); } inline void QuadratureAxisymmetric::compute_dN() { // most components remain zero, and are not written m_B.fill(0.0); #pragma omp parallel { xt::xtensor J = xt::empty({2, 2}); xt::xtensor Jinv = xt::empty({2, 2}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto N = xt::adapt(&m_N(q, 0), xt::xshape()); // J(i,j) += dNxi(m,i) * x(m,j); J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) + dNxi(3, 0) * x(3, 0); J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) + dNxi(3, 0) * x(3, 1); J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) + dNxi(3, 1) * x(3, 0); J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) + dNxi(3, 1) * x(3, 1); double Jdet = detail::tensor<2>::inv(J, Jinv); // radius for computation of volume double rq = N(0) * x(0, 1) + N(1) * x(1, 1) + N(2) * x(2, 1) + N(3) * x(3, 1); // dNx(m,i) += Jinv(i,j) * dNxi(m,j) for (size_t m = 0; m < m_nne; ++m) { // B(m, r, r, r) = dNdx(m,1) B(m, 0, 0, 0) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); // B(m, r, z, z) = dNdx(m,1) B(m, 0, 2, 2) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); // B(m, t, t, r) B(m, 1, 1, 0) = 1.0 / rq * N(m); // B(m, z, r, r) = dNdx(m,0) B(m, 2, 0, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); // B(m, z, z, z) = dNdx(m,0) B(m, 2, 2, 2) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); } m_vol(e, q) = m_w(q) * Jdet * 2.0 * M_PI * rq; } } } } +inline xt::xtensor QuadratureAxisymmetric::GradN() const +{ + return m_dNx; // empty: but function is private to cannot be called +} + +inline xt::xtensor QuadratureAxisymmetric::B() const +{ + return m_B; +} + inline void QuadratureAxisymmetric::gradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += B(m,i,j,k) * u(m,perm(k)) // (where perm(0) = 1, perm(2) = 0) gradu(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); gradu(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); gradu(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); gradu(0, 2) = B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0); gradu(2, 0) = B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1); } } } inline void QuadratureAxisymmetric::gradN_vector_T( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += B(m,i,j,k) * u(m,perm(k)) // (where perm(0) = 1, perm(2) = 0) gradu(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); gradu(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); gradu(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); gradu(2, 0) = B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0); gradu(0, 2) = B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1); } } } inline void QuadratureAxisymmetric::symGradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += B(m,i,j,k) * u(m,perm(k)) // eps(j,i) = 0.5 * (gradu(i,j) + gradu(j,i)) // (where perm(0) = 1, perm(2) = 0) eps(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); eps(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); eps(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); eps(2, 0) = 0.5 * (B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0) + B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1)); eps(0, 2) = eps(2, 0); } } } inline void QuadratureAxisymmetric::int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, {m_nelem, m_nip})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { M(m * m_ndim + 0, n * m_ndim + 0) += N(m) * rho * N(n) * vol; M(m * m_ndim + 1, n * m_ndim + 1) += N(m) * rho * N(n) * vol; } } } } } inline void QuadratureAxisymmetric::int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); // f(m,i) += B(m,i,j,perm(k)) * sig(i,j) * dV // (where perm(0) = 1, perm(2) = 0) for (size_t m = 0; m < m_nne; ++m) { f(m, 0) += vol * (B(m, 2, 2, 2) * sig(2, 2) + B(m, 0, 2, 2) * sig(0, 2)); f(m, 1) += vol * (B(m, 0, 0, 0) * sig(0, 0) + B(m, 1, 1, 0) * sig(1, 1) + B(m, 2, 0, 0) * sig(2, 0)); } } } } inline void QuadratureAxisymmetric::int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim, m_tdim, m_tdim})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto C = xt::adapt(&qtensor(e, q, 0, 0, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); // K(m*m_ndim+perm(c), n*m_ndim+perm(f)) = B(m,a,b,c) * C(a,b,d,e) * B(n,e,d,f) * vol; // (where perm(0) = 1, perm(2) = 0) for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 0, 0, 0) * C(0, 0, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 0, 0, 0) * C(0, 0, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 1, 1, 0) * C(1, 1, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 1, 1, 0) * C(1, 1, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 2, 2, 2) * C(2, 2, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 2, 2, 2) * C(2, 2, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 0, 2, 2) * C(0, 2, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 0, 2, 2) * C(0, 2, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 2, 0, 0) * C(2, 0, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 2, 0, 0) * C(2, 0, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 0, 2) * B(n, 2, 0, 0) * vol; } } } } } -inline xt::xtensor -QuadratureAxisymmetric::GradN_vector(const xt::xtensor& elemvec) const -{ - xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); - this->gradN_vector(elemvec, qtensor); - return qtensor; -} - -inline xt::xtensor -QuadratureAxisymmetric::GradN_vector_T(const xt::xtensor& elemvec) const -{ - xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); - this->gradN_vector_T(elemvec, qtensor); - return qtensor; -} - -inline xt::xtensor -QuadratureAxisymmetric::SymGradN_vector(const xt::xtensor& elemvec) const -{ - xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); - this->symGradN_vector(elemvec, qtensor); - return qtensor; -} - -inline xt::xtensor -QuadratureAxisymmetric::Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const -{ - xt::xtensor elemmat = xt::empty({m_nelem, m_nne * m_ndim, m_nne * m_ndim}); - this->int_N_scalar_NT_dV(qscalar, elemmat); - return elemmat; -} - -inline xt::xtensor -QuadratureAxisymmetric::Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - this->int_gradN_dot_tensor2_dV(qtensor, elemvec); - return elemvec; -} - -inline xt::xtensor -QuadratureAxisymmetric::Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const -{ - xt::xtensor elemmat = xt::empty({m_nelem, m_ndim * m_nne, m_ndim * m_nne}); - this->int_gradN_dot_tensor4_dot_gradNT_dV(qtensor, elemmat); - return elemmat; -} - } // namespace Quad4 } // namespace Element } // namespace GooseFEM #endif