diff --git a/src/GooseFEM/ElementQuad4.cpp b/src/GooseFEM/ElementQuad4.cpp index ede3d7b..67aa9b8 100644 --- a/src/GooseFEM/ElementQuad4.cpp +++ b/src/GooseFEM/ElementQuad4.cpp @@ -1,644 +1,644 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_ELEMENTQUAD4_CPP #define GOOSEFEM_ELEMENTQUAD4_CPP // ------------------------------------------------------------------------------------------------- #include "ElementQuad4.h" // =================================== GooseFEM::Element::Quad4 ==================================== namespace GooseFEM { namespace Element { namespace Quad4 { // ================================ GooseFEM::Element::Quad4::Gauss ================================ namespace Gauss { // --------------------------------- number of integration points ---------------------------------- inline size_t nip() { return 4; } // ----------------------- integration point coordinates (local coordinates) ----------------------- inline ArrD xi() { size_t nip = 4; size_t ndim = 2; ArrD xi({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; } // ----------------------------------- integration point weights ----------------------------------- inline ArrD w() { size_t nip = 4; ArrD w({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================ GooseFEM::Element::Quad4::Nodal ================================ namespace Nodal { // --------------------------------- number of integration points ---------------------------------- inline size_t nip() { return 4; } // ----------------------- integration point coordinates (local coordinates) ----------------------- inline ArrD xi() { size_t nip = 4; size_t ndim = 2; ArrD xi({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; } // ----------------------------------- integration point weights ----------------------------------- inline ArrD w() { size_t nip = 4; ArrD w({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= // ------------------------------------------ constructor ------------------------------------------ inline Quadrature::Quadrature(const ArrD &x, const ArrD &xi, const ArrD &w) : m_x(x), m_w(w), m_xi(xi) { // check input assert( m_x.ndim() == 3 ); // shape: [nelem, nne, ndim] assert( m_x.shape(1) == m_nne ); // number of nodes per element assert( m_x.shape(2) == m_ndim ); // number of dimensions // extract number of elements m_nelem = m_x.shape(0); // integration scheme // - default if ( m_w.size() == 0 and m_xi.size() == 0 ) { m_nip = Gauss::nip(); m_xi = Gauss::xi(); m_w = Gauss::w(); } // - input else if ( m_w.size() > 0 and m_xi.size() > 0 ) { m_nip = m_w.size(); } // - unknown else { throw std::runtime_error("Input integration point coordinates and weights"); } // check input assert( m_xi.ndim() == 2 ); // shape: [nip, ndim] assert( m_xi.shape(0) == m_nip ); // number of integration points assert( m_xi.shape(1) == m_ndim ); // number of dimensions assert( m_w .ndim() == 1 ); // shape: [nip] assert( m_w .size() == m_nip ); // number of integration points // allocate arrays // - shape functions m_N.resize({m_nip,m_nne}); // - shape function gradients in local coordinates m_dNxi.resize({m_nip,m_nne,m_ndim}); // - shape function gradients in global coordinates m_dNx.resize({m_nelem,m_nip,m_nne,m_ndim}); // - integration point volume m_vol.resize({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 ) { m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1)); m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0)); m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1)); m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0)); m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1)); m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0)); m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1)); m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0)); } // compute the shape function gradients, based on "x" compute_dN(); } // -------------------------------------- integration volume --------------------------------------- inline ArrD Quadrature::dV() const { return m_vol; } // ---------------------------- integration volume per tensor-component ---------------------------- inline ArrD Quadrature::dV(size_t ncomp) const { ArrD out = ArrD::Zero({m_nelem, m_nip, ncomp}); #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 i = 0 ; i < ncomp ; ++i ) out(e,k,i) = m_vol(e,k); 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 ArrD &x) { // check input assert( x.ndim() == 3 ); // shape: [nelem, nne, ndim] assert( x.shape(0) == m_nelem ); // number of elements assert( x.shape(1) == m_nne ); // number of nodes per element assert( x.shape(2) == m_ndim ); // number of dimensions assert( x.size() == m_x.size() ); // total number of components (redundant) // update positions m_x.setCopy(x.begin(), x.end()); // update the shape function gradients for the new "x" compute_dN(); } // ------------------------ shape function gradients in global coordinates ------------------------- inline void Quadrature::compute_dN() { #pragma omp parallel { // intermediate quantities and local views double Jdet, w, vol; - cppmat::tiny::matrix dNx; - cppmat::view::matrix dNxi, x; + cppmat::tiny::matrix dNx; + cppmat::view::matrix dNxi, x; cppmat::cartesian2d::tensor2 J, Jinv; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias x.setMap(&m_x(e)); // nodal positions // loop over integration points for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias dNxi.setMap(&m_dNxi( k)); // shape function gradients (local coordinates) vol = m_vol (e,k); // volume w = m_w ( k); // weight factor // - Jacobian (loops unrolled for efficiency) // J(i,j) += dNxi(m,i) * xe(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 Jdet = J.det(); Jinv = J.inv(); // - integration point volume vol = w * Jdet; // - 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); } // - copy to matrix: shape function gradients (global coordinates) std::copy(dNx.begin(), dNx.end(), m_dNx.item(e,k)); // - copy to matrix: integration point volume m_vol(e,k) = vol; } } } // #pragma omp parallel } // ------------------- dyadic product "qtensor(i,j) = dNdx(m,i) * elemvec(m,j)" -------------------- template inline ArrD Quadrature::gradN_vector(const ArrD &elemvec) const { // check input assert( elemvec.ndim() == 3 ); // shape: [nelem, nne, ndim] assert( elemvec.shape(0) == m_nelem ); // number of elements assert( elemvec.shape(1) == m_nne ); // number of nodes per element assert( elemvec.shape(2) == m_ndim ); // number of dimensions // temporary tensor, to deduce size of the output T dummytensor; size_t tdim = dummytensor.size(); // zero-initialize output: matrix of tensors ArrD qtensor = ArrD::Zero({m_nelem, m_nip, tdim}); #pragma omp parallel { // intermediate quantities and local views T gradu; - cppmat::view::matrix dNx, u; + cppmat::view::matrix dNx, u; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias u.setMap(&elemvec(e)); // element vector (e.g. nodal displacements) // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias dNx.setMap(&m_dNx(e,k)); // shape function gradients (global coordinates) // - evaluate dyadic product (loops unrolled for efficiency) // gradu(i,j) += dNx(m,i) * ue(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); // - copy resulting integration point tensor std::copy(gradu.begin(), gradu.end(), qtensor.item(e,k)); } } } // #pragma omp parallel return qtensor; } // ---------------------------------- transpose of "GradN_vector" ---------------------------------- template inline ArrD Quadrature::gradN_vector_T(const ArrD &elemvec) const { // check input assert( elemvec.ndim() == 3 ); // shape: [nelem, nne, ndim] assert( elemvec.shape(0) == m_nelem ); // number of elements assert( elemvec.shape(1) == m_nne ); // number of nodes per element assert( elemvec.shape(2) == m_ndim ); // number of dimensions // temporary tensor, to deduce size of the output T dummytensor; size_t tdim = dummytensor.size(); // zero-initialize output: matrix of tensors ArrD qtensor = ArrD::Zero({m_nelem, m_nip, tdim}); #pragma omp parallel { // intermediate quantities and local views T gradu; - cppmat::view::matrix dNx, u; + cppmat::view::matrix dNx, u; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias u.setMap(&elemvec(e)); // element vector (e.g. nodal displacements) // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias dNx.setMap(&m_dNx(e,k)); // shape function gradients (global coordinates) // - evaluate dyadic product (loops unrolled for efficiency) // gradu(j,i) += dNx(m,i) * ue(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); // - copy resulting integration point tensor std::copy(gradu.begin(), gradu.end(), qtensor.item(e,k)); } } } // #pragma omp parallel return qtensor; } // ------------------------------- symmetric part of "GradN_vector" -------------------------------- template inline ArrD Quadrature::symGradN_vector(const ArrD &elemvec) const { // check input assert( elemvec.ndim() == 3 ); // shape: [nelem, nne, ndim] assert( elemvec.shape(0) == m_nelem ); // number of elements assert( elemvec.shape(1) == m_nne ); // number of nodes per element assert( elemvec.shape(2) == m_ndim ); // number of dimensions // temporary tensor, to deduce size of the output T dummytensor; size_t tdim = dummytensor.size(); // zero-initialize output: matrix of tensors ArrD qtensor = ArrD::Zero({m_nelem, m_nip, tdim}); #pragma omp parallel { // intermediate quantities and local views T eps; cppmat::cartesian2d::tensor2 gradu; - cppmat::view::matrix dNx, u; + cppmat::view::matrix dNx, u; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias u.setMap(&elemvec(e)); // element vector (e.g. nodal displacements) // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias dNx.setMap(&m_dNx(e,k)); // shape function gradients (global coordinates) // - evaluate dyadic product (loops unrolled for efficiency) // gradu(i,j) += dNx(m,i) * ue(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); // - symmetrize (loops unrolled for efficiency) // eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) ) eps(0,0) = gradu(0,0); eps(0,1) = .5 * ( gradu(0,1) + gradu(1,0) ); eps(1,0) = eps (0,1); eps(1,1) = gradu(1,1); // - copy resulting integration point tensor std::copy(eps.begin(), eps.end(), qtensor.item(e,k)); } } } // #pragma omp parallel return qtensor; } // ------- scalar product "elemmat(m*ndim+i,n*ndim+i) = N(m) * qscalar * N(n)"; for all "i" -------- inline ArrD Quadrature::int_N_scalar_NT_dV(const ArrD &qscalar) const { // check input assert( qscalar.ndim() == 2 ); // shape: [nelem, nip] assert( qscalar.shape(0) == m_nelem ); // number of elements assert( qscalar.shape(1) == m_nip ); // number of integration points // zero-initialize: matrix of matrices ArrD elemmat = ArrD::Zero({m_nelem, m_nne*m_ndim, m_nne*m_ndim}); #pragma omp parallel { // intermediate quantities and local views - cppmat::tiny::matrix M; - cppmat::view::vector N; + cppmat::tiny::matrix M; + cppmat::view::vector N; double rho, vol; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // zero-initialize (e.g. mass matrix) M.setZero(); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias N.setMap(&m_N( k)); // shape functions vol = m_vol (e,k); // integration point volume rho = qscalar(e,k); // integration point scalar (e.g. density) // - 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; } } } // copy result to element matrix std::copy(M.begin(), M.end(), elemmat.item(e)); } } // #pragma omp parallel return elemmat; } // ------------ integral of dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" ------------ template inline ArrD Quadrature::int_gradN_dot_tensor2_dV(const ArrD &qtensor) const { #ifndef NDEBUG // dummy variable T tmp; // check input assert( qtensor.ndim() == 3 ); // shape: [nelem, nip, #tensor-components] assert( qtensor.shape(0) == m_nelem ); // number of elements assert( qtensor.shape(1) == m_nip ); // number of integration points assert( qtensor.shape(2) == tmp.size() ); // tensor dimensions #endif // zero-initialize output: matrix of vectors ArrD elemvec = ArrD::Zero({m_nelem, m_nne, m_ndim}); #pragma omp parallel { // intermediate quantities and local views - cppmat::view::matrix dNx; - cppmat::tiny::matrix f; + cppmat::view::matrix dNx; + cppmat::tiny::matrix f; double vol; T sig; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // zero-initialize (e.g. nodal force) f.setZero(); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias dNx.setMap(&m_dNx (e,k)); // shape function gradients (global coordinates) sig.setMap(&qtensor(e,k)); // integration point tensor (e.g. stress) vol = m_vol (e,k); // integration point volume // - evaluate dot product, and assemble (loops partly unrolled for efficiency) // f(m,j) += dNdx(m,i) * sig(i,j) * dV; 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; } } // copy result to element vector std::copy(f.begin(), f.end(), elemvec.item(e)); } } // #pragma omp parallel return elemvec; } // ---------------------- wrappers with default storage (no template needed) ----------------------- inline ArrD Quadrature::gradN_vector(const ArrD &elemvec) const { return gradN_vector>(elemvec); } // ------------------------------------------------------------------------------------------------- inline ArrD Quadrature::gradN_vector_T(const ArrD &elemvec) const { return gradN_vector_T>(elemvec); } // ------------------------------------------------------------------------------------------------- inline ArrD Quadrature::symGradN_vector(const ArrD &elemvec) const { return symGradN_vector>(elemvec); } // ------------------------------------------------------------------------------------------------- inline ArrD Quadrature::int_gradN_dot_tensor2_dV(const ArrD &qtensor) const { assert( qtensor.ndim() == 3 ); // shape: [nelem, nip, #tensor-components] if ( qtensor.shape(2) == m_ndim*m_ndim ) return int_gradN_dot_tensor2_dV>(qtensor); else if ( qtensor.shape(2) == (m_ndim+1)*m_ndim/2 ) return int_gradN_dot_tensor2_dV>(qtensor); else throw std::runtime_error("assert: qtensor.shape(2) == 4 or qtensor.shape(2) == 3"); } // ------------------------------------------------------------------------------------------------- inline ArrD Quadrature::int_gradN_dot_tensor2s_dV(const ArrD &qtensor) const { return int_gradN_dot_tensor2_dV>(qtensor); } // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/ElementQuad4.h b/src/GooseFEM/ElementQuad4.h index 5c1d72d..3021a7c 100644 --- a/src/GooseFEM/ElementQuad4.h +++ b/src/GooseFEM/ElementQuad4.h @@ -1,125 +1,125 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_ELEMENTQUAD4_H #define GOOSEFEM_ELEMENTQUAD4_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =================================== GooseFEM::Element::Quad4 ==================================== namespace GooseFEM { namespace Element { namespace Quad4 { // ================================ GooseFEM::Element::Quad4::Gauss ================================ namespace Gauss { inline size_t nip(); // number of integration points inline ArrD xi(); // integration point coordinates (local coordinates) inline ArrD w(); // integration point weights } // ================================ GooseFEM::Element::Quad4::Nodal ================================ namespace Nodal { inline size_t nip(); // number of integration points inline ArrD xi(); // integration point coordinates (local coordinates) inline ArrD 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 + size_t m_nne=4; // number of nodes per element + size_t m_ndim=2; // number of dimensions // data arrays ArrD m_x; // nodal positions stored per element [nelem, nne, ndim] ArrD m_w; // weight of each integration point [nip] ArrD m_xi; // local coordinate of each integration point [nip, ndim] ArrD m_N; // shape functions [nip, nne] ArrD m_dNxi; // shape function gradients w.r.t. local coordinate [nip, nne, ndim] ArrD m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] ArrD m_vol; // integration point volume [nelem, nip] private: // compute "vol" and "dNdx" based on current "x" void compute_dN(); public: // convention: // "elemmat" - matrices stored per element - ArrD - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - ArrD - [nelem, nne, ndim] // "qtensor" - integration point tensor - ArrD - [nelem, nip, #tensor-components] // "qscalar" - integration point scalar - ArrD - [nelem, nip] // // alias: // T2 = cppmat::cartesian2d::tensor2 - #tensor-components = 4 // T2s = cppmat::cartesian2d::tensor2s - #tensor-components = 3 // vT2 = cppmat::view::cartesian2d::tensor2 - #tensor-components = 4 // vT2s = cppmat::view::cartesian2d::tensor2s - #tensor-components = 3 // constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature(){}; Quadrature(const ArrD &x, const ArrD &xi=ArrD(), const ArrD &w=ArrD()); // update the nodal positions (shape of "x" should match the earlier definition) void update_x(const ArrD &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 // return integration volume ArrD dV() const; // returns: qscalar ArrD dV(size_t ncomp) const; // returns: qtensor (same volume per tensor-component) // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part // - allow template (e.g. T2/T2s, or higher dimensional tensors) template ArrD gradN_vector (const ArrD &elemvec) const; // returns: qtensor template ArrD gradN_vector_T (const ArrD &elemvec) const; // returns: qtensor template ArrD symGradN_vector(const ArrD &elemvec) const; // returns: qtensor // - default template ArrD gradN_vector (const ArrD &elemvec) const; // template: T2 ArrD gradN_vector_T (const ArrD &elemvec) const; // template: T2 ArrD symGradN_vector(const ArrD &elemvec) const; // template: T2s // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV" ArrD int_N_scalar_NT_dV(const ArrD &qscalar) const; // returns: elemmat // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" // - allow template (e.g. vT2/vT2s, or higher dimensional tensors) template ArrD int_gradN_dot_tensor2_dV(const ArrD &qtensor) const; // returns: elemvec // - default template ArrD int_gradN_dot_tensor2_dV (const ArrD &qtensor) const; // template: vT2/vT2s (auto-select) ArrD int_gradN_dot_tensor2s_dV(const ArrD &qtensor) const; // template: vT2s }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif