diff --git a/include/xGooseFEM/ElementHex8.h b/include/xGooseFEM/ElementHex8.h index 4c4b25b..f65efbe 100644 --- a/include/xGooseFEM/ElementHex8.h +++ b/include/xGooseFEM/ElementHex8.h @@ -1,132 +1,161 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_ELEMENTHEX8_H #define XGOOSEFEM_ELEMENTHEX8_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" -// ==================================== GooseFEM::Element::Hex8 ==================================== +// ================================================================================================= namespace xGooseFEM { namespace Element { namespace Hex8 { -// ======================================== tensor algebra ========================================= +// ================================================================================================= using T2 = xt::xtensor_fixed>; inline double inv(const T2 &A, T2 &Ainv); -// ================================ GooseFEM::Element::Hex8::Gauss ================================= +// ================================================================================================= 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 } -// ================================ GooseFEM::Element::Hex8::Nodal ================================= +// ================================================================================================= 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 } // ================================================================================================= // ------------------------------------------ quadrature ------------------------------------------- class Quadrature { -private: - - // dimensions (flexible) - size_t m_nelem; // number of elements - size_t m_nip; // number of integration points - - // dimensions (fixed for this element type) - static const size_t m_nne=8; // number of nodes per element - static const size_t m_ndim=3; // number of dimensions - - // data arrays - 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 gradients w.r.t. local coordinate [nip, nne, ndim] - xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] - xt::xtensor m_vol; // integration point volume [nelem, nip] - -private: - - // compute "vol" and "dNdx" based on current "x" - void compute_dN(); - public: - // convention: + // fixed dimensions: + // ndim = 3 - number of dimensions + // nne = 8 - 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) : Quadrature(x, Gauss::xi(), Gauss::w()){}; + 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 dimensions - size_t nelem() const; // number of elements - size_t nne() const; // number of nodes per element - size_t ndim() const; // number of dimension - size_t nip() const; // number of integration points + + size_t nelem() const { return m_nelem; }; // number of elements + size_t nne() const { return m_nne; }; // number of nodes per element + size_t ndim() const { return m_ndim; }; // number of dimension + size_t nip() const { return m_nip; }; // number of integration points // return integration volume - // - in-place + void dV(xt::xtensor &qscalar) const; - void dV(xt::xtensor &qtensor) const; - // - return qscalar/qtensor + + void dV(xt::xtensor &qtensor) const; // same volume for all tensor components + + // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part + + 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 dV() const; + xt::xtensor dVtensor() const; - // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part - // - in-place - 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; - // - return qtensor - xt::xtensor gradN_vector (const xt::xtensor &elemvec) const; - xt::xtensor gradN_vector_T (const xt::xtensor &elemvec) const; + 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; - // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV" - // - in-place - void int_N_scalar_NT_dV(const xt::xtensor &qscalar, xt::xtensor &elemmat) const; - // - return elemmat xt::xtensor int_N_scalar_NT_dV(const xt::xtensor &qscalar) const; - // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" - // - in-place - void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const; - // - return elemvec 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: + + // dimensions (flexible) + size_t m_nelem; // number of elements + size_t m_nip; // number of integration points + + // dimensions (fixed for this element type) + static const size_t m_nne=8; // number of nodes per element + static const size_t m_ndim=3; // number of dimensions + + // data arrays + 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 gradients w.r.t. local coordinate [nip, nne, ndim] + xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] + xt::xtensor m_vol; // integration point volume [nelem, nip] + }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementHex8.hpp b/include/xGooseFEM/ElementHex8.hpp index 1abddbc..f3b4e2a 100644 --- a/include/xGooseFEM/ElementHex8.hpp +++ b/include/xGooseFEM/ElementHex8.hpp @@ -1,714 +1,676 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_ELEMENTHEX8_CPP #define XGOOSEFEM_ELEMENTHEX8_CPP // ------------------------------------------------------------------------------------------------- #include "ElementHex8.h" -// ==================================== GooseFEM::Element::Hex8 ==================================== +// ================================================================================================= namespace xGooseFEM { namespace Element { namespace Hex8 { -// ======================================== tensor algebra ========================================= +// ================================================================================================= inline double inv(const T2 &A, T2 &Ainv) { // compute determinant - double det = ( A[0] * A[4] * A[8] + - A[1] * A[5] * A[6] + - A[2] * A[3] * A[7] ) - - ( A[2] * A[4] * A[6] + - A[1] * A[3] * A[8] + - A[0] * A[5] * A[7] ); + double det = ( A(0,0) * A(1,1) * A(2,2) + + A(0,1) * A(1,2) * A(2,0) + + A(0,2) * A(1,0) * A(2,1) ) - + ( A(0,2) * A(1,1) * A(2,0) + + A(0,1) * A(1,0) * A(2,2) + + A(0,0) * A(1,2) * A(2,1) ); // compute inverse - Ainv[0] = (A[4]*A[8]-A[5]*A[7]) / det; - Ainv[1] = (A[2]*A[7]-A[1]*A[8]) / det; - Ainv[2] = (A[1]*A[5]-A[2]*A[4]) / det; - Ainv[3] = (A[5]*A[6]-A[3]*A[8]) / det; - Ainv[4] = (A[0]*A[8]-A[2]*A[6]) / det; - Ainv[5] = (A[2]*A[3]-A[0]*A[5]) / det; - Ainv[6] = (A[3]*A[7]-A[4]*A[6]) / det; - Ainv[7] = (A[1]*A[6]-A[0]*A[7]) / det; - Ainv[8] = (A[0]*A[4]-A[1]*A[3]) / det; + Ainv(0,0) = (A(1,1)*A(2,2)-A(1,2)*A(2,1)) / det; + Ainv(0,1) = (A(0,2)*A(2,1)-A(0,1)*A(2,2)) / det; + Ainv(0,2) = (A(0,1)*A(1,2)-A(0,2)*A(1,1)) / det; + Ainv(1,0) = (A(1,2)*A(2,0)-A(1,0)*A(2,2)) / det; + Ainv(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0)) / det; + Ainv(1,2) = (A(0,2)*A(1,0)-A(0,0)*A(1,2)) / det; + Ainv(2,0) = (A(1,0)*A(2,1)-A(1,1)*A(2,0)) / det; + Ainv(2,1) = (A(0,1)*A(2,0)-A(0,0)*A(2,1)) / det; + Ainv(2,2) = (A(0,0)*A(1,1)-A(0,1)*A(1,0)) / det; return det; } -// ================================ GooseFEM::Element::Hex8::Gauss ================================= +// ================================================================================================= namespace Gauss { -// --------------------------------- number of integration points ---------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t nip() { return 8; } -// ----------------------- integration point coordinates (local coordinates) ----------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor xi() { size_t nip = 8; size_t ndim = 3; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); xi(0,2) = -1./std::sqrt(3.); xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); xi(1,2) = -1./std::sqrt(3.); xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); xi(2,2) = -1./std::sqrt(3.); xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); xi(3,2) = -1./std::sqrt(3.); xi(4,0) = -1./std::sqrt(3.); xi(4,1) = -1./std::sqrt(3.); xi(4,2) = +1./std::sqrt(3.); xi(5,0) = +1./std::sqrt(3.); xi(5,1) = -1./std::sqrt(3.); xi(5,2) = +1./std::sqrt(3.); xi(6,0) = +1./std::sqrt(3.); xi(6,1) = +1./std::sqrt(3.); xi(6,2) = +1./std::sqrt(3.); xi(7,0) = -1./std::sqrt(3.); xi(7,1) = +1./std::sqrt(3.); xi(7,2) = +1./std::sqrt(3.); return xi; } -// ----------------------------------- integration point weights ----------------------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor w() { size_t nip = 8; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; w(4) = 1.; w(5) = 1.; w(6) = 1.; w(7) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } -// ================================ GooseFEM::Element::Hex8::Nodal ================================ +// ================================================================================================= namespace Nodal { -// --------------------------------- number of integration points ---------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t nip() { return 8; } -// ----------------------- integration point coordinates (local coordinates) ----------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor xi() { size_t nip = 8; size_t ndim = 3; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1.; xi(0,1) = -1.; xi(0,2) = -1.; xi(1,0) = +1.; xi(1,1) = -1.; xi(1,2) = -1.; xi(2,0) = +1.; xi(2,1) = +1.; xi(2,2) = -1.; xi(3,0) = -1.; xi(3,1) = +1.; xi(3,2) = -1.; xi(4,0) = -1.; xi(4,1) = -1.; xi(4,2) = +1.; xi(5,0) = +1.; xi(5,1) = -1.; xi(5,2) = +1.; xi(6,0) = +1.; xi(6,1) = +1.; xi(6,2) = +1.; xi(7,0) = -1.; xi(7,1) = +1.; xi(7,2) = +1.; return xi; } -// ----------------------------------- integration point weights ----------------------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor w() { size_t nip = 8; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; w(4) = 1.; w(5) = 1.; w(6) = 1.; w(7) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= -// ------------------------------------------ constructor ------------------------------------------ - -inline Quadrature::Quadrature(const xt::xtensor &x) : m_x(x) -{ - assert( m_x.shape()[1] == m_nne ); - assert( m_x.shape()[2] == m_ndim ); - - // integration scheme - m_nip = Gauss::nip(); - m_xi = Gauss::xi(); - m_w = Gauss::w(); - - // extract number of elements - m_nelem = m_x.shape()[0]; - - // allocate arrays - m_N = xt::empty({ m_nip, m_nne }); - m_dNxi = xt::empty({ m_nip, m_nne, m_ndim}); - m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim}); - m_vol = xt::empty({m_nelem, m_nip }); - - // shape functions - for ( size_t k = 0 ; k < m_nip ; ++k ) - { - m_N(k,0) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,1) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,2) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,3) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,4) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); - m_N(k,5) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); - m_N(k,6) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); - m_N(k,7) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); - } - - // shape function gradients in local coordinates - for ( size_t k = 0 ; k < m_nip ; ++k ) - { - // - dN / dxi_0 - m_dNxi(k,0,0) = -.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,1,0) = +.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,2,0) = +.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,3,0) = -.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,4,0) = -.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); - m_dNxi(k,5,0) = +.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); - m_dNxi(k,6,0) = +.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); - m_dNxi(k,7,0) = -.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); - // - dN / dxi_1 - m_dNxi(k,0,1) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,1,1) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,2,1) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,3,1) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,4,1) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); - m_dNxi(k,5,1) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); - m_dNxi(k,6,1) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); - m_dNxi(k,7,1) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); - // - dN / dxi_2 - m_dNxi(k,0,2) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,1,2) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,2,2) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); - m_dNxi(k,3,2) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); - m_dNxi(k,4,2) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,5,2) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,6,2) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); - m_dNxi(k,7,2) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); - } - - // compute the shape function gradients, based on "x" - compute_dN(); -} - -// ------------------------------------------ constructor ------------------------------------------ - inline Quadrature::Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w) : m_x(x), m_w(w), m_xi(xi) { + // check input assert( m_x.shape()[1] == m_nne ); assert( m_x.shape()[2] == m_ndim ); // extract number of elements and number of integration points m_nelem = m_x.shape()[0]; m_nip = m_w.size(); + // check input assert( m_xi.shape()[0] == m_nip ); assert( m_xi.shape()[1] == m_ndim ); assert( m_w .size() == m_nip ); // allocate arrays m_N = xt::empty({ m_nip, m_nne }); m_dNxi = xt::empty({ m_nip, m_nne, m_ndim}); m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim}); m_vol = xt::empty({m_nelem, m_nip }); // shape functions - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { - m_N(k,0) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,1) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,2) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,3) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); - m_N(k,4) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); - m_N(k,5) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); - m_N(k,6) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); - m_N(k,7) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); + m_N(q,0) = .125 * (1.-m_xi(q,0)) * (1.-m_xi(q,1)) * (1.-m_xi(q,2)); + m_N(q,1) = .125 * (1.+m_xi(q,0)) * (1.-m_xi(q,1)) * (1.-m_xi(q,2)); + m_N(q,2) = .125 * (1.+m_xi(q,0)) * (1.+m_xi(q,1)) * (1.-m_xi(q,2)); + m_N(q,3) = .125 * (1.-m_xi(q,0)) * (1.+m_xi(q,1)) * (1.-m_xi(q,2)); + m_N(q,4) = .125 * (1.-m_xi(q,0)) * (1.-m_xi(q,1)) * (1.+m_xi(q,2)); + m_N(q,5) = .125 * (1.+m_xi(q,0)) * (1.-m_xi(q,1)) * (1.+m_xi(q,2)); + m_N(q,6) = .125 * (1.+m_xi(q,0)) * (1.+m_xi(q,1)) * (1.+m_xi(q,2)); + m_N(q,7) = .125 * (1.-m_xi(q,0)) * (1.+m_xi(q,1)) * (1.+m_xi(q,2)); } // shape function gradients in local coordinates - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - dN / dxi_0 - m_dNxi(k,0,0) = -.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,1,0) = +.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,2,0) = +.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,3,0) = -.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); - m_dNxi(k,4,0) = -.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); - m_dNxi(k,5,0) = +.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); - m_dNxi(k,6,0) = +.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); - m_dNxi(k,7,0) = -.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); + m_dNxi(q,0,0) = -.125*(1.-m_xi(q,1))*(1.-m_xi(q,2)); + m_dNxi(q,1,0) = +.125*(1.-m_xi(q,1))*(1.-m_xi(q,2)); + m_dNxi(q,2,0) = +.125*(1.+m_xi(q,1))*(1.-m_xi(q,2)); + m_dNxi(q,3,0) = -.125*(1.+m_xi(q,1))*(1.-m_xi(q,2)); + m_dNxi(q,4,0) = -.125*(1.-m_xi(q,1))*(1.+m_xi(q,2)); + m_dNxi(q,5,0) = +.125*(1.-m_xi(q,1))*(1.+m_xi(q,2)); + m_dNxi(q,6,0) = +.125*(1.+m_xi(q,1))*(1.+m_xi(q,2)); + m_dNxi(q,7,0) = -.125*(1.+m_xi(q,1))*(1.+m_xi(q,2)); // - dN / dxi_1 - m_dNxi(k,0,1) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,1,1) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,2,1) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,3,1) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); - m_dNxi(k,4,1) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); - m_dNxi(k,5,1) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); - m_dNxi(k,6,1) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); - m_dNxi(k,7,1) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); + m_dNxi(q,0,1) = -.125*(1.-m_xi(q,0))*(1.-m_xi(q,2)); + m_dNxi(q,1,1) = -.125*(1.+m_xi(q,0))*(1.-m_xi(q,2)); + m_dNxi(q,2,1) = +.125*(1.+m_xi(q,0))*(1.-m_xi(q,2)); + m_dNxi(q,3,1) = +.125*(1.-m_xi(q,0))*(1.-m_xi(q,2)); + m_dNxi(q,4,1) = -.125*(1.-m_xi(q,0))*(1.+m_xi(q,2)); + m_dNxi(q,5,1) = -.125*(1.+m_xi(q,0))*(1.+m_xi(q,2)); + m_dNxi(q,6,1) = +.125*(1.+m_xi(q,0))*(1.+m_xi(q,2)); + m_dNxi(q,7,1) = +.125*(1.-m_xi(q,0))*(1.+m_xi(q,2)); // - dN / dxi_2 - m_dNxi(k,0,2) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,1,2) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,2,2) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); - m_dNxi(k,3,2) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); - m_dNxi(k,4,2) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,5,2) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); - m_dNxi(k,6,2) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); - m_dNxi(k,7,2) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); + m_dNxi(q,0,2) = -.125*(1.-m_xi(q,0))*(1.-m_xi(q,1)); + m_dNxi(q,1,2) = -.125*(1.+m_xi(q,0))*(1.-m_xi(q,1)); + m_dNxi(q,2,2) = -.125*(1.+m_xi(q,0))*(1.+m_xi(q,1)); + m_dNxi(q,3,2) = -.125*(1.-m_xi(q,0))*(1.+m_xi(q,1)); + m_dNxi(q,4,2) = +.125*(1.-m_xi(q,0))*(1.-m_xi(q,1)); + m_dNxi(q,5,2) = +.125*(1.+m_xi(q,0))*(1.-m_xi(q,1)); + m_dNxi(q,6,2) = +.125*(1.+m_xi(q,0))*(1.+m_xi(q,1)); + m_dNxi(q,7,2) = +.125*(1.-m_xi(q,0))*(1.+m_xi(q,1)); } // compute the shape function gradients, based on "x" compute_dN(); } -// --------------------------- integration volume (per tensor-component) --------------------------- +// ------------------------------------------------------------------------------------------------- inline void Quadrature::dV(xt::xtensor &qscalar) const { assert( qscalar.shape()[0] == m_nelem ); assert( qscalar.shape()[1] == m_nip ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t k = 0 ; k < m_nip ; ++k ) - qscalar(e,k) = m_vol(e,k); + for ( size_t q = 0 ; q < m_nip ; ++q ) + qscalar(e,q) = m_vol(e,q); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::dV(xt::xtensor &qtensor) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) - qtensor(e,k,i,j) = m_vol(e,k); + qtensor(e,q,i,j) = m_vol(e,q); } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor Quadrature::dV() const -{ - xt::xtensor out = xt::empty({m_nelem, m_nip}); - - this->dV(out); - - return out; -} - -// ------------------------------------------------------------------------------------------------- - -inline xt::xtensor Quadrature::dVtensor() const -{ - xt::xtensor out = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); - - this->dV(out); - - return out; -} - -// -------------------------------------- number of elements --------------------------------------- - -inline size_t Quadrature::nelem() const -{ - return m_nelem; -} - -// ---------------------------------- number of nodes per element ---------------------------------- - -inline size_t Quadrature::nne() const -{ - return m_nne; -} - -// ------------------------------------- number of dimensions -------------------------------------- - -inline size_t Quadrature::ndim() const -{ - return m_ndim; -} - -// --------------------------------- number of integration points ---------------------------------- - -inline size_t Quadrature::nip() const -{ - return m_nip; -} - -// --------------------------------------- update positions ---------------------------------------- - inline void Quadrature::update_x(const xt::xtensor &x) { assert( x.shape()[0] == m_nelem ); assert( x.shape()[1] == m_nne ); assert( x.shape()[2] == m_ndim ); assert( x.size() == m_x.size() ); // update positions xt::noalias(m_x) = x; // update the shape function gradients for the new "x" compute_dN(); } -// ------------------------ shape function gradients in global coordinates ------------------------- +// ------------------------------------------------------------------------------------------------- inline void Quadrature::compute_dN() { // loop over all elements (in parallel) #pragma omp parallel { // allocate local variables T2 J, Jinv; #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias nodal positions auto x = xt::adapt(&m_x(e,0,0), xt::xshape()); // loop over integration points - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNxi = xt::adapt(&m_dNxi( k,0,0), xt::xshape()); - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); + auto dNxi = xt::adapt(&m_dNxi( q,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); // - zero-initialize J.fill(0.0); // - Jacobian for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) J(i,j) += dNxi(m,i) * x(m,j); // - determinant and inverse of the Jacobian double Jdet = inv(J, Jinv); // - shape function gradients wrt global coordinates (loops partly unrolled for efficiency) // dNx(m,i) += Jinv(i,j) * dNxi(m,j); for ( size_t m = 0 ; m < m_nne ; ++m ) { dNx(m,0) = Jinv(0,0) * dNxi(m,0) + Jinv(0,1) * dNxi(m,1) + Jinv(0,2) * dNxi(m,2); dNx(m,1) = Jinv(1,0) * dNxi(m,0) + Jinv(1,1) * dNxi(m,1) + Jinv(1,2) * dNxi(m,2); dNx(m,2) = Jinv(2,0) * dNxi(m,0) + Jinv(2,1) * dNxi(m,1) + Jinv(2,2) * dNxi(m,2); } // - integration point volume - m_vol(e,k) = m_w(k) * Jdet; + m_vol(e,q) = m_w(q) * Jdet; } } } } -// ------------------- dyadic product "qtensor(i,j) = dNdx(m,i) * elemvec(m,j)" -------------------- +// ------------------------------------------------------------------------------------------------- inline void Quadrature::gradN_vector( const xt::xtensor &elemvec, xt::xtensor &qtensor) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); // zero-initialize output: matrix of tensors qtensor.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias element vector (e.g. nodal displacements) auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto gradu = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); // - evaluate dyadic product for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) gradu(i,j) += dNx(m,i) * u(m,j); } } } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const -{ - xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); - - this->gradN_vector(elemvec, qtensor); - - return qtensor; -} - -// ---------------------------------- transpose of "gradN_vector" ---------------------------------- - inline void Quadrature::gradN_vector_T( const xt::xtensor &elemvec, xt::xtensor &qtensor) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); // zero-initialize output: matrix of tensors qtensor.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias element vector (e.g. nodal displacements) auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto gradu = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); // - evaluate transpose of dyadic product for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) gradu(j,i) += dNx(m,i) * u(m,j); } } } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const -{ - xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); - - this->gradN_vector_T(elemvec, qtensor); - - return qtensor; -} - -// ------------------------------- symmetric part of "gradN_vector" -------------------------------- - inline void Quadrature::symGradN_vector( const xt::xtensor &elemvec, xt::xtensor &qtensor) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); // zero-initialize output: matrix of tensors qtensor.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias element vector (e.g. nodal displacements) auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto eps = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto eps = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); // - evaluate symmetrized dyadic product for ( size_t m = 0 ; m < m_nne ; ++m ) { for ( size_t i = 0 ; i < m_ndim ; ++i ) { for ( size_t j = 0 ; j < m_ndim ; ++j ) { eps(i,j) += dNx(m,i) * u(m,j) / 2.; eps(j,i) += dNx(m,i) * u(m,j) / 2.; } } } } } } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const -{ - xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); - - this->symGradN_vector(elemvec, qtensor); - - return qtensor; -} - -// ------- scalar product "elemmat(m*ndim+i,n*ndim+i) = N(m) * qscalar * N(n)"; for all "i" -------- - inline void Quadrature::int_N_scalar_NT_dV( const xt::xtensor &qscalar, xt::xtensor &elemmat) const { assert( qscalar.shape()[0] == m_nelem ); assert( qscalar.shape()[1] == m_nip ); assert( elemmat.shape()[0] == m_nelem ); assert( elemmat.shape()[1] == m_nne*m_ndim ); assert( elemmat.shape()[2] == m_nne*m_ndim ); // zero-initialize: matrix of matrices elemmat.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias (e.g. mass matrix) auto M = xt::adapt(&elemmat(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto N = xt::adapt(&m_N(k,0), xt::xshape()); - auto& vol = m_vol (e,k); - auto& rho = qscalar(e,k); + auto N = xt::adapt(&m_N(q,0), xt::xshape()); + auto& vol = m_vol (e,q); + auto& rho = qscalar(e,q); // - evaluate scalar product, for all dimensions, and assemble // 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; M(m*m_ndim+2, n*m_ndim+2) += N(m) * rho * N(n) * vol; } } } } } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor Quadrature::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; -} - -// ------------ integral of dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" ------------ - inline void Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nip ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); // zero-initialize output: matrix of vectors elemvec.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias (e.g. nodal force) auto f = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto sig = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); - auto& vol = m_vol(e,k); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto sig = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); + auto& vol = m_vol(e,q); // - evaluate dot product, and assemble for ( size_t m = 0 ; m < m_nne ; ++m ) { f(m,0) += ( dNx(m,0) * sig(0,0) + dNx(m,1) * sig(1,0) + dNx(m,2) * sig(2,0) ) * vol; f(m,1) += ( dNx(m,0) * sig(0,1) + dNx(m,1) * sig(1,1) + dNx(m,2) * sig(2,1) ) * vol; f(m,2) += ( dNx(m,0) * sig(0,2) + dNx(m,1) * sig(1,2) + dNx(m,2) * sig(2,2) ) * vol; } } } } // ------------------------------------------------------------------------------------------------- -inline xt::xtensor Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const +inline void Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor, + xt::xtensor &elemmat) const +{ + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nip ); + assert( qtensor.shape()[2] == m_ndim ); + assert( qtensor.shape()[3] == m_ndim ); + assert( qtensor.shape()[4] == m_ndim ); + assert( qtensor.shape()[5] == m_ndim ); + + assert( elemmat.shape()[0] == m_nelem ); + assert( elemmat.shape()[1] == m_nne*m_ndim ); + assert( elemmat.shape()[2] == m_nne*m_ndim ); + + // zero-initialize output: matrix of vector + elemmat.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias (e.g. nodal force) + auto K = xt::adapt(&elemmat(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ){ + + // - alias + auto dNx = xt::adapt(&m_dNx(e,q,0,0), xt::xshape()); + auto C = xt::adapt(&qtensor(e,q,0,0,0,0), xt::xshape()); + auto& vol = m_vol(e,q); + + // - evaluate dot product, and assemble + for ( size_t m = 0 ; m < m_nne ; ++m ) + for ( size_t n = 0 ; n < m_nne ; ++n ) + for ( size_t i = 0 ; i < m_ndim ; ++i ) + for ( size_t j = 0 ; j < m_ndim ; ++j ) + for ( size_t k = 0 ; k < m_ndim ; ++k ) + for ( size_t l = 0 ; l < m_ndim ; ++l ) + K(m*m_ndim+j, n*m_ndim+k) += dNx(m,i) * C(i,j,k,l) * dNx(n,l) * vol; + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor Quadrature::dV() const +{ + xt::xtensor out = xt::empty({m_nelem, m_nip}); + + this->dV(out); + + return out; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor Quadrature::dVtensor() const +{ + xt::xtensor out = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); + + this->dV(out); + + return out; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const +{ + xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); + + this->gradN_vector(elemvec, qtensor); + + return qtensor; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const +{ + xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); + + this->gradN_vector_T(elemvec, qtensor); + + return qtensor; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const +{ + xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); + + this->symGradN_vector(elemvec, qtensor); + + return qtensor; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor Quadrature::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 Quadrature::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 Quadrature::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 ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementQuad4.h b/include/xGooseFEM/ElementQuad4.h index 1bac6d1..4311538 100644 --- a/include/xGooseFEM/ElementQuad4.h +++ b/include/xGooseFEM/ElementQuad4.h @@ -1,156 +1,161 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_ELEMENTQUAD4_H #define XGOOSEFEM_ELEMENTQUAD4_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" -// =================================== GooseFEM::Element::Quad4 ==================================== +// ================================================================================================= namespace xGooseFEM { namespace Element { namespace Quad4 { -// ======================================== tensor algebra ========================================= +// ================================================================================================= using T2 = xt::xtensor_fixed>; inline double inv(const T2 &A, T2 &Ainv); -// ================================ GooseFEM::Element::Quad4::Gauss ================================ +// ================================================================================================= 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 } -// ================================ GooseFEM::Element::Quad4::Nodal ================================ +// ================================================================================================= 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 } // ================================================================================================= // ------------------------------------------ quadrature ------------------------------------------- class Quadrature { -private: - - // dimensions (flexible) - size_t m_nelem; // number of elements - size_t m_nip; // number of integration points - - // dimensions (fixed for this element type) - static const size_t m_nne=4; // number of nodes per element - static const size_t m_ndim=2; // number of dimensions - - // data arrays - 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 gradients w.r.t. local coordinate [nip, nne, ndim] - xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] - xt::xtensor m_vol; // integration point volume [nelem, nip] - -private: - - // compute "vol" and "dNdx" based on current "x" - void compute_dN(); - public: - // convention: + // 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) : Quadrature(x, Gauss::xi(), Gauss::w()){}; 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 dimensions - size_t nelem() const; // number of elements - size_t nne() const; // number of nodes per element - size_t ndim() const; // number of dimension - size_t nip() const; // number of integration points + size_t nelem() const { return m_nelem; }; // number of elements + size_t nne() const { return m_nne; }; // number of nodes per element + size_t ndim() const { return m_ndim; }; // number of dimension + size_t nip() const { return m_nip; }; // number of integration points // return integration volume void dV(xt::xtensor &qscalar) const; + void dV(xt::xtensor &qtensor) const; // same volume for all tensor components // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part - void gradN_vector (const xt::xtensor &elemvec, + void gradN_vector(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; - void gradN_vector_T (const xt::xtensor &elemvec, + 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 dV() const; xt::xtensor dVtensor() const; 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: + + // dimensions (flexible) + size_t m_nelem; // number of elements + size_t m_nip; // number of integration points + + // dimensions (fixed for this element type) + static const size_t m_nne=4; // number of nodes per element + static const size_t m_ndim=2; // number of dimensions + + // data arrays + 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 gradients w.r.t. local coordinate [nip, nne, ndim] + xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] + xt::xtensor m_vol; // integration point volume [nelem, nip] + }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementQuad4.hpp b/include/xGooseFEM/ElementQuad4.hpp index 1550a14..2e8637a 100644 --- a/include/xGooseFEM/ElementQuad4.hpp +++ b/include/xGooseFEM/ElementQuad4.hpp @@ -1,702 +1,625 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_ELEMENTQUAD4_CPP #define XGOOSEFEM_ELEMENTQUAD4_CPP // ------------------------------------------------------------------------------------------------- #include "ElementQuad4.h" // ================================================================================================= namespace xGooseFEM { namespace Element { namespace Quad4 { // ================================================================================================= inline double inv(const T2 &A, T2 &Ainv) { // compute determinant - double det = A[0] * A[3] - A[1] * A[2]; + double det = A(0,0) * A(1,1) - A(0,1) * A(1,0); // compute inverse - Ainv[0] = A[3] / det; - Ainv[1] = -1. * A[1] / det; - Ainv[2] = -1. * A[2] / det; - Ainv[3] = A[0] / det; + Ainv(0,0) = A(1,1) / det; + Ainv(0,1) = -1. * A(0,1) / det; + Ainv(1,0) = -1. * A(1,0) / det; + Ainv(1,1) = A(0,0) / det; return det; } // ================================================================================================= namespace Gauss { // ------------------------------------------------------------------------------------------------- inline size_t nip() { return 4; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor xi() { size_t nip = 4; size_t ndim = 2; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); return xi; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor w() { size_t nip = 4; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= namespace Nodal { // ------------------------------------------------------------------------------------------------- inline size_t nip() { return 4; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor xi() { size_t nip = 4; size_t ndim = 2; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1.; xi(0,1) = -1.; xi(1,0) = +1.; xi(1,1) = -1.; xi(2,0) = +1.; xi(2,1) = +1.; xi(3,0) = -1.; xi(3,1) = +1.; return xi; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor w() { size_t nip = 4; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= -// ------------------------------------------------------------------------------------------------- - -inline Quadrature::Quadrature(const xt::xtensor &x) : m_x(x) -{ - assert( m_x.shape()[1] == m_nne ); - assert( m_x.shape()[2] == m_ndim ); - - // integration scheme - m_nip = Gauss::nip(); - m_xi = Gauss::xi(); - m_w = Gauss::w(); - - // extract number of elements - m_nelem = m_x.shape()[0]; - - // allocate arrays - m_N = xt::empty({ m_nip, m_nne }); - m_dNxi = xt::empty({ m_nip, m_nne, m_ndim}); - m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim}); - m_vol = xt::empty({m_nelem, m_nip }); - - // shape functions - for ( size_t k = 0 ; k < m_nip ; ++k ) - { - m_N(k,0) = .25 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)); - m_N(k,1) = .25 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)); - m_N(k,2) = .25 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)); - m_N(k,3) = .25 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)); - } - - // shape function gradients in local coordinates - for ( size_t k = 0 ; k < m_nip ; ++k ) - { - // - dN / dxi_0 - m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1)); - m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1)); - m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1)); - m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1)); - // - dN / dxi_1 - m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0)); - m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0)); - m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0)); - m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0)); - } - - // compute the shape function gradients, based on "x" - compute_dN(); -} - -// ------------------------------------------------------------------------------------------------- - inline Quadrature::Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w) : m_x(x), m_w(w), m_xi(xi) { + // check input assert( m_x.shape()[1] == m_nne ); assert( m_x.shape()[2] == m_ndim ); // extract number of elements and number of integration points m_nelem = m_x.shape()[0]; m_nip = m_w.size(); + // check input assert( m_xi.shape()[0] == m_nip ); assert( m_xi.shape()[1] == m_ndim ); assert( m_w .size() == m_nip ); // allocate arrays m_N = xt::empty({ m_nip, m_nne }); m_dNxi = xt::empty({ m_nip, m_nne, m_ndim}); m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim}); m_vol = xt::empty({m_nelem, m_nip }); // shape functions - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { - m_N(k,0) = .25 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)); - m_N(k,1) = .25 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)); - m_N(k,2) = .25 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)); - m_N(k,3) = .25 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)); + m_N(q,0) = .25 * (1.-m_xi(q,0)) * (1.-m_xi(q,1)); + m_N(q,1) = .25 * (1.+m_xi(q,0)) * (1.-m_xi(q,1)); + m_N(q,2) = .25 * (1.+m_xi(q,0)) * (1.+m_xi(q,1)); + m_N(q,3) = .25 * (1.-m_xi(q,0)) * (1.+m_xi(q,1)); } // shape function gradients in local coordinates - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - dN / dxi_0 - m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1)); - m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1)); - m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1)); - m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1)); + m_dNxi(q,0,0) = -.25*(1.-m_xi(q,1)); + m_dNxi(q,1,0) = +.25*(1.-m_xi(q,1)); + m_dNxi(q,2,0) = +.25*(1.+m_xi(q,1)); + m_dNxi(q,3,0) = -.25*(1.+m_xi(q,1)); // - dN / dxi_1 - m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0)); - m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0)); - m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0)); - m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0)); + m_dNxi(q,0,1) = -.25*(1.-m_xi(q,0)); + m_dNxi(q,1,1) = -.25*(1.+m_xi(q,0)); + m_dNxi(q,2,1) = +.25*(1.+m_xi(q,0)); + m_dNxi(q,3,1) = +.25*(1.-m_xi(q,0)); } // compute the shape function gradients, based on "x" compute_dN(); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::dV(xt::xtensor &qscalar) const { assert( qscalar.shape()[0] == m_nelem ); assert( qscalar.shape()[1] == m_nip ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t k = 0 ; k < m_nip ; ++k ) - qscalar(e,k) = m_vol(e,k); + for ( size_t q = 0 ; q < m_nip ; ++q ) + qscalar(e,q) = m_vol(e,q); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::dV(xt::xtensor &qtensor) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) - qtensor(e,k,i,j) = m_vol(e,k); -} - -// ------------------------------------------------------------------------------------------------- - -inline size_t Quadrature::nelem() const -{ - return m_nelem; -} - -// ------------------------------------------------------------------------------------------------- - -inline size_t Quadrature::nne() const -{ - return m_nne; -} - -// ------------------------------------------------------------------------------------------------- - -inline size_t Quadrature::ndim() const -{ - return m_ndim; -} - -// ------------------------------------------------------------------------------------------------- - -inline size_t Quadrature::nip() const -{ - return m_nip; + qtensor(e,q,i,j) = m_vol(e,q); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::update_x(const xt::xtensor &x) { assert( x.shape()[0] == m_nelem ); assert( x.shape()[1] == m_nne ); assert( x.shape()[2] == m_ndim ); assert( x.size() == m_x.size() ); // update positions xt::noalias(m_x) = x; // update the shape function gradients for the new "x" compute_dN(); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::compute_dN() { // loop over all elements (in parallel) #pragma omp parallel { // allocate local variables T2 J, Jinv; #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias nodal positions auto x = xt::adapt(&m_x(e,0,0), xt::xshape()); // loop over integration points - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNxi = xt::adapt(&m_dNxi( k,0,0), xt::xshape()); - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); + auto dNxi = xt::adapt(&m_dNxi( q,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); // - Jacobian (loops unrolled for efficiency) // 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); // - determinant and inverse of the Jacobian double Jdet = inv(J, Jinv); // - shape function gradients wrt global coordinates (loops partly unrolled for efficiency) // dNx(m,i) += Jinv(i,j) * dNxi(m,j); for ( size_t m = 0 ; m < m_nne ; ++m ) { dNx(m,0) = Jinv(0,0) * dNxi(m,0) + Jinv(0,1) * dNxi(m,1); dNx(m,1) = Jinv(1,0) * dNxi(m,0) + Jinv(1,1) * dNxi(m,1); } // - integration point volume - m_vol(e,k) = m_w(k) * Jdet; + m_vol(e,q) = m_w(q) * Jdet; } } } } // ------------------------------------------------------------------------------------------------- inline void Quadrature::gradN_vector( const xt::xtensor &elemvec, xt::xtensor &qtensor) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); // zero-initialize output: matrix of tensors qtensor.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias element vector (e.g. nodal displacements) auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto gradu = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); // - evaluate dyadic product (loops unrolled for efficiency) // gradu(i,j) += dNx(m,i) * u(m,j) gradu(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0); gradu(0,1) = dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1); gradu(1,0) = dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0); gradu(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1); } } } // ------------------------------------------------------------------------------------------------- inline void Quadrature::gradN_vector_T( const xt::xtensor &elemvec, xt::xtensor &qtensor) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); // zero-initialize output: matrix of tensors qtensor.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias element vector (e.g. nodal displacements) auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto gradu = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); // - evaluate transpose of dyadic product (loops unrolled for efficiency) // gradu(j,i) += dNx(m,i) * u(m,j) gradu(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0); gradu(1,0) = dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1); gradu(0,1) = dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0); gradu(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1); } } } // ------------------------------------------------------------------------------------------------- inline void Quadrature::symGradN_vector( const xt::xtensor &elemvec, xt::xtensor &qtensor) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); // zero-initialize output: matrix of tensors qtensor.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias element vector (e.g. nodal displacements) auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto eps = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto eps = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); // - evaluate symmetrized dyadic product (loops unrolled for efficiency) // grad(i,j) += dNx(m,i) * u(m,j) // eps (j,i) = 0.5 * ( grad(i,j) + grad(j,i) ) eps(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0); eps(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1); eps(0,1) = ( dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1) + - dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0) ) / 2.; + dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0) ) * 0.5; eps(1,0) = eps(0,1); } } } -// ------- scalar product "elemmat(m*ndim+i,n*ndim+i) = N(m) * qscalar * N(n)"; for all "i" -------- +// ------------------------------------------------------------------------------------------------- inline void Quadrature::int_N_scalar_NT_dV( const xt::xtensor &qscalar, xt::xtensor &elemmat) const { assert( qscalar.shape()[0] == m_nelem ); assert( qscalar.shape()[1] == m_nip ); assert( elemmat.shape()[0] == m_nelem ); assert( elemmat.shape()[1] == m_nne*m_ndim ); assert( elemmat.shape()[2] == m_nne*m_ndim ); // zero-initialize: matrix of matrices elemmat.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias (e.g. mass matrix) auto M = xt::adapt(&elemmat(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto N = xt::adapt(&m_N(k,0), xt::xshape()); - auto& vol = m_vol (e,k); - auto& rho = qscalar(e,k); + auto N = xt::adapt(&m_N(q,0), xt::xshape()); + auto& vol = m_vol (e,q); + auto& rho = qscalar(e,q); // - evaluate scalar product, for all dimensions, and assemble // 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; } } } } } -// ------------ integral of dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" ------------ +// ------------------------------------------------------------------------------------------------- inline void Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nip ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); // zero-initialize output: matrix of vectors elemvec.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias (e.g. nodal force) auto f = xt::adapt(&elemvec(e,0,0), xt::xshape()); // loop over all integration points in element "e" - for ( size_t k = 0 ; k < m_nip ; ++k ) + for ( size_t q = 0 ; q < m_nip ; ++q ) { // - alias - auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); - auto sig = xt::adapt(&qtensor(e,k,0,0), xt::xshape()); - auto& vol = m_vol(e,k); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto sig = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); + auto& vol = m_vol(e,q); // - evaluate dot product, and assemble for ( size_t m = 0 ; m < m_nne ; ++m ) { f(m,0) += ( dNx(m,0) * sig(0,0) + dNx(m,1) * sig(1,0) ) * vol; f(m,1) += ( dNx(m,0) * sig(0,1) + dNx(m,1) * sig(1,1) ) * vol; } } } } // ------------------------------------------------------------------------------------------------- inline void Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor, xt::xtensor &elemmat) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nip ); assert( qtensor.shape()[2] == m_ndim ); assert( qtensor.shape()[3] == m_ndim ); assert( qtensor.shape()[4] == m_ndim ); assert( qtensor.shape()[5] == m_ndim ); assert( elemmat.shape()[0] == m_nelem ); assert( elemmat.shape()[1] == m_nne*m_ndim ); assert( elemmat.shape()[2] == m_nne*m_ndim ); // zero-initialize output: matrix of vector elemmat.fill(0.0); // loop over all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias (e.g. nodal force) auto K = xt::adapt(&elemmat(e,0,0), xt::xshape()); // loop over all integration points in element "e" for ( size_t q = 0 ; q < m_nip ; ++q ){ // - alias auto dNx = xt::adapt(&m_dNx(e,q,0,0), xt::xshape()); auto C = xt::adapt(&qtensor(e,q,0,0,0,0), xt::xshape()); auto& vol = m_vol(e,q); // - evaluate dot product, and assemble for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t n = 0 ; n < m_nne ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) for ( size_t k = 0 ; k < m_ndim ; ++k ) for ( size_t l = 0 ; l < m_ndim ; ++l ) K(m*m_ndim+j, n*m_ndim+k) += dNx(m,i) * C(i,j,k,l) * dNx(n,l) * vol; } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::dV() const { xt::xtensor out = xt::empty({m_nelem, m_nip}); this->dV(out); return out; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::dVtensor() const { xt::xtensor out = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->dV(out); return out; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->gradN_vector(elemvec, qtensor); return qtensor; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->gradN_vector_T(elemvec, qtensor); return qtensor; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->symGradN_vector(elemvec, qtensor); return qtensor; } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::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 Quadrature::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 Quadrature::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 ... // ================================================================================================= #endif