diff --git a/docs/GooseFEM.svg b/docs/GooseFEM.svg new file mode 100644 index 0000000..e786884 --- /dev/null +++ b/docs/GooseFEM.svg @@ -0,0 +1,293 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + GooseFEM + + + diff --git a/include/xGooseFEM/Dynamics.h b/include/xGooseFEM/Dynamics.h index 8578a00..66afee2 100644 --- a/include/xGooseFEM/Dynamics.h +++ b/include/xGooseFEM/Dynamics.h @@ -1,64 +1,67 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_DYNAMICS_H #define XGOOSEFEM_DYNAMICS_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // ======================================= GooseFEM::Dynamics ======================================= namespace xGooseFEM { namespace Dynamics { // ------------------------------------------ dummy class ------------------------------------------ class Geometry { public: + // destructor + virtual ~Geometry() {}; + // solve for DOF-accelerations [ndof] virtual xt::xtensor solve_A() { return xt::empty({0}); }; virtual xt::xtensor solve_V() { return xt::empty({0}); }; // return nodal vectors [nnode, ndim] virtual xt::xtensor u() const { return xt::empty({0,0}); }; virtual xt::xtensor v() const { return xt::empty({0,0}); }; virtual xt::xtensor a() const { return xt::empty({0,0}); }; // return DOF values [ndof] virtual xt::xtensor dofs_u() const { return xt::empty({0}); }; virtual xt::xtensor dofs_v() const { return xt::empty({0}); }; virtual xt::xtensor dofs_a() const { return xt::empty({0}); }; // overwrite nodal vectors [nnode, ndim] virtual void set_u(const xt::xtensor &nodevec) { UNUSED(nodevec); return; }; // overwrite nodal vectors, reconstructed from DOF values [ndof] virtual void set_u(const xt::xtensor &dofval) { UNUSED(dofval); return; }; virtual void set_v(const xt::xtensor &dofval) { UNUSED(dofval); return; }; virtual void set_a(const xt::xtensor &dofval) { UNUSED(dofval); return; }; }; // ------------------------------------ evaluate one time step ------------------------------------- inline void Verlet (Geometry &geometry, double dt, size_t nstep=1); inline void velocityVerlet(Geometry &geometry, double dt, size_t nstep=1); namespace Overdamped { inline void forwardEuler (Geometry &geometry, double dt, size_t nstep=1); } // ------------------------------------------------------------------------------------------------- }} // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementQuad4.h b/include/xGooseFEM/ElementQuad4.h index 099dde1..1bac6d1 100644 --- a/include/xGooseFEM/ElementQuad4.h +++ b/include/xGooseFEM/ElementQuad4.h @@ -1,132 +1,156 @@ /* ================================================================================================= (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: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - [nelem, nip] // constructor: integration point coordinates and weights are optional (default: Gauss) + Quadrature() = default; + Quadrature(const xt::xtensor &x); + Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w); // update the nodal positions (shape of "x" should match the earlier definition) + void update_x(const xt::xtensor &x); // return 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 - // - in-place + void dV(xt::xtensor &qscalar) const; void dV(xt::xtensor &qtensor) const; // same volume for all tensor components - // - return qscalar/qtensor + + // 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; + }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementQuad4.hpp b/include/xGooseFEM/ElementQuad4.hpp index 8bf8611..1550a14 100644 --- a/include/xGooseFEM/ElementQuad4.hpp +++ b/include/xGooseFEM/ElementQuad4.hpp @@ -1,642 +1,702 @@ /* ================================================================================================= (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" -// =================================== GooseFEM::Element::Quad4 ==================================== +// ================================================================================================= namespace xGooseFEM { namespace Element { namespace Quad4 { -// ======================================== tensor algebra ========================================= +// ================================================================================================= inline double inv(const T2 &A, T2 &Ainv) { // compute determinant double det = A[0] * A[3] - A[1] * A[2]; // compute inverse Ainv[0] = A[3] / det; Ainv[1] = -1. * A[1] / det; Ainv[2] = -1. * A[2] / det; Ainv[3] = A[0] / det; return det; } -// ================================ GooseFEM::Element::Quad4::Gauss ================================ +// ================================================================================================= namespace Gauss { -// --------------------------------- number of integration points ---------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t nip() { return 4; } -// ----------------------- integration point coordinates (local coordinates) ----------------------- +// ------------------------------------------------------------------------------------------------- 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; } -// ----------------------------------- integration point weights ----------------------------------- +// ------------------------------------------------------------------------------------------------- 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; } // ------------------------------------------------------------------------------------------------- } -// ================================ GooseFEM::Element::Quad4::Nodal ================================ +// ================================================================================================= namespace Nodal { -// --------------------------------- number of integration points ---------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t nip() { return 4; } -// ----------------------- integration point coordinates (local coordinates) ----------------------- +// ------------------------------------------------------------------------------------------------- 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; } -// ----------------------------------- integration point weights ----------------------------------- +// ------------------------------------------------------------------------------------------------- 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; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= -// ------------------------------------------ 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) = .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(); } -// ------------------------------------------ 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) { 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(); 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 ) { 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(); } -// --------------------------- 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); } // ------------------------------------------------------------------------------------------------- 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 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 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 ) { // - alias auto dNxi = xt::adapt(&m_dNxi( k,0,0), xt::xshape()); auto dNx = xt::adapt(&m_dNx (e,k,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; } } } } -// ------------------- 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 ) { // - alias auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e,k,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 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 ) { // - alias auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e,k,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 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 ) { // - alias auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape()); auto eps = xt::adapt(&qtensor(e,k,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.; eps(1,0) = eps(0,1); } } } -// ------------------------------------------------------------------------------------------------- - -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 ) { // - alias auto N = xt::adapt(&m_N(k,0), xt::xshape()); auto& vol = m_vol (e,k); auto& rho = qscalar(e,k); // - 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; } } } } } -// ------------------------------------------------------------------------------------------------- - -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 ) { // - 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); // - 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 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/GooseFEM.h b/include/xGooseFEM/GooseFEM.h index af1a7c8..d885316 100644 --- a/include/xGooseFEM/GooseFEM.h +++ b/include/xGooseFEM/GooseFEM.h @@ -1,104 +1,111 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_H #define XGOOSEFEM_H // ================================================================================================= #define _USE_MATH_DEFINES // to use "M_PI" from "math.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace xt::placeholders; // ================================================================================================= #define GOOSEFEM_WORLD_VERSION 0 #define GOOSEFEM_MAJOR_VERSION 1 #define GOOSEFEM_MINOR_VERSION 0 #define GOOSEFEM_VERSION_AT_LEAST(x,y,z) \ (GOOSEFEM_WORLD_VERSION>x || (GOOSEFEM_WORLD_VERSION>=x && \ (GOOSEFEM_MAJOR_VERSION>y || (GOOSEFEM_MAJOR_VERSION>=y && \ GOOSEFEM_MINOR_VERSION>=z)))) #define GOOSEFEM_VERSION(x,y,z) \ (GOOSEFEM_WORLD_VERSION==x && \ GOOSEFEM_MAJOR_VERSION==y && \ GOOSEFEM_MINOR_VERSION==z) // ================================================================================================= // dummy operation that can be use to suppress the "unused parameter" warnings #define UNUSED(p) ( (void)(p) ) // ================================================================================================= // alias Eigen sparse matrices namespace xGooseFEM { typedef Eigen::SparseMatrix SpMatD; typedef Eigen::SparseMatrix SpMatS; + typedef Eigen::Triplet TripD; } // ================================================================================================= #include "Mesh.h" #include "MeshTri3.h" #include "MeshQuad4.h" #include "MeshHex8.h" #include "Element.h" #include "ElementQuad4.h" #include "ElementHex8.h" #include "Vector.h" +#include "VectorPartitioned.h" +#include "MatrixPartitioned.h" #include "MatrixDiagonal.h" +#include "MatrixDiagonalPartitioned.h" #include "Iterate.h" #include "Dynamics.h" #include "Mesh.hpp" #include "MeshTri3.hpp" #include "MeshQuad4.hpp" #include "MeshHex8.hpp" #include "Element.hpp" #include "ElementQuad4.hpp" #include "ElementHex8.hpp" #include "Vector.hpp" +#include "VectorPartitioned.hpp" +#include "MatrixPartitioned.hpp" #include "MatrixDiagonal.hpp" +#include "MatrixDiagonalPartitioned.hpp" #include "Iterate.hpp" #include "Dynamics.hpp" // ================================================================================================= #endif diff --git a/include/xGooseFEM/MatrixDiagonal.h b/include/xGooseFEM/MatrixDiagonal.h index 862c400..0390489 100644 --- a/include/xGooseFEM/MatrixDiagonal.h +++ b/include/xGooseFEM/MatrixDiagonal.h @@ -1,147 +1,81 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_MATRIXDIAGONAL_H #define XGOOSEFEM_MATRIXDIAGONAL_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =========================================== GooseFEM ============================================ namespace xGooseFEM { // ------------------------------------------------------------------------------------------------- class MatrixDiagonal { public: // constructor MatrixDiagonal() = default; - MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs); - MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs, - const xt::xtensor &iip); - - // index operators: access plain storage - double& operator[](size_t i); - const double& operator[](size_t i) const; - - // index operators: access using matrix indices - double& operator()(size_t i); - const double& operator()(size_t i) const; - double& operator()(size_t i, size_t j); - const double& operator()(size_t i, size_t j) const; - // dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs - size_t nnu() const; // number of unknown DOFs - size_t nnp() const; // number of prescribed DOFs // DOF lists xt::xtensor dofs() const; // DOFs - xt::xtensor iiu() const; // unknown DOFs - xt::xtensor iip() const; // prescribed DOFs // product: b_i = A_ij * x_j - // b = A * x - xt::xtensor dot (const xt::xtensor &x ) const; - // b_u = A_uu * x_u + A_up * x_p - xt::xtensor dot_u(const xt::xtensor &x_u, const xt::xtensor &x_p) const; - // b_p = A_pu * x_u + A_pp * x_p - xt::xtensor dot_p(const xt::xtensor &x_u, const xt::xtensor &x_p) const; + xt::xtensor dot(const xt::xtensor &x) const; // assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] // WARNING: ignores any off-diagonal terms void assemble(const xt::xtensor &elemmat); - // set matrix components from externally assembled object - void set (const xt::xtensor &A ) const; // diagonal [ndof] - void set_uu(const xt::xtensor &A_uu) const; // diagonal [nnu] - void set_pp(const xt::xtensor &A_pp) const; // diagonal [nnp] - - // solve - // x = A \ b - xt::xtensor solve (const xt::xtensor &b ); - // x = assembly{ A_uu \ ( b_u - A_up * x_p ) ; x_p } - xt::xtensor solve (const xt::xtensor &b , const xt::xtensor &x_p); - // x_u = A_uu \ ( b_u - A_up * x_p ) - xt::xtensor solve_u(const xt::xtensor &b_u, const xt::xtensor &x_p); - - // return (sub-)matrix as diagonal matrix (column) - // assembly{ A_uu ; A_pp } - xt::xtensor asDiagonal () const; - // A_uu - xt::xtensor asDiagonal_uu() const; - // A_pp - xt::xtensor asDiagonal_pp() const; - - // return (sub-)matrix as sparse matrix - // assembly{ A_uu ; A_pp } - SpMatD asSparse () const; - // A_uu - SpMatD asSparse_uu() const; - // empty - SpMatD asSparse_up() const; - // empty - SpMatD asSparse_pu() const; - // A_pp - SpMatD asSparse_pp() const; - - // return (sub-)matrix as dense matrix - // assembly{ A_uu ; A_pp } - xt::xtensor asDense () const; - // A_uu - xt::xtensor asDense_uu() const; - // empty - xt::xtensor asDense_up() const; - // empty - xt::xtensor asDense_pu() const; - // A_pp - xt::xtensor asDense_pp() const; + // solve: x = A \ b + xt::xtensor solve(const xt::xtensor &b); + + // return matrix as diagonal matrix (column) + xt::xtensor asDiagonal() const; private: // the diagonal matrix (not-partitioned), and its inverse (re-used to solve different RHS) xt::xtensor m_data; xt::xtensor m_inv; // signal changes to data compare to the last inverse bool m_change=false; // bookkeeping - xt::xtensor m_conn; // connectivity [nelem, nne ] - xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] - xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] - xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] + xt::xtensor m_conn; // connectivity [nelem, nne ] + xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] // dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs - size_t m_nnu; // number of unknown DOFs - size_t m_nnp; // number of prescribed DOFs // compute inverse (automatically evaluated by "solve") void factorize(); }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/MatrixDiagonal.hpp b/include/xGooseFEM/MatrixDiagonal.hpp index 8223f52..4792cd4 100644 --- a/include/xGooseFEM/MatrixDiagonal.hpp +++ b/include/xGooseFEM/MatrixDiagonal.hpp @@ -1,378 +1,163 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_MATRIXDIAGONAL_CPP #define XGOOSEFEM_MATRIXDIAGONAL_CPP // ------------------------------------------------------------------------------------------------- #include "MatrixDiagonal.h" -// =========================================== GooseFEM ============================================ +// ================================================================================================= namespace xGooseFEM { -// ------------------------------------------ constructor ------------------------------------------ +// ------------------------------------------------------------------------------------------------- inline MatrixDiagonal::MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs) : m_conn(conn), m_dofs(dofs) { // mesh dimensions m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; - // list with prescribed and unknown DOFs - m_iip = xt::empty({0}); - m_iiu = xt::unique(dofs); - // dimensions of the system m_ndof = xt::amax(m_dofs)[0] + 1; - m_nnp = m_iip.size(); - m_nnu = m_iiu.size(); // allocate matrix and its inverse m_data = xt::empty({m_ndof}); m_inv = xt::empty({m_ndof}); // check consistency - assert( xt::amax(m_conn)[0] + 1 == m_nnode ); - assert( m_ndof <= m_nnode * m_ndim ); -} - -// ------------------------------------------ constructor ------------------------------------------ - -inline MatrixDiagonal::MatrixDiagonal(const xt::xtensor &conn, - const xt::xtensor &dofs, const xt::xtensor &iip) : - m_conn(conn), m_dofs(dofs), m_iip(iip) -{ - // mesh dimensions - m_nelem = m_conn.shape()[0]; - m_nne = m_conn.shape()[1]; - m_nnode = m_dofs.shape()[0]; - m_ndim = m_dofs.shape()[1]; - - // list with unknown DOFs - m_iiu = xt::setdiff1d(dofs, iip); - - // dimensions of the system - m_ndof = xt::amax(m_dofs)[0] + 1; - m_nnp = m_iip.size(); - m_nnu = m_iiu.size(); - - // allocate matrix and its inverse - m_data = xt::empty({m_ndof}); - m_inv = xt::empty({m_ndof}); - - // check consistency - assert( xt::amax(m_conn)[0] + 1 == m_nnode ); - assert( xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0] ); - assert( m_ndof <= m_nnode * m_ndim ); -} - -// ---------------------------------------- index operator ----------------------------------------- - -inline double& MatrixDiagonal::operator[](size_t i) -{ - m_change = true; - - return m_data[i]; -} - -// ---------------------------------------- index operator ----------------------------------------- - -inline const double& MatrixDiagonal::operator[](size_t i) const -{ - return m_data[i]; + assert( xt::amax(m_conn)[0] + 1 == m_nnode ); + assert( m_ndof <= m_nnode * m_ndim ); } -// ---------------------------------------- index operator ----------------------------------------- - -inline double& MatrixDiagonal::operator()(size_t i) -{ - m_change = true; - - return m_data[i]; -} - -// ---------------------------------------- index operator ----------------------------------------- - -inline const double& MatrixDiagonal::operator()(size_t i) const -{ - return m_data[i]; -} - -// ---------------------------------------- index operator ----------------------------------------- - -inline double& MatrixDiagonal::operator()(size_t i, size_t j) -{ - assert( i == j ); - - UNUSED(j); - - m_change = true; - - return m_data[i]; -} - -// ---------------------------------------- index operator ----------------------------------------- - -inline const double& MatrixDiagonal::operator()(size_t i, size_t j) const -{ - assert( i == j ); - - UNUSED(j); - - return m_data[i]; -} - -// -------------------------------------- number of elements --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t MatrixDiagonal::nelem() const { return m_nelem; } -// ---------------------------------- number of nodes per element ---------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t MatrixDiagonal::nne() const { return m_nne; } -// ---------------------------------------- number of nodes ---------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t MatrixDiagonal::nnode() const { return m_nnode; } -// ------------------------------------- number of dimensions -------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t MatrixDiagonal::ndim() const { return m_ndim; } -// ---------------------------------------- number of DOFs ----------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t MatrixDiagonal::ndof() const { return m_ndof; } -// ------------------------------------ number of unknown DOFs ------------------------------------- - -inline size_t MatrixDiagonal::nnu() const -{ - return m_nnu; -} - -// ----------------------------------- number of prescribed DOFs ----------------------------------- - -inline size_t MatrixDiagonal::nnp() const -{ - return m_nnp; -} - -// ------------------------------------------ return DOFs ------------------------------------------ +// ------------------------------------------------------------------------------------------------- inline xt::xtensor MatrixDiagonal::dofs() const { return m_dofs; } -// -------------------------------------- return unknown DOFs -------------------------------------- - -inline xt::xtensor MatrixDiagonal::iiu() const -{ - return m_iiu; -} - -// ------------------------------------ return prescribed DOFs ------------------------------------- - -inline xt::xtensor MatrixDiagonal::iip() const -{ - return m_iip; -} - -// --------------------------------------- b_i = A_ij * x_j ---------------------------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor MatrixDiagonal::dot(const xt::xtensor &x) const { // check input assert( x.size() == m_ndof ); // compute product return m_data * x; } -// --------------------------------------- b_i = A_ij * x_j ---------------------------------------- - -inline xt::xtensor MatrixDiagonal::dot_u( - const xt::xtensor &x_u, const xt::xtensor &x_p) const -{ - // suppress warning - UNUSED(x_p); - - // check input - assert( x_u.size() == m_nnu ); - assert( x_p.size() == m_nnp ); - - // allocate output - xt::xtensor b_u = xt::empty({m_nnu}); - - // compute product - #pragma omp parallel for - for ( size_t i = 0 ; i < m_nnu ; ++i ) - b_u(i) = m_data(m_iiu(i)) * x_u(i); - - return b_u; -} - -// --------------------------------------- b_i = A_ij * x_j ---------------------------------------- - -inline xt::xtensor MatrixDiagonal::dot_p( - const xt::xtensor &x_u, const xt::xtensor &x_p) const -{ - // suppress warning - UNUSED(x_u); - - // check input - assert( x_u.size() == m_nnu ); - assert( x_p.size() == m_nnp ); - - // allocate output - xt::xtensor b_p = xt::empty({m_nnp}); - - // compute product - #pragma omp parallel for - for ( size_t i = 0 ; i < m_nnp ; ++i ) - b_p(i) = m_data(m_iip(i)) * x_p(i); - - return b_p; -} - // ----------------------------- assemble matrices stored per element ------------------------------ inline void MatrixDiagonal::assemble(const xt::xtensor &elemmat) { // check input assert( elemmat.shape()[0] == m_nelem ); assert( elemmat.shape()[1] == m_nne*m_ndim ); assert( elemmat.shape()[2] == m_nne*m_ndim ); assert( Element::isDiagonal(elemmat) ); // zero-initialize matrix m_data.fill(0.0); // assemble for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) m_data(m_dofs(m_conn(e,m),i)) += elemmat(e,m*m_ndim+i,m*m_ndim+i); // signal change m_change = true; } // ------------------------------------------------------------------------------------------------- inline void MatrixDiagonal::factorize() { // skip for unchanged "m_data" if ( ! m_change ) return; - // invert if needed + // invert #pragma omp parallel for for ( size_t i = 0 ; i < m_ndof ; ++i ) m_inv(i) = 1. / m_data(i); // reset signal m_change = false; } -// --------------------------------------- solve: A * x = b ---------------------------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor MatrixDiagonal::solve(const xt::xtensor &b) { // check input - assert( m_nnp == 0 ); assert( b.size() == m_ndof ); - // compute inverse - this->factorize(); - - // solve - xt::xtensor x = m_inv * b; - - return x; -} - -// --------------------------------------- solve: A * x = b ---------------------------------------- - -inline xt::xtensor MatrixDiagonal::solve( - const xt::xtensor &b, const xt::xtensor &x_p) -{ - // check input - assert( x_p.size() == m_nnp ); - assert( b .size() == m_ndof ); - - // compute inverse + // factorise (if needed) this->factorize(); // solve xt::xtensor x = m_inv * b; - // set prescribed DOFs - #pragma omp parallel for - for ( size_t i = 0 ; i < m_nnp ; ++i ) - x(m_iip(i)) = x_p(i); - + // return output return x; } -// --------------------------------------- solve: A * x = b ---------------------------------------- - -inline xt::xtensor MatrixDiagonal::solve_u( - const xt::xtensor &b_u, const xt::xtensor &x_p) -{ - // suppress warning - UNUSED(x_p); - - // check input - assert( x_p.size() == m_nnp ); - assert( b_u.size() == m_nnu ); - - // compute inverse - this->factorize(); - - // allocate output - xt::xtensor x_u = xt::empty({m_nnu}); - - // solve - #pragma omp parallel for - for ( size_t i = 0 ; i < m_nnu ; ++i ) - x_u(i) = m_inv(m_iiu(i)) * b_u(i); - - return x_u; -} - -// ----------------------------------- return as diagonal matrix ----------------------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor MatrixDiagonal::asDiagonal() const { return m_data; } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/MatrixDiagonal.h b/include/xGooseFEM/MatrixDiagonalPartitioned.h similarity index 53% copy from include/xGooseFEM/MatrixDiagonal.h copy to include/xGooseFEM/MatrixDiagonalPartitioned.h index 862c400..c372bb9 100644 --- a/include/xGooseFEM/MatrixDiagonal.h +++ b/include/xGooseFEM/MatrixDiagonalPartitioned.h @@ -1,147 +1,99 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef XGOOSEFEM_MATRIXDIAGONAL_H -#define XGOOSEFEM_MATRIXDIAGONAL_H +#ifndef XGOOSEFEM_MATRIXDIAGONALPARTITIONED_H +#define XGOOSEFEM_MATRIXDIAGONALPARTITIONED_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =========================================== GooseFEM ============================================ namespace xGooseFEM { // ------------------------------------------------------------------------------------------------- -class MatrixDiagonal +class MatrixDiagonalPartitioned { public: - // constructor - MatrixDiagonal() = default; - - MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs); + // default constructor + MatrixDiagonalPartitioned() = default; - MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs, + // constructor + MatrixDiagonalPartitioned(const xt::xtensor &conn, const xt::xtensor &dofs, const xt::xtensor &iip); - // index operators: access plain storage - double& operator[](size_t i); - const double& operator[](size_t i) const; - - // index operators: access using matrix indices - double& operator()(size_t i); - const double& operator()(size_t i) const; - double& operator()(size_t i, size_t j); - const double& operator()(size_t i, size_t j) const; - // dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs size_t nnu() const; // number of unknown DOFs size_t nnp() const; // number of prescribed DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // unknown DOFs xt::xtensor iip() const; // prescribed DOFs // product: b_i = A_ij * x_j // b = A * x xt::xtensor dot (const xt::xtensor &x ) const; // b_u = A_uu * x_u + A_up * x_p xt::xtensor dot_u(const xt::xtensor &x_u, const xt::xtensor &x_p) const; // b_p = A_pu * x_u + A_pp * x_p xt::xtensor dot_p(const xt::xtensor &x_u, const xt::xtensor &x_p) const; // assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] // WARNING: ignores any off-diagonal terms void assemble(const xt::xtensor &elemmat); - // set matrix components from externally assembled object - void set (const xt::xtensor &A ) const; // diagonal [ndof] - void set_uu(const xt::xtensor &A_uu) const; // diagonal [nnu] - void set_pp(const xt::xtensor &A_pp) const; // diagonal [nnp] - - // solve - // x = A \ b - xt::xtensor solve (const xt::xtensor &b ); - // x = assembly{ A_uu \ ( b_u - A_up * x_p ) ; x_p } - xt::xtensor solve (const xt::xtensor &b , const xt::xtensor &x_p); - // x_u = A_uu \ ( b_u - A_up * x_p ) - xt::xtensor solve_u(const xt::xtensor &b_u, const xt::xtensor &x_p); - - // return (sub-)matrix as diagonal matrix (column) - // assembly{ A_uu ; A_pp } - xt::xtensor asDiagonal () const; - // A_uu - xt::xtensor asDiagonal_uu() const; - // A_pp - xt::xtensor asDiagonal_pp() const; - - // return (sub-)matrix as sparse matrix - // assembly{ A_uu ; A_pp } - SpMatD asSparse () const; - // A_uu - SpMatD asSparse_uu() const; - // empty - SpMatD asSparse_up() const; - // empty - SpMatD asSparse_pu() const; - // A_pp - SpMatD asSparse_pp() const; - - // return (sub-)matrix as dense matrix - // assembly{ A_uu ; A_pp } - xt::xtensor asDense () const; - // A_uu - xt::xtensor asDense_uu() const; - // empty - xt::xtensor asDense_up() const; - // empty - xt::xtensor asDense_pu() const; - // A_pp - xt::xtensor asDense_pp() const; + // solve: x_u = A_uu \ ( b_u - A_up * x_p ) + xt::xtensor solve(const xt::xtensor &b_u, const xt::xtensor &x_p); + + // return matrix as diagonal matrix + xt::xtensor asDiagonal() const; private: - // the diagonal matrix (not-partitioned), and its inverse (re-used to solve different RHS) - xt::xtensor m_data; - xt::xtensor m_inv; + // the diagonal matrix, and its inverse (re-used to solve different RHS) + xt::xtensor m_data_uu; + xt::xtensor m_data_pp; + xt::xtensor m_inv_uu; // signal changes to data compare to the last inverse bool m_change=false; // bookkeeping - xt::xtensor m_conn; // connectivity [nelem, nne ] - xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] - xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] - xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] + xt::xtensor m_conn; // connectivity [nelem, nne ] + xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] + xt::xtensor m_part; // DOF-numbers per node, renumbered [nnode, ndim] + xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] + xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] // dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs size_t m_nnu; // number of unknown DOFs size_t m_nnp; // number of prescribed DOFs // compute inverse (automatically evaluated by "solve") void factorize(); }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/MatrixDiagonalPartitioned.hpp b/include/xGooseFEM/MatrixDiagonalPartitioned.hpp new file mode 100644 index 0000000..aa19c00 --- /dev/null +++ b/include/xGooseFEM/MatrixDiagonalPartitioned.hpp @@ -0,0 +1,297 @@ +/* ================================================================================================= + +(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM + +================================================================================================= */ + +#ifndef XGOOSEFEM_MATRIXDIAGONALPARTITIONED_CPP +#define XGOOSEFEM_MATRIXDIAGONALPARTITIONED_CPP + +// ------------------------------------------------------------------------------------------------- + +#include "MatrixDiagonalPartitioned.h" + +// ================================================================================================= + +namespace xGooseFEM { + +// ------------------------------------------------------------------------------------------------- + +inline MatrixDiagonalPartitioned::MatrixDiagonalPartitioned(const xt::xtensor &conn, + const xt::xtensor &dofs, const xt::xtensor &iip) : + m_conn(conn), m_dofs(dofs), m_iip(iip) +{ + // mesh dimensions + m_nelem = m_conn.shape()[0]; + m_nne = m_conn.shape()[1]; + m_nnode = m_dofs.shape()[0]; + m_ndim = m_dofs.shape()[1]; + + // list with unknown DOFs + m_iiu = xt::setdiff1d(dofs, iip); + + // dimensions of the system + m_ndof = xt::amax(m_dofs)[0] + 1; + m_nnp = m_iip.size(); + m_nnu = m_iiu.size(); + + // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) + m_part = Mesh::reorder(m_dofs, m_iip, "end"); + + // allocate matrix and its inverse + m_data_uu = xt::empty({m_nnu}); + m_data_pp = xt::empty({m_nnp}); + m_inv_uu = xt::empty({m_nnu}); + + // check consistency + assert( xt::amax(m_conn)[0] + 1 == m_nnode ); + assert( xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0] ); + assert( m_ndof <= m_nnode * m_ndim ); +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::nelem() const +{ + return m_nelem; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::nne() const +{ + return m_nne; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::nnode() const +{ + return m_nnode; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::ndim() const +{ + return m_ndim; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::ndof() const +{ + return m_ndof; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::nnu() const +{ + return m_nnu; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixDiagonalPartitioned::nnp() const +{ + return m_nnp; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::dofs() const +{ + return m_dofs; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::iiu() const +{ + return m_iiu; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::iip() const +{ + return m_iip; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::dot(const xt::xtensor &x) const +{ + // check input + assert( x.size() == m_ndof ); + + // allocate output + xt::xtensor b = xt::empty({m_ndof}); + + // compute product + // - + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnu ; ++i ) + b(m_iiu(i)) = m_data_uu(i) * x(m_iiu(i)); + // - + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnp ; ++i ) + b(m_iip(i)) = m_data_pp(i) * x(m_iip(i)); + + // return output + return b; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::dot_u( + const xt::xtensor &x_u, const xt::xtensor &x_p) const +{ + // suppress warning + UNUSED(x_p); + + // check input + assert( x_u.size() == m_nnu ); + assert( x_p.size() == m_nnp ); + + // allocate output + xt::xtensor b_u = xt::empty({m_nnu}); + + // compute product + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnu ; ++i ) + b_u(i) = m_data_uu(i) * x_u(i); + + return b_u; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::dot_p( + const xt::xtensor &x_u, const xt::xtensor &x_p) const +{ + // suppress warning + UNUSED(x_u); + + // check input + assert( x_u.size() == m_nnu ); + assert( x_p.size() == m_nnp ); + + // allocate output + xt::xtensor b_p = xt::empty({m_nnp}); + + // compute product + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnp ; ++i ) + b_p(i) = m_data_pp(i) * x_p(i); + + return b_p; +} + +// ------------------------------------------------------------------------------------------------- + +inline void MatrixDiagonalPartitioned::assemble(const xt::xtensor &elemmat) +{ + // check input + assert( elemmat.shape()[0] == m_nelem ); + assert( elemmat.shape()[1] == m_nne*m_ndim ); + assert( elemmat.shape()[2] == m_nne*m_ndim ); + assert( Element::isDiagonal(elemmat) ); + + // zero-initialize matrix + m_data_uu.fill(0.0); + m_data_pp.fill(0.0); + + // assemble + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + for ( size_t m = 0 ; m < m_nne ; ++m ) + { + for ( size_t i = 0 ; i < m_ndim ; ++i ) + { + size_t di = m_part(m_conn(e,m),i); + + if ( di < m_nnu ) + m_data_uu(di ) += elemmat(e,m*m_ndim+i,m*m_ndim+i); + else + m_data_pp(di-m_nnu) += elemmat(e,m*m_ndim+i,m*m_ndim+i); + } + } + } + + // signal change + m_change = true; +} + +// ------------------------------------------------------------------------------------------------- + +inline void MatrixDiagonalPartitioned::factorize() +{ + // skip for unchanged "m_data" + if ( ! m_change ) return; + + // invert + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnu ; ++i ) + m_inv_uu(i) = 1. / m_data_uu(i); + + // reset signal + m_change = false; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::solve( + const xt::xtensor &b_u, const xt::xtensor &x_p) +{ + // suppress warning + UNUSED(x_p); + + // check input + assert( b_u.shape()[0] == m_nnu ); + assert( x_p.shape()[0] == m_nnp ); + + // factorise (if needed) + this->factorize(); + + // allocate output + xt::xtensor x_u = xt::empty({m_nnu}); + + // solve + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnu ; ++i ) + x_u(i) = m_inv_uu(i) * b_u(i); + + // return output + return x_u; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixDiagonalPartitioned::asDiagonal() const +{ + // allocate output + xt::xtensor out = xt::zeros({m_ndof}); + + // assemble output + // - + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnu ; ++i ) + out(m_iiu(i)) = m_data_uu(i); + // - + #pragma omp parallel for + for ( size_t i = 0 ; i < m_nnp ; ++i ) + out(m_iip(i)) = m_data_pp(i); + + // return output + return out; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace ... + +// ================================================================================================= + +#endif diff --git a/include/xGooseFEM/MatrixPartitioned.h b/include/xGooseFEM/MatrixPartitioned.h new file mode 100644 index 0000000..c7d7644 --- /dev/null +++ b/include/xGooseFEM/MatrixPartitioned.h @@ -0,0 +1,97 @@ +/* ================================================================================================= + +(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM + +================================================================================================= */ + +#ifndef XGOOSEFEM_MATRIXPARTITIONED_H +#define XGOOSEFEM_MATRIXPARTITIONED_H + +// ------------------------------------------------------------------------------------------------- + +#include "GooseFEM.h" + +// =========================================== GooseFEM ============================================ + +namespace xGooseFEM { + +// ------------------------------------------------------------------------------------------------- + +class MatrixPartitioned +{ +public: + + // default constructor + MatrixPartitioned() = default; + + // constructor + MatrixPartitioned(const xt::xtensor &conn, const xt::xtensor &dofs, + const xt::xtensor &iip); + + // dimensions + size_t nelem() const; // number of elements + size_t nne() const; // number of nodes per element + size_t nnode() const; // number of nodes + size_t ndim() const; // number of dimensions + size_t ndof() const; // number of DOFs + size_t nnu() const; // number of unknown DOFs + size_t nnp() const; // number of prescribed DOFs + + // DOF lists + xt::xtensor dofs() const; // DOFs + xt::xtensor iiu() const; // unknown DOFs + xt::xtensor iip() const; // prescribed DOFs + + // assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] + void assemble(const xt::xtensor &elemmat); + + // solve: x_u = A_uu \ ( b_u - A_up * x_p ) + xt::xtensor solve(const xt::xtensor &b_u, const xt::xtensor &x_p); + +private: + + // the matrix + Eigen::SparseMatrix m_data_uu; + Eigen::SparseMatrix m_data_up; + Eigen::SparseMatrix m_data_pu; + Eigen::SparseMatrix m_data_pp; + + // the matrix to assemble + std::vector m_trip_uu; + std::vector m_trip_up; + std::vector m_trip_pu; + std::vector m_trip_pp; + + // solver (re-used to solve different RHS) + Eigen::SimplicialLDLT> m_solver; + + // signal changes to data compare to the last inverse + bool m_change=false; + + // bookkeeping + xt::xtensor m_conn; // connectivity [nelem, nne ] + xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] + xt::xtensor m_part; // DOF-numbers per node, renumbered [nnode, ndim] + xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] + xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] + + // dimensions + size_t m_nelem; // number of elements + size_t m_nne; // number of nodes per element + size_t m_nnode; // number of nodes + size_t m_ndim; // number of dimensions + size_t m_ndof; // number of DOFs + size_t m_nnu; // number of unknown DOFs + size_t m_nnp; // number of prescribed DOFs + + // compute inverse (automatically evaluated by "solve") + void factorize(); +}; + +// ------------------------------------------------------------------------------------------------- + +} // namespace ... + +// ================================================================================================= + +#endif diff --git a/include/xGooseFEM/MatrixPartitioned.hpp b/include/xGooseFEM/MatrixPartitioned.hpp new file mode 100644 index 0000000..880e2ad --- /dev/null +++ b/include/xGooseFEM/MatrixPartitioned.hpp @@ -0,0 +1,238 @@ +/* ================================================================================================= + +(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM + +================================================================================================= */ + +#ifndef XGOOSEFEM_MATRIXPARTITIONED_CPP +#define XGOOSEFEM_MATRIXPARTITIONED_CPP + +// ------------------------------------------------------------------------------------------------- + +#include "MatrixPartitioned.h" + +#include + +// ================================================================================================= + +namespace xGooseFEM { + +// ------------------------------------------------------------------------------------------------- + +inline MatrixPartitioned::MatrixPartitioned(const xt::xtensor &conn, + const xt::xtensor &dofs, const xt::xtensor &iip) : + m_conn(conn), m_dofs(dofs), m_iip(iip) +{ + // mesh dimensions + m_nelem = m_conn.shape()[0]; + m_nne = m_conn.shape()[1]; + m_nnode = m_dofs.shape()[0]; + m_ndim = m_dofs.shape()[1]; + + // list with unknown DOFs + m_iiu = xt::setdiff1d(dofs, iip); + + // dimensions of the system + m_ndof = xt::amax(m_dofs)[0] + 1; + m_nnp = m_iip.size(); + m_nnu = m_iiu.size(); + + // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) + m_part = Mesh::reorder(m_dofs, m_iip, "end"); + + // allocate triplet list + m_trip_uu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); + m_trip_up.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); + m_trip_pu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); + m_trip_pp.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim); + + // allocate sparse matrices + m_data_uu.resize(m_nnu,m_nnu); + m_data_up.resize(m_nnu,m_nnp); + m_data_pu.resize(m_nnp,m_nnu); + m_data_pp.resize(m_nnp,m_nnp); + + // check consistency + assert( xt::amax(m_conn)[0] + 1 == m_nnode ); + assert( xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0] ); + assert( m_ndof <= m_nnode * m_ndim ); +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::nelem() const +{ + return m_nelem; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::nne() const +{ + return m_nne; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::nnode() const +{ + return m_nnode; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::ndim() const +{ + return m_ndim; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::ndof() const +{ + return m_ndof; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::nnu() const +{ + return m_nnu; +} + +// ------------------------------------------------------------------------------------------------- + +inline size_t MatrixPartitioned::nnp() const +{ + return m_nnp; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixPartitioned::dofs() const +{ + return m_dofs; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixPartitioned::iiu() const +{ + return m_iiu; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixPartitioned::iip() const +{ + return m_iip; +} + +// ------------------------------------------------------------------------------------------------- + +inline void MatrixPartitioned::assemble(const xt::xtensor &elemmat) +{ + // check input + assert( elemmat.shape()[0] == m_nelem ); + assert( elemmat.shape()[1] == m_nne*m_ndim ); + assert( elemmat.shape()[2] == m_nne*m_ndim ); + + // initialize triplet list + m_trip_uu.clear(); + m_trip_up.clear(); + m_trip_pu.clear(); + m_trip_pp.clear(); + + // assemble to triplet list + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + for ( size_t m = 0 ; m < m_nne ; ++m ) + { + for ( size_t i = 0 ; i < m_ndim ; ++i ) + { + size_t di = m_part(m_conn(e,m),i); + + for ( size_t n = 0 ; n < m_nne ; ++n ) + { + for ( size_t j = 0 ; j < m_ndim ; ++j ) + { + size_t dj = m_part(m_conn(e,n),j); + + if ( di < m_nnu and dj < m_nnu ) + m_trip_uu.push_back(TripD(di ,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); + else if ( di < m_nnu ) + m_trip_up.push_back(TripD(di ,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); + else if ( dj < m_nnu ) + m_trip_pu.push_back(TripD(di-m_nnu,dj ,elemmat(e,m*m_ndim+i,n*m_ndim+j))); + else + m_trip_pp.push_back(TripD(di-m_nnu,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j))); + } + } + } + } + } + + // convert to sparse matrix + m_data_uu.setFromTriplets(m_trip_uu.begin(), m_trip_uu.end()); + m_data_up.setFromTriplets(m_trip_up.begin(), m_trip_up.end()); + m_data_pu.setFromTriplets(m_trip_pu.begin(), m_trip_pu.end()); + m_data_pp.setFromTriplets(m_trip_pp.begin(), m_trip_pp.end()); + + // signal change + m_change = true; +} + +// ------------------------------------------------------------------------------------------------- + +inline void MatrixPartitioned::factorize() +{ + // skip for unchanged "m_data" + if ( ! m_change ) return; + + // factorise + m_solver.compute(m_data_uu); + + // reset signal + m_change = false; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor MatrixPartitioned::solve( + const xt::xtensor &b_u, const xt::xtensor &x_p) +{ + // check input + assert( b_u.shape()[0] == m_nnu ); + assert( x_p.shape()[0] == m_nnp ); + + // factorise (if needed) + this->factorize(); + + // copy of input as Eigen objects + // - allocate + Eigen::VectorXd B_u(m_nnu,1); + Eigen::VectorXd X_p(m_nnp,1); + // - copy + std::copy(b_u.begin(), b_u.end(), B_u.data()); + std::copy(x_p.begin(), x_p.end(), X_p.data()); + + // solve + Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_data_up * X_p)); + + // copy of output from Eigen object + // - allocate + xt::xtensor x_u = xt::empty({m_nnu}); + // - copy + std::copy(X_u.data(), X_u.data()+m_nnu, x_u.begin()); + + // return output + return x_u; +} + +// ------------------------------------------------------------------------------------------------- + +} // namespace ... + +// ================================================================================================= + +#endif diff --git a/include/xGooseFEM/Vector.h b/include/xGooseFEM/Vector.h index 8849f94..2e85114 100644 --- a/include/xGooseFEM/Vector.h +++ b/include/xGooseFEM/Vector.h @@ -1,192 +1,126 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_VECTOR_H #define XGOOSEFEM_VECTOR_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =========================================== GooseFEM ============================================ namespace xGooseFEM { // ------------------------------------------------------------------------------------------------- /* - "nodevec" - nodal vectors - [nnode, ndim] - "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] - "dofval" - DOF values - [ndof] - "dofval_u" - DOF values (Unknown) "== dofval[iiu]" - [nnu] - "dofval_p" - DOF values (Prescribed) "== dofval[iiu]" - [nnp] + "nodevec" - nodal vectors - [nnode, ndim] + "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] + "dofval" - DOF values - [ndof] */ class Vector { public: // constructor Vector() = default; Vector(const xt::xtensor &conn, const xt::xtensor &dofs); - Vector(const xt::xtensor &conn, const xt::xtensor &dofs, - const xt::xtensor &iip); - // dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs - size_t nnu() const; // number of unknown DOFs - size_t nnp() const; // number of prescribed DOFs // DOF lists xt::xtensor dofs() const; // DOFs - xt::xtensor iiu() const; // unknown DOFs - xt::xtensor iip() const; // prescribed DOFs // convert to "dofval" (overwrite entries that occur more than once) -- (auto allocation below) - void asDofs (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, - xt::xtensor &dofval) const; - - void asDofs (const xt::xtensor &nodevec, + void asDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const; - void asDofs (const xt::xtensor &elemvec, + void asDofs(const xt::xtensor &elemvec, xt::xtensor &dofval) const; - void asDofs_u(const xt::xtensor &nodevec, - xt::xtensor &dofval_u) const; - - void asDofs_u(const xt::xtensor &elemvec, - xt::xtensor &dofval_u) const; - - void asDofs_p(const xt::xtensor &nodevec, - xt::xtensor &dofval_p) const; - - void asDofs_p(const xt::xtensor &elemvec, - xt::xtensor &dofval_p) const; - // convert to "nodevec" (overwrite entries that occur more than once) -- (auto allocation below) - void asNode(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, - xt::xtensor &nodevec) const; - void asNode(const xt::xtensor &dofval, xt::xtensor &nodevec) const; - void asNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const; // convert to "elemvec" (overwrite entries that occur more than once) -- (auto allocation below) - void asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, - xt::xtensor &elemvec) const; - void asElement(const xt::xtensor &dofval, xt::xtensor &elemvec) const; void asElement(const xt::xtensor &nodevec, xt::xtensor &elemvec) const; - // convert to "dofval" (overwrite entries that occur more than once) - - xt::xtensor asDofs (const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const; - - xt::xtensor asDofs (const xt::xtensor &nodevec) const; + // assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) - xt::xtensor asDofs (const xt::xtensor &elemvec) const; + void assembleDofs(const xt::xtensor &nodevec, + xt::xtensor &dofval) const; - xt::xtensor asDofs_u(const xt::xtensor &nodevec) const; + void assembleDofs(const xt::xtensor &elemvec, + xt::xtensor &dofval) const; - xt::xtensor asDofs_u(const xt::xtensor &elemvec) const; + // assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) - xt::xtensor asDofs_p(const xt::xtensor &nodevec) const; + void assembleNode(const xt::xtensor &elemvec, + xt::xtensor &nodevec) const; - xt::xtensor asDofs_p(const xt::xtensor &elemvec) const; + // auto allocation of the functions above - // convert to "nodevec" (overwrite entries that occur more than once) + xt::xtensor asDofs(const xt::xtensor &nodevec) const; - xt::xtensor asNode(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const; + xt::xtensor asDofs(const xt::xtensor &elemvec) const; xt::xtensor asNode(const xt::xtensor &dofval) const; xt::xtensor asNode(const xt::xtensor &elemvec) const; - // convert to "elemvec" (overwrite entries that occur more than once) - - xt::xtensor asElement(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const; - xt::xtensor asElement(const xt::xtensor &dofval) const; xt::xtensor asElement(const xt::xtensor &nodevec) const; - // assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) - - void assembleDofs (const xt::xtensor &nodevec, xt::xtensor &dofval ) const; - void assembleDofs (const xt::xtensor &elemvec, xt::xtensor &dofval ) const; - void assembleDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval_u) const; - void assembleDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval_u) const; - void assembleDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval_p) const; - void assembleDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval_p) const; - - // assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) - - void assembleNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const; + xt::xtensor assembleDofs(const xt::xtensor &nodevec) const; - // assemble "dofval" (adds entries that occur more that once) - - xt::xtensor assembleDofs (const xt::xtensor &nodevec) const; - xt::xtensor assembleDofs (const xt::xtensor &elemvec) const; - xt::xtensor assembleDofs_u(const xt::xtensor &nodevec) const; - xt::xtensor assembleDofs_u(const xt::xtensor &elemvec) const; - xt::xtensor assembleDofs_p(const xt::xtensor &nodevec) const; - xt::xtensor assembleDofs_p(const xt::xtensor &elemvec) const; - - // assemble "nodevec" (adds entries that occur more that once) + xt::xtensor assembleDofs(const xt::xtensor &elemvec) const; xt::xtensor assembleNode(const xt::xtensor &elemvec) const; private: // bookkeeping - xt::xtensor m_conn; // connectivity [nelem, nne ] - xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] - xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] - xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] - - // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) - xt::xtensor m_part; + xt::xtensor m_conn; // connectivity [nelem, nne ] + xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] // dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs - size_t m_nnu; // number of unknown DOFs - size_t m_nnp; // number of prescribed DOFs }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/Vector.hpp b/include/xGooseFEM/Vector.hpp index 4e3ad9a..e99bc65 100644 --- a/include/xGooseFEM/Vector.hpp +++ b/include/xGooseFEM/Vector.hpp @@ -1,723 +1,337 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_VECTOR_CPP #define XGOOSEFEM_VECTOR_CPP // ------------------------------------------------------------------------------------------------- #include "Vector.h" -// =========================================== GooseFEM ============================================ +// ================================================================================================= namespace xGooseFEM { -// ------------------------------------------ constructor ------------------------------------------ +// ------------------------------------------------------------------------------------------------- inline Vector::Vector(const xt::xtensor &conn, const xt::xtensor &dofs) : m_conn(conn), m_dofs(dofs) { // mesh dimensions m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; - // list with prescribed and unknown DOFs - m_iip = xt::empty({0}); - m_iiu = xt::unique(dofs); - - // dimensions of the system - m_ndof = xt::amax(m_dofs)[0] + 1; - m_nnp = m_iip.size(); - m_nnu = m_iiu.size(); - - // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) - m_part = Mesh::reorder(m_dofs, m_iip, "end"); - - // check consistency - assert( xt::amax(m_conn)[0] + 1 == m_nnode ); - assert( m_ndof <= m_nnode * m_ndim ); -} - -// ------------------------------------------ constructor ------------------------------------------ - -inline Vector::Vector(const xt::xtensor &conn, const xt::xtensor &dofs, - const xt::xtensor &iip) : m_conn(conn), m_dofs(dofs), m_iip(iip) -{ - // mesh dimensions - m_nelem = m_conn.shape()[0]; - m_nne = m_conn.shape()[1]; - m_nnode = m_dofs.shape()[0]; - m_ndim = m_dofs.shape()[1]; - - // list with unknown DOFs - m_iiu = xt::setdiff1d(dofs, iip); - // dimensions of the system m_ndof = xt::amax(m_dofs)[0] + 1; - m_nnp = m_iip.size(); - m_nnu = m_iiu.size(); - - // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) - m_part = Mesh::reorder(m_dofs, m_iip, "end"); // check consistency - assert( xt::amax(m_conn)[0] + 1 == m_nnode ); - assert( xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0] ); - assert( m_ndof <= m_nnode * m_ndim ); + assert( xt::amax(m_conn)[0] + 1 == m_nnode ); + assert( m_ndof <= m_nnode * m_ndim ); } -// -------------------------------------- number of elements --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t Vector::nelem() const { return m_nelem; } -// ---------------------------------- number of nodes per element ---------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t Vector::nne() const { return m_nne; } -// ---------------------------------------- number of nodes ---------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t Vector::nnode() const { return m_nnode; } -// ------------------------------------- number of dimensions -------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t Vector::ndim() const { return m_ndim; } -// ---------------------------------------- number of DOFs ----------------------------------------- +// ------------------------------------------------------------------------------------------------- inline size_t Vector::ndof() const { return m_ndof; } -// ------------------------------------ number of unknown DOFs ------------------------------------- - -inline size_t Vector::nnu() const -{ - return m_nnu; -} - -// ----------------------------------- number of prescribed DOFs ----------------------------------- - -inline size_t Vector::nnp() const -{ - return m_nnp; -} - -// ------------------------------------------ return DOFs ------------------------------------------ +// ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::dofs() const { return m_dofs; } -// -------------------------------------- return unknown DOFs -------------------------------------- - -inline xt::xtensor Vector::iiu() const -{ - return m_iiu; -} - -// ------------------------------------ return prescribed DOFs ------------------------------------- - -inline xt::xtensor Vector::iip() const -{ - return m_iip; -} - -// ------------------- assemble: dofval[iiu] = dofval_u; dofval[iip] = dofval_p -------------------- - -inline void Vector::asDofs(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p, xt::xtensor &dofval) const -{ - assert( dofval_u.size() == m_nnu ); - assert( dofval_p.size() == m_nnp ); - assert( dofval.size() == m_ndof ); - - #pragma omp parallel for - for ( size_t i = 0 ; i < m_nnu ; ++i ) dofval(m_iiu(i)) = dofval_u(i); - - #pragma omp parallel for - for ( size_t i = 0 ; i < m_nnp ; ++i ) dofval(m_iip(i)) = dofval_p(i); -} - -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->asDofs(dofval_u, dofval_p, dofval); - - return dofval; -} - -// --------------------------------------- nodevec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::asDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(n,i)) = nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->asDofs(nodevec, dofval); - - return dofval; -} - -// -------------------------------------- nodevec -> dofval_u -------------------------------------- - -inline void Vector::asDofs_u(const xt::xtensor &nodevec, - xt::xtensor &dofval_u) const -{ - assert( nodevec.shape()[0] == m_nnode ); - assert( nodevec.shape()[1] == m_ndim ); - assert( dofval_u.size() == m_nnu ); - - #pragma omp parallel for - for ( size_t n = 0 ; n < m_nnode ; ++n ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(n,i) < m_nnu ) - dofval_u(m_part(n,i)) = nodevec(n,i); -} - -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_u(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval_u = xt::empty({m_nnu}); - - this->asDofs_u(nodevec, dofval_u); - - return dofval_u; -} - -// -------------------------------------- nodevec -> dofval_p -------------------------------------- - -inline void Vector::asDofs_p(const xt::xtensor &nodevec, - xt::xtensor &dofval_p) const -{ - assert( nodevec.shape()[0] == m_nnode ); - assert( nodevec.shape()[1] == m_ndim ); - assert( dofval_p.size() == m_nnp ); - - #pragma omp parallel for - for ( size_t n = 0 ; n < m_nnode ; ++n ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(n,i) >= m_nnu ) - dofval_p(m_part(n,i)-m_nnu) = nodevec(n,i); -} - -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_p(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval_p = xt::empty({m_nnp}); - - this->asDofs_p(nodevec, dofval_p); - - return dofval_p; -} - -// --------------------------------------- elemvec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::asDofs(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(m_conn(e,m),i)) = elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->asDofs(elemvec, dofval); - - return dofval; -} - -// -------------------------------------- elemvec -> dofval_u -------------------------------------- - -inline void Vector::asDofs_u(const xt::xtensor &elemvec, - xt::xtensor &dofval_u) const -{ - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - assert( dofval_u.size() == m_nnu ); - - #pragma omp parallel for - for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(m_conn(e,m),i) < m_nnu ) - dofval_u(m_part(m_conn(e,m),i)) = elemvec(e,m,i); -} - -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_u(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval_u = xt::empty({m_nnu}); - - this->asDofs_u(elemvec, dofval_u); - - return dofval_u; -} - -// -------------------------------------- elemvec -> dofval_p -------------------------------------- - -inline void Vector::asDofs_p(const xt::xtensor &elemvec, - xt::xtensor &dofval_p) const -{ - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - assert( dofval_p.size() == m_nnp ); - - #pragma omp parallel for - for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(m_conn(e,m),i) >= m_nnu ) - dofval_p(m_part(m_conn(e,m),i)-m_nnu) = elemvec(e,m,i); -} - -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_p(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval_p = xt::empty({m_nnp}); - - this->asDofs_p(elemvec, dofval_p); - - return dofval_p; -} - -// --------------------------------------- dofval -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::asNode(const xt::xtensor &dofval, xt::xtensor &nodevec) const { assert( dofval.size() == m_ndof ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(n,i) = dofval(m_dofs(n,i)); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asNode(const xt::xtensor &dofval) const -{ - xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - - this->asNode(dofval, nodevec); - - return nodevec; -} - -// --------------------------------------- dofval -> nodevec --------------------------------------- - -inline void Vector::asNode(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p, xt::xtensor &nodevec) const -{ - assert( dofval_u.size() == m_nnu ); - assert( dofval_p.size() == m_nnp ); - assert( nodevec.shape()[0] == m_nnode ); - assert( nodevec.shape()[1] == m_ndim ); - - #pragma omp parallel for - for ( size_t n = 0 ; n < m_nnode ; ++n ) { - for ( size_t i = 0 ; i < m_ndim ; ++i ) { - if ( m_part(n,i) < m_nnu ) nodevec(n,i) = dofval_u(m_part(n,i) ); - else nodevec(n,i) = dofval_p(m_part(n,i)-m_nnu); - } - } -} - -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asNode(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const -{ - xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - - this->asNode(dofval_u, dofval_p, nodevec); - - return nodevec; -} - -// -------------------------------------- elemvec -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::asNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(m_conn(e,m),i) = elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asNode(const xt::xtensor &elemvec) const -{ - xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - - this->asNode(elemvec, nodevec); - - return nodevec; -} - -// --------------------------------------- dofval -> elemvec --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::asElement(const xt::xtensor &dofval, xt::xtensor &elemvec) const { assert( dofval.size() == m_ndof ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) elemvec(e,m,i) = dofval(m_dofs(m_conn(e,m),i)); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asElement(const xt::xtensor &dofval) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - - this->asElement(dofval, elemvec); - - return elemvec; -} - -// --------------------------------------- dofval -> elemvec --------------------------------------- - -inline void Vector::asElement(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p, xt::xtensor &elemvec) const -{ - assert( dofval_u.size() == m_nnu ); - assert( dofval_p.size() == m_nnp ); - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - - #pragma omp parallel for - for ( size_t e = 0 ; e < m_nelem ; ++e ) { - for ( size_t m = 0 ; m < m_nne ; ++m ) { - for ( size_t i = 0 ; i < m_ndim ; ++i ) { - if ( m_part(m_conn(e,m),i) Vector::asElement(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - - this->asElement(dofval_u, dofval_p, elemvec); - - return elemvec; -} - -// -------------------------------------- nodevec -> elemvec --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::asElement(const xt::xtensor &nodevec, xt::xtensor &elemvec) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) elemvec(e,m,i) = nodevec(m_conn(e,m),i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asElement(const xt::xtensor &nodevec) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - - this->asElement(nodevec, elemvec); - - return elemvec; -} - -// --------------------------------------- nodevec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- inline void Vector::assembleDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_ndof ); dofval.fill(0.0); for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(n,i)) += nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->assembleDofs(nodevec, dofval); - - return dofval; -} - -// -------------------------------------- nodevec -> dofval_u -------------------------------------- - -inline void Vector::assembleDofs_u(const xt::xtensor &nodevec, - xt::xtensor &dofval_u) const -{ - assert( nodevec.shape()[0] == m_nnode ); - assert( nodevec.shape()[1] == m_ndim ); - assert( dofval_u.size() == m_nnu ); - - dofval_u.fill(0.0); - - for ( size_t n = 0 ; n < m_nnode ; ++n ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(n,i) < m_nnu ) - dofval_u(m_part(n,i)) += nodevec(n,i); -} - -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::assembleDofs_u(const xt::xtensor &nodevec) const +inline void Vector::assembleDofs(const xt::xtensor &elemvec, + xt::xtensor &dofval) const { - xt::xtensor dofval_u = xt::empty({m_nnu}); + assert( elemvec.shape()[0] == m_nelem ); + assert( elemvec.shape()[1] == m_nne ); + assert( elemvec.shape()[2] == m_ndim ); + assert( dofval.size() == m_ndof ); - this->assembleDofs_u(nodevec, dofval_u); + dofval.fill(0.0); - return dofval_u; + for ( size_t e = 0 ; e < m_nelem ; ++e ) + for ( size_t m = 0 ; m < m_nne ; ++m ) + for ( size_t i = 0 ; i < m_ndim ; ++i ) + dofval(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); } -// -------------------------------------- nodevec -> dofval_p -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_p(const xt::xtensor &nodevec, - xt::xtensor &dofval_p) const +inline void Vector::assembleNode(const xt::xtensor &elemvec, + xt::xtensor &nodevec) const { + assert( elemvec.shape()[0] == m_nelem ); + assert( elemvec.shape()[1] == m_nne ); + assert( elemvec.shape()[2] == m_ndim ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); - assert( dofval_p.size() == m_nnp ); - dofval_p.fill(0.0); + nodevec.fill(0.0); - for ( size_t n = 0 ; n < m_nnode ; ++n ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(n,i) >= m_nnu ) - dofval_p(m_part(n,i)-m_nnu) += nodevec(n,i); + for ( size_t e = 0 ; e < m_nelem ; ++e ) + for ( size_t m = 0 ; m < m_nne ; ++m ) + for ( size_t i = 0 ; i < m_ndim ; ++i ) + nodevec(m_conn(e,m),i) += elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::assembleDofs_p(const xt::xtensor &nodevec) const +inline xt::xtensor Vector::asDofs(const xt::xtensor &nodevec) const { - xt::xtensor dofval_p = xt::empty({m_nnp}); + xt::xtensor dofval = xt::empty({m_ndof}); - this->assembleDofs_p(nodevec, dofval_p); + this->asDofs(nodevec, dofval); - return dofval_p; + return dofval; } -// --------------------------------------- elemvec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs(const xt::xtensor &elemvec, - xt::xtensor &dofval) const +inline xt::xtensor Vector::asDofs(const xt::xtensor &elemvec) const { - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - assert( dofval.size() == m_ndof ); + xt::xtensor dofval = xt::empty({m_ndof}); - dofval.fill(0.0); + this->asDofs(elemvec, dofval); - for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - dofval(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); + return dofval; } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::assembleDofs(const xt::xtensor &elemvec) const +inline xt::xtensor Vector::asNode(const xt::xtensor &dofval) const { - xt::xtensor dofval = xt::empty({m_ndof}); + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - this->assembleDofs(elemvec, dofval); + this->asNode(dofval, nodevec); - return dofval; + return nodevec; } -// -------------------------------------- elemvec -> dofval_u -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_u(const xt::xtensor &elemvec, - xt::xtensor &dofval_u) const +inline xt::xtensor Vector::asNode(const xt::xtensor &elemvec) const { - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - assert( dofval_u.size() == m_nnu ); + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - dofval_u.fill(0.0); + this->asNode(elemvec, nodevec); - for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(m_conn(e,m),i) < m_nnu ) - dofval_u(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); + return nodevec; } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::assembleDofs_u(const xt::xtensor &elemvec) const +inline xt::xtensor Vector::asElement(const xt::xtensor &dofval) const { - xt::xtensor dofval_u = xt::empty({m_nnu}); + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - this->assembleDofs_u(elemvec, dofval_u); + this->asElement(dofval, elemvec); - return dofval_u; + return elemvec; } -// -------------------------------------- elemvec -> dofval_p -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_p(const xt::xtensor &elemvec, - xt::xtensor &dofval_p) const +inline xt::xtensor Vector::asElement(const xt::xtensor &nodevec) const { - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - assert( dofval_p.size() == m_nnp ); + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - dofval_p.fill(0.0); + this->asElement(nodevec, elemvec); - for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - if ( m_part(m_conn(e,m),i) >= m_nnu ) - dofval_p(m_dofs(m_conn(e,m),i)-m_nnu) += elemvec(e,m,i); + return elemvec; } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::assembleDofs_p(const xt::xtensor &elemvec) const +inline xt::xtensor Vector::assembleDofs(const xt::xtensor &nodevec) const { - xt::xtensor dofval_p = xt::empty({m_nnp}); + xt::xtensor dofval = xt::empty({m_ndof}); - this->assembleDofs_p(elemvec, dofval_p); + this->assembleDofs(nodevec, dofval); - return dofval_p; + return dofval; } -// -------------------------------------- elemvec -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleNode(const xt::xtensor &elemvec, - xt::xtensor &nodevec) const +inline xt::xtensor Vector::assembleDofs(const xt::xtensor &elemvec) const { - assert( elemvec.shape()[0] == m_nelem ); - assert( elemvec.shape()[1] == m_nne ); - assert( elemvec.shape()[2] == m_ndim ); - assert( nodevec.shape()[0] == m_nnode ); - assert( nodevec.shape()[1] == m_ndim ); + xt::xtensor dofval = xt::empty({m_ndof}); - nodevec.fill(0.0); + this->assembleDofs(elemvec, dofval); - for ( size_t e = 0 ; e < m_nelem ; ++e ) - for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - nodevec(m_conn(e,m),i) += elemvec(e,m,i); + return dofval; } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleNode(const xt::xtensor &elemvec) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->assembleNode(elemvec, nodevec); return nodevec; } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/Vector.h b/include/xGooseFEM/VectorPartitioned.h similarity index 81% copy from include/xGooseFEM/Vector.h copy to include/xGooseFEM/VectorPartitioned.h index 8849f94..540a242 100644 --- a/include/xGooseFEM/Vector.h +++ b/include/xGooseFEM/VectorPartitioned.h @@ -1,192 +1,194 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef XGOOSEFEM_VECTOR_H -#define XGOOSEFEM_VECTOR_H +#ifndef XGOOSEFEM_VECTORPARTITIONED_H +#define XGOOSEFEM_VECTORPARTITIONED_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =========================================== GooseFEM ============================================ namespace xGooseFEM { // ------------------------------------------------------------------------------------------------- /* "nodevec" - nodal vectors - [nnode, ndim] "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] "dofval" - DOF values - [ndof] "dofval_u" - DOF values (Unknown) "== dofval[iiu]" - [nnu] "dofval_p" - DOF values (Prescribed) "== dofval[iiu]" - [nnp] */ -class Vector +class VectorPartitioned { public: // constructor - Vector() = default; + VectorPartitioned() = default; - Vector(const xt::xtensor &conn, const xt::xtensor &dofs); - - Vector(const xt::xtensor &conn, const xt::xtensor &dofs, + VectorPartitioned(const xt::xtensor &conn, const xt::xtensor &dofs, const xt::xtensor &iip); // dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs size_t nnu() const; // number of unknown DOFs size_t nnp() const; // number of prescribed DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // unknown DOFs xt::xtensor iip() const; // prescribed DOFs // convert to "dofval" (overwrite entries that occur more than once) -- (auto allocation below) void asDofs (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &dofval) const; void asDofs (const xt::xtensor &nodevec, xt::xtensor &dofval) const; void asDofs (const xt::xtensor &elemvec, xt::xtensor &dofval) const; void asDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval_u) const; void asDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval_u) const; void asDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval_p) const; void asDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval_p) const; // convert to "nodevec" (overwrite entries that occur more than once) -- (auto allocation below) void asNode(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &nodevec) const; void asNode(const xt::xtensor &dofval, xt::xtensor &nodevec) const; - void asNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const; // convert to "elemvec" (overwrite entries that occur more than once) -- (auto allocation below) void asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &elemvec) const; void asElement(const xt::xtensor &dofval, xt::xtensor &elemvec) const; void asElement(const xt::xtensor &nodevec, xt::xtensor &elemvec) const; - // convert to "dofval" (overwrite entries that occur more than once) + // assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) + + void assembleDofs (const xt::xtensor &nodevec, + xt::xtensor &dofval ) const; + + void assembleDofs (const xt::xtensor &elemvec, + xt::xtensor &dofval ) const; + + void assembleDofs_u(const xt::xtensor &nodevec, + xt::xtensor &dofval_u) const; + + void assembleDofs_u(const xt::xtensor &elemvec, + xt::xtensor &dofval_u) const; + + void assembleDofs_p(const xt::xtensor &nodevec, + xt::xtensor &dofval_p) const; + + void assembleDofs_p(const xt::xtensor &elemvec, + xt::xtensor &dofval_p) const; + + // assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) + + void assembleNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const; + + // auto allocation of the functions above - xt::xtensor asDofs (const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const; + xt::xtensor asDofs (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const; xt::xtensor asDofs (const xt::xtensor &nodevec) const; xt::xtensor asDofs (const xt::xtensor &elemvec) const; xt::xtensor asDofs_u(const xt::xtensor &nodevec) const; xt::xtensor asDofs_u(const xt::xtensor &elemvec) const; xt::xtensor asDofs_p(const xt::xtensor &nodevec) const; xt::xtensor asDofs_p(const xt::xtensor &elemvec) const; - // convert to "nodevec" (overwrite entries that occur more than once) - - xt::xtensor asNode(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const; + xt::xtensor asNode(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const; xt::xtensor asNode(const xt::xtensor &dofval) const; xt::xtensor asNode(const xt::xtensor &elemvec) const; - // convert to "elemvec" (overwrite entries that occur more than once) - - xt::xtensor asElement(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const; + xt::xtensor asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const; xt::xtensor asElement(const xt::xtensor &dofval) const; xt::xtensor asElement(const xt::xtensor &nodevec) const; - // assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) - - void assembleDofs (const xt::xtensor &nodevec, xt::xtensor &dofval ) const; - void assembleDofs (const xt::xtensor &elemvec, xt::xtensor &dofval ) const; - void assembleDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval_u) const; - void assembleDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval_u) const; - void assembleDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval_p) const; - void assembleDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval_p) const; - - // assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) - - void assembleNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const; - - // assemble "dofval" (adds entries that occur more that once) - xt::xtensor assembleDofs (const xt::xtensor &nodevec) const; + xt::xtensor assembleDofs (const xt::xtensor &elemvec) const; + xt::xtensor assembleDofs_u(const xt::xtensor &nodevec) const; + xt::xtensor assembleDofs_u(const xt::xtensor &elemvec) const; + xt::xtensor assembleDofs_p(const xt::xtensor &nodevec) const; - xt::xtensor assembleDofs_p(const xt::xtensor &elemvec) const; - // assemble "nodevec" (adds entries that occur more that once) + xt::xtensor assembleDofs_p(const xt::xtensor &elemvec) const; xt::xtensor assembleNode(const xt::xtensor &elemvec) const; private: // bookkeeping xt::xtensor m_conn; // connectivity [nelem, nne ] xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) xt::xtensor m_part; // dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs size_t m_nnu; // number of unknown DOFs size_t m_nnp; // number of prescribed DOFs }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/Vector.hpp b/include/xGooseFEM/VectorPartitioned.hpp similarity index 64% copy from include/xGooseFEM/Vector.hpp copy to include/xGooseFEM/VectorPartitioned.hpp index 4e3ad9a..c0fe780 100644 --- a/include/xGooseFEM/Vector.hpp +++ b/include/xGooseFEM/VectorPartitioned.hpp @@ -1,723 +1,703 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef XGOOSEFEM_VECTOR_CPP -#define XGOOSEFEM_VECTOR_CPP +#ifndef XGOOSEFEM_VECTORPARTITIONED_CPP +#define XGOOSEFEM_VECTORPARTITIONED_CPP // ------------------------------------------------------------------------------------------------- -#include "Vector.h" +#include "VectorPartitioned.h" -// =========================================== GooseFEM ============================================ +// ================================================================================================= namespace xGooseFEM { -// ------------------------------------------ constructor ------------------------------------------ - -inline Vector::Vector(const xt::xtensor &conn, const xt::xtensor &dofs) : - m_conn(conn), m_dofs(dofs) -{ - // mesh dimensions - m_nelem = m_conn.shape()[0]; - m_nne = m_conn.shape()[1]; - m_nnode = m_dofs.shape()[0]; - m_ndim = m_dofs.shape()[1]; - - // list with prescribed and unknown DOFs - m_iip = xt::empty({0}); - m_iiu = xt::unique(dofs); - - // dimensions of the system - m_ndof = xt::amax(m_dofs)[0] + 1; - m_nnp = m_iip.size(); - m_nnu = m_iiu.size(); - - // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) - m_part = Mesh::reorder(m_dofs, m_iip, "end"); - - // check consistency - assert( xt::amax(m_conn)[0] + 1 == m_nnode ); - assert( m_ndof <= m_nnode * m_ndim ); -} - -// ------------------------------------------ constructor ------------------------------------------ +// ------------------------------------------------------------------------------------------------- -inline Vector::Vector(const xt::xtensor &conn, const xt::xtensor &dofs, - const xt::xtensor &iip) : m_conn(conn), m_dofs(dofs), m_iip(iip) +inline VectorPartitioned::VectorPartitioned(const xt::xtensor &conn, + const xt::xtensor &dofs, const xt::xtensor &iip) : + m_conn(conn), m_dofs(dofs), m_iip(iip) { // mesh dimensions m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; // list with unknown DOFs m_iiu = xt::setdiff1d(dofs, iip); // dimensions of the system m_ndof = xt::amax(m_dofs)[0] + 1; m_nnp = m_iip.size(); m_nnu = m_iiu.size(); // DOFs per node, such that iiu = arange(nnu), iip = nnu + arange(nnp) m_part = Mesh::reorder(m_dofs, m_iip, "end"); // check consistency - assert( xt::amax(m_conn)[0] + 1 == m_nnode ); - assert( xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0] ); - assert( m_ndof <= m_nnode * m_ndim ); + assert( xt::amax(m_conn)[0] + 1 == m_nnode ); + assert( xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0] ); + assert( m_ndof <= m_nnode * m_ndim ); } -// -------------------------------------- number of elements --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::nelem() const +inline size_t VectorPartitioned::nelem() const { return m_nelem; } -// ---------------------------------- number of nodes per element ---------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::nne() const +inline size_t VectorPartitioned::nne() const { return m_nne; } -// ---------------------------------------- number of nodes ---------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::nnode() const +inline size_t VectorPartitioned::nnode() const { return m_nnode; } -// ------------------------------------- number of dimensions -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::ndim() const +inline size_t VectorPartitioned::ndim() const { return m_ndim; } -// ---------------------------------------- number of DOFs ----------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::ndof() const +inline size_t VectorPartitioned::ndof() const { return m_ndof; } -// ------------------------------------ number of unknown DOFs ------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::nnu() const +inline size_t VectorPartitioned::nnu() const { return m_nnu; } -// ----------------------------------- number of prescribed DOFs ----------------------------------- +// ------------------------------------------------------------------------------------------------- -inline size_t Vector::nnp() const +inline size_t VectorPartitioned::nnp() const { return m_nnp; } -// ------------------------------------------ return DOFs ------------------------------------------ +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::dofs() const +inline xt::xtensor VectorPartitioned::dofs() const { return m_dofs; } -// -------------------------------------- return unknown DOFs -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::iiu() const +inline xt::xtensor VectorPartitioned::iiu() const { return m_iiu; } -// ------------------------------------ return prescribed DOFs ------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::iip() const +inline xt::xtensor VectorPartitioned::iip() const { return m_iip; } // ------------------- assemble: dofval[iiu] = dofval_u; dofval[iip] = dofval_p -------------------- -inline void Vector::asDofs(const xt::xtensor &dofval_u, +inline void VectorPartitioned::asDofs(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &dofval) const { assert( dofval_u.size() == m_nnu ); assert( dofval_p.size() == m_nnp ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t i = 0 ; i < m_nnu ; ++i ) dofval(m_iiu(i)) = dofval_u(i); #pragma omp parallel for for ( size_t i = 0 ; i < m_nnp ; ++i ) dofval(m_iip(i)) = dofval_p(i); } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::asDofs(const xt::xtensor &dofval_u, +inline xt::xtensor VectorPartitioned::asDofs(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const { xt::xtensor dofval = xt::empty({m_ndof}); this->asDofs(dofval_u, dofval_p, dofval); return dofval; } -// --------------------------------------- nodevec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asDofs(const xt::xtensor &nodevec, +inline void VectorPartitioned::asDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(n,i)) = nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->asDofs(nodevec, dofval); - - return dofval; -} - -// -------------------------------------- nodevec -> dofval_u -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asDofs_u(const xt::xtensor &nodevec, +inline void VectorPartitioned::asDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval_u) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval_u.size() == m_nnu ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) < m_nnu ) dofval_u(m_part(n,i)) = nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_u(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval_u = xt::empty({m_nnu}); - - this->asDofs_u(nodevec, dofval_u); - - return dofval_u; -} - -// -------------------------------------- nodevec -> dofval_p -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asDofs_p(const xt::xtensor &nodevec, +inline void VectorPartitioned::asDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval_p) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval_p.size() == m_nnp ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) >= m_nnu ) dofval_p(m_part(n,i)-m_nnu) = nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_p(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval_p = xt::empty({m_nnp}); - - this->asDofs_p(nodevec, dofval_p); - - return dofval_p; -} - -// --------------------------------------- elemvec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asDofs(const xt::xtensor &elemvec, +inline void VectorPartitioned::asDofs(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(m_conn(e,m),i)) = elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->asDofs(elemvec, dofval); - - return dofval; -} - -// -------------------------------------- elemvec -> dofval_u -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asDofs_u(const xt::xtensor &elemvec, +inline void VectorPartitioned::asDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval_u) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval_u.size() == m_nnu ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) < m_nnu ) dofval_u(m_part(m_conn(e,m),i)) = elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_u(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval_u = xt::empty({m_nnu}); - - this->asDofs_u(elemvec, dofval_u); - - return dofval_u; -} - -// -------------------------------------- elemvec -> dofval_p -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asDofs_p(const xt::xtensor &elemvec, +inline void VectorPartitioned::asDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval_p) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval_p.size() == m_nnp ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) >= m_nnu ) dofval_p(m_part(m_conn(e,m),i)-m_nnu) = elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asDofs_p(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval_p = xt::empty({m_nnp}); - - this->asDofs_p(elemvec, dofval_p); - - return dofval_p; -} - -// --------------------------------------- dofval -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asNode(const xt::xtensor &dofval, +inline void VectorPartitioned::asNode(const xt::xtensor &dofval, xt::xtensor &nodevec) const { assert( dofval.size() == m_ndof ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(n,i) = dofval(m_dofs(n,i)); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asNode(const xt::xtensor &dofval) const -{ - xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - - this->asNode(dofval, nodevec); - - return nodevec; -} - -// --------------------------------------- dofval -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asNode(const xt::xtensor &dofval_u, +inline void VectorPartitioned::asNode(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &nodevec) const { assert( dofval_u.size() == m_nnu ); assert( dofval_p.size() == m_nnp ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) { for ( size_t i = 0 ; i < m_ndim ; ++i ) { if ( m_part(n,i) < m_nnu ) nodevec(n,i) = dofval_u(m_part(n,i) ); else nodevec(n,i) = dofval_p(m_part(n,i)-m_nnu); } } } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asNode(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const -{ - xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - - this->asNode(dofval_u, dofval_p, nodevec); - - return nodevec; -} - -// -------------------------------------- elemvec -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asNode(const xt::xtensor &elemvec, +inline void VectorPartitioned::asNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(m_conn(e,m),i) = elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asNode(const xt::xtensor &elemvec) const -{ - xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); - - this->asNode(elemvec, nodevec); - - return nodevec; -} - -// --------------------------------------- dofval -> elemvec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asElement(const xt::xtensor &dofval, +inline void VectorPartitioned::asElement(const xt::xtensor &dofval, xt::xtensor &elemvec) const { assert( dofval.size() == m_ndof ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) elemvec(e,m,i) = dofval(m_dofs(m_conn(e,m),i)); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asElement(const xt::xtensor &dofval) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - - this->asElement(dofval, elemvec); - - return elemvec; -} - -// --------------------------------------- dofval -> elemvec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asElement(const xt::xtensor &dofval_u, +inline void VectorPartitioned::asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &elemvec) const { assert( dofval_u.size() == m_nnu ); assert( dofval_p.size() == m_nnp ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { for ( size_t m = 0 ; m < m_nne ; ++m ) { for ( size_t i = 0 ; i < m_ndim ; ++i ) { if ( m_part(m_conn(e,m),i) Vector::asElement(const xt::xtensor &dofval_u, - const xt::xtensor &dofval_p) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - - this->asElement(dofval_u, dofval_p, elemvec); - - return elemvec; -} - -// -------------------------------------- nodevec -> elemvec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::asElement(const xt::xtensor &nodevec, +inline void VectorPartitioned::asElement(const xt::xtensor &nodevec, xt::xtensor &elemvec) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) elemvec(e,m,i) = nodevec(m_conn(e,m),i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::asElement(const xt::xtensor &nodevec) const -{ - xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); - - this->asElement(nodevec, elemvec); - - return elemvec; -} - -// --------------------------------------- nodevec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs(const xt::xtensor &nodevec, +inline void VectorPartitioned::assembleDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_ndof ); dofval.fill(0.0); for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(n,i)) += nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->assembleDofs(nodevec, dofval); - - return dofval; -} - -// -------------------------------------- nodevec -> dofval_u -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_u(const xt::xtensor &nodevec, +inline void VectorPartitioned::assembleDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval_u) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval_u.size() == m_nnu ); dofval_u.fill(0.0); for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) < m_nnu ) dofval_u(m_part(n,i)) += nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs_u(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval_u = xt::empty({m_nnu}); - - this->assembleDofs_u(nodevec, dofval_u); - - return dofval_u; -} - -// -------------------------------------- nodevec -> dofval_p -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_p(const xt::xtensor &nodevec, +inline void VectorPartitioned::assembleDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval_p) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval_p.size() == m_nnp ); dofval_p.fill(0.0); for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) >= m_nnu ) dofval_p(m_part(n,i)-m_nnu) += nodevec(n,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs_p(const xt::xtensor &nodevec) const -{ - xt::xtensor dofval_p = xt::empty({m_nnp}); - - this->assembleDofs_p(nodevec, dofval_p); - - return dofval_p; -} - -// --------------------------------------- elemvec -> dofval --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs(const xt::xtensor &elemvec, +inline void VectorPartitioned::assembleDofs(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_ndof ); dofval.fill(0.0); for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval = xt::empty({m_ndof}); - - this->assembleDofs(elemvec, dofval); - - return dofval; -} - -// -------------------------------------- elemvec -> dofval_u -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_u(const xt::xtensor &elemvec, +inline void VectorPartitioned::assembleDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval_u) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval_u.size() == m_nnu ); dofval_u.fill(0.0); for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) < m_nnu ) dofval_u(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs_u(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval_u = xt::empty({m_nnu}); - - this->assembleDofs_u(elemvec, dofval_u); - - return dofval_u; -} - -// -------------------------------------- elemvec -> dofval_p -------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleDofs_p(const xt::xtensor &elemvec, +inline void VectorPartitioned::assembleDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval_p) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval_p.size() == m_nnp ); dofval_p.fill(0.0); for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) >= m_nnu ) dofval_p(m_dofs(m_conn(e,m),i)-m_nnu) += elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- - -inline xt::xtensor Vector::assembleDofs_p(const xt::xtensor &elemvec) const -{ - xt::xtensor dofval_p = xt::empty({m_nnp}); - - this->assembleDofs_p(elemvec, dofval_p); - - return dofval_p; -} - -// -------------------------------------- elemvec -> nodevec --------------------------------------- +// ------------------------------------------------------------------------------------------------- -inline void Vector::assembleNode(const xt::xtensor &elemvec, +inline void VectorPartitioned::assembleNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); nodevec.fill(0.0); for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(m_conn(e,m),i) += elemvec(e,m,i); } -// ------------------------------------------- interface ------------------------------------------- +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asDofs(const xt::xtensor &nodevec) const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + + this->asDofs(nodevec, dofval); + + return dofval; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asDofs_u(const xt::xtensor &nodevec) const +{ + xt::xtensor dofval_u = xt::empty({m_nnu}); + + this->asDofs_u(nodevec, dofval_u); + + return dofval_u; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asDofs_p(const xt::xtensor &nodevec) const +{ + xt::xtensor dofval_p = xt::empty({m_nnp}); + + this->asDofs_p(nodevec, dofval_p); + + return dofval_p; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asDofs(const xt::xtensor &elemvec) const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + + this->asDofs(elemvec, dofval); + + return dofval; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asDofs_u(const xt::xtensor &elemvec) const +{ + xt::xtensor dofval_u = xt::empty({m_nnu}); + + this->asDofs_u(elemvec, dofval_u); + + return dofval_u; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asDofs_p(const xt::xtensor &elemvec) const +{ + xt::xtensor dofval_p = xt::empty({m_nnp}); + + this->asDofs_p(elemvec, dofval_p); + + return dofval_p; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asNode(const xt::xtensor &dofval) const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + + this->asNode(dofval, nodevec); + + return nodevec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asNode(const xt::xtensor &dofval_u, + const xt::xtensor &dofval_p) const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + + this->asNode(dofval_u, dofval_p, nodevec); + + return nodevec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asNode(const xt::xtensor &elemvec) const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + + this->asNode(elemvec, nodevec); + + return nodevec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asElement(const xt::xtensor &dofval) const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + + this->asElement(dofval, elemvec); + + return elemvec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asElement(const xt::xtensor &dofval_u, + const xt::xtensor &dofval_p) const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + + this->asElement(dofval_u, dofval_p, elemvec); + + return elemvec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::asElement(const xt::xtensor &nodevec) const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + + this->asElement(nodevec, elemvec); + + return elemvec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::assembleDofs( + const xt::xtensor &nodevec) const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + + this->assembleDofs(nodevec, dofval); + + return dofval; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::assembleDofs_u( + const xt::xtensor &nodevec) const +{ + xt::xtensor dofval_u = xt::empty({m_nnu}); + + this->assembleDofs_u(nodevec, dofval_u); + + return dofval_u; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::assembleDofs_p( + const xt::xtensor &nodevec) const +{ + xt::xtensor dofval_p = xt::empty({m_nnp}); + + this->assembleDofs_p(nodevec, dofval_p); + + return dofval_p; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::assembleDofs( + const xt::xtensor &elemvec) const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + + this->assembleDofs(elemvec, dofval); + + return dofval; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::assembleDofs_u( + const xt::xtensor &elemvec) const +{ + xt::xtensor dofval_u = xt::empty({m_nnu}); + + this->assembleDofs_u(elemvec, dofval_u); + + return dofval_u; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor VectorPartitioned::assembleDofs_p( + const xt::xtensor &elemvec) const +{ + xt::xtensor dofval_p = xt::empty({m_nnp}); + + this->assembleDofs_p(elemvec, dofval_p); + + return dofval_p; +} + +// ------------------------------------------------------------------------------------------------- -inline xt::xtensor Vector::assembleNode(const xt::xtensor &elemvec) const +inline xt::xtensor VectorPartitioned::assembleNode( + const xt::xtensor &elemvec) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->assembleNode(elemvec, nodevec); return nodevec; } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif