diff --git a/include/GooseFEM/ElementHex8.h b/include/GooseFEM/ElementHex8.h index 2256da2..bcf9213 100644 --- a/include/GooseFEM/ElementHex8.h +++ b/include/GooseFEM/ElementHex8.h @@ -1,136 +1,148 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTHEX8_H #define GOOSEFEM_ELEMENTHEX8_H #include "config.h" namespace GooseFEM { namespace Element { namespace Hex8 { template inline double inv(const T& A, T& Ainv); namespace Gauss { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // namespace Gauss 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 } // namespace Nodal class Quadrature { public: // Fixed dimensions: // ndim = 3 - number of dimensions // nne = 8 - number of nodes per element // // Naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - [nelem, nip] // Constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature() = default; Quadrature(const xt::xtensor& x); Quadrature( const xt::xtensor& x, 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 shape function gradients xt::xtensor GradN() const; // Convert "qscalar" to "qtensor" of certain rank template void asTensor(const xt::xtensor& qscalar, xt::xtensor& qtensor) const; // Return integration volume xt::xtensor dV() const; // Dyadic product (and its transpose and symmetric part) // qtensor(i,j) += dNdx(m,i) * elemvec(m,j) void gradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void gradN_vector_T(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void symGradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; // Integral of the scalar product // elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV void int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const; // Integral of the dot product // elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV void int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const; // Integral of the dot product // elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV void int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const; // Auto-allocation of the functions above xt::xtensor GradN_vector(const xt::xtensor& elemvec) const; xt::xtensor GradN_vector_T(const xt::xtensor& elemvec) const; xt::xtensor SymGradN_vector(const xt::xtensor& elemvec) const; xt::xtensor Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const; xt::xtensor Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const; xt::xtensor Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const; // Convert "qscalar" to "qtensor" of certain rank template xt::xtensor AsTensor(const xt::xtensor& qscalar) const; xt::xarray AsTensor(size_t rank, const xt::xtensor& qscalar) const; + template + xt::xtensor AllocateQtensor() const; + + template + xt::xtensor AllocateQtensor(double val) const; + + xt::xarray AllocateQtensor(size_t rank) const; + xt::xarray AllocateQtensor(size_t rank, double val) const; + + xt::xtensor AllocateQscalar() const; + xt::xtensor AllocateQscalar(double val) const; + private: // Compute "vol" and "dNdx" based on current "x" void compute_dN(); private: // Dimensions (flexible) size_t m_nelem; // number of elements size_t m_nip; // number of integration points // Dimensions (fixed for this element type) static const size_t m_nne = 8; // number of nodes per element static const size_t m_ndim = 3; // number of dimensions // Data arrays xt::xtensor m_x; // nodal positions stored per element [nelem, nne, ndim] xt::xtensor m_w; // weight of each integration point [nip] xt::xtensor m_xi; // local coordinate of each integration point [nip, ndim] xt::xtensor m_N; // shape functions [nip, nne] xt::xtensor m_dNxi; // shape function grad. wrt local coor. [nip, nne, ndim] xt::xtensor m_dNx; // shape function grad. wrt global coor. [nelem, nip, nne, ndim] xt::xtensor m_vol; // integration point volume [nelem, nip] }; } // namespace Hex8 } // namespace Element } // namespace GooseFEM #include "ElementHex8.hpp" #endif diff --git a/include/GooseFEM/ElementHex8.hpp b/include/GooseFEM/ElementHex8.hpp index 7bbbb49..5698cce 100644 --- a/include/GooseFEM/ElementHex8.hpp +++ b/include/GooseFEM/ElementHex8.hpp @@ -1,590 +1,638 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTHEX8_HPP #define GOOSEFEM_ELEMENTHEX8_HPP #include "ElementHex8.h" namespace GooseFEM { namespace Element { namespace Hex8 { template inline double inv(const T& A, T& Ainv) { double det = (A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) + A(0, 2) * A(1, 0) * A(2, 1)) - (A(0, 2) * A(1, 1) * A(2, 0) + A(0, 1) * A(1, 0) * A(2, 2) + A(0, 0) * A(1, 2) * A(2, 1)); Ainv(0, 0) = (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)) / det; Ainv(0, 1) = (A(0, 2) * A(2, 1) - A(0, 1) * A(2, 2)) / det; Ainv(0, 2) = (A(0, 1) * A(1, 2) - A(0, 2) * A(1, 1)) / det; Ainv(1, 0) = (A(1, 2) * A(2, 0) - A(1, 0) * A(2, 2)) / det; Ainv(1, 1) = (A(0, 0) * A(2, 2) - A(0, 2) * A(2, 0)) / det; Ainv(1, 2) = (A(0, 2) * A(1, 0) - A(0, 0) * A(1, 2)) / det; Ainv(2, 0) = (A(1, 0) * A(2, 1) - A(1, 1) * A(2, 0)) / det; Ainv(2, 1) = (A(0, 1) * A(2, 0) - A(0, 0) * A(2, 1)) / det; Ainv(2, 2) = (A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0)) / det; return det; } namespace Gauss { inline size_t nip() { return 8; } inline xt::xtensor xi() { size_t nip = 8; size_t ndim = 3; xt::xtensor xi = xt::empty({nip, ndim}); xi(0, 0) = -1.0 / std::sqrt(3.0); xi(0, 1) = -1.0 / std::sqrt(3.0); xi(0, 2) = -1.0 / std::sqrt(3.0); xi(1, 0) = +1.0 / std::sqrt(3.0); xi(1, 1) = -1.0 / std::sqrt(3.0); xi(1, 2) = -1.0 / std::sqrt(3.0); xi(2, 0) = +1.0 / std::sqrt(3.0); xi(2, 1) = +1.0 / std::sqrt(3.0); xi(2, 2) = -1.0 / std::sqrt(3.0); xi(3, 0) = -1.0 / std::sqrt(3.0); xi(3, 1) = +1.0 / std::sqrt(3.0); xi(3, 2) = -1.0 / std::sqrt(3.0); xi(4, 0) = -1.0 / std::sqrt(3.0); xi(4, 1) = -1.0 / std::sqrt(3.0); xi(4, 2) = +1.0 / std::sqrt(3.0); xi(5, 0) = +1.0 / std::sqrt(3.0); xi(5, 1) = -1.0 / std::sqrt(3.0); xi(5, 2) = +1.0 / std::sqrt(3.0); xi(6, 0) = +1.0 / std::sqrt(3.0); xi(6, 1) = +1.0 / std::sqrt(3.0); xi(6, 2) = +1.0 / std::sqrt(3.0); xi(7, 0) = -1.0 / std::sqrt(3.0); xi(7, 1) = +1.0 / std::sqrt(3.0); xi(7, 2) = +1.0 / std::sqrt(3.0); return xi; } inline xt::xtensor w() { size_t nip = 8; xt::xtensor w = xt::empty({nip}); w(0) = 1.0; w(1) = 1.0; w(2) = 1.0; w(3) = 1.0; w(4) = 1.0; w(5) = 1.0; w(6) = 1.0; w(7) = 1.0; return w; } } // namespace Gauss namespace Nodal { inline size_t nip() { return 8; } inline xt::xtensor xi() { size_t nip = 8; size_t ndim = 3; xt::xtensor xi = xt::empty({nip, ndim}); xi(0, 0) = -1.0; xi(0, 1) = -1.0; xi(0, 2) = -1.0; xi(1, 0) = +1.0; xi(1, 1) = -1.0; xi(1, 2) = -1.0; xi(2, 0) = +1.0; xi(2, 1) = +1.0; xi(2, 2) = -1.0; xi(3, 0) = -1.0; xi(3, 1) = +1.0; xi(3, 2) = -1.0; xi(4, 0) = -1.0; xi(4, 1) = -1.0; xi(4, 2) = +1.0; xi(5, 0) = +1.0; xi(5, 1) = -1.0; xi(5, 2) = +1.0; xi(6, 0) = +1.0; xi(6, 1) = +1.0; xi(6, 2) = +1.0; xi(7, 0) = -1.0; xi(7, 1) = +1.0; xi(7, 2) = +1.0; return xi; } inline xt::xtensor w() { size_t nip = 8; xt::xtensor w = xt::empty({nip}); w(0) = 1.0; w(1) = 1.0; w(2) = 1.0; w(3) = 1.0; w(4) = 1.0; w(5) = 1.0; w(6) = 1.0; w(7) = 1.0; return w; } } // namespace Nodal inline Quadrature::Quadrature(const xt::xtensor& x) : Quadrature(x, Gauss::xi(), Gauss::w()) { } inline Quadrature::Quadrature( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w) : m_x(x), m_w(w), m_xi(xi) { GOOSEFEM_ASSERT(m_x.shape(1) == m_nne); GOOSEFEM_ASSERT(m_x.shape(2) == m_ndim); m_nelem = m_x.shape(0); m_nip = m_w.size(); GOOSEFEM_ASSERT(m_xi.shape(0) == m_nip); GOOSEFEM_ASSERT(m_xi.shape(1) == m_ndim); GOOSEFEM_ASSERT(m_w.size() == m_nip); 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 q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.125 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_N(q, 1) = 0.125 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_N(q, 2) = 0.125 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_N(q, 3) = 0.125 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_N(q, 4) = 0.125 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)) * (1.0 + m_xi(q, 2)); m_N(q, 5) = 0.125 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)) * (1.0 + m_xi(q, 2)); m_N(q, 6) = 0.125 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)) * (1.0 + m_xi(q, 2)); m_N(q, 7) = 0.125 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)) * (1.0 + m_xi(q, 2)); } // shape function gradients in local coordinates for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.125 * (1.0 - m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 1, 0) = +0.125 * (1.0 - m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 2, 0) = +0.125 * (1.0 + m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 3, 0) = -0.125 * (1.0 + m_xi(q, 1)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 4, 0) = -0.125 * (1.0 - m_xi(q, 1)) * (1.0 + m_xi(q, 2)); m_dNxi(q, 5, 0) = +0.125 * (1.0 - m_xi(q, 1)) * (1.0 + m_xi(q, 2)); m_dNxi(q, 6, 0) = +0.125 * (1.0 + m_xi(q, 1)) * (1.0 + m_xi(q, 2)); m_dNxi(q, 7, 0) = -0.125 * (1.0 + m_xi(q, 1)) * (1.0 + m_xi(q, 2)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.125 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 1, 1) = -0.125 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 2, 1) = +0.125 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 3, 1) = +0.125 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 2)); m_dNxi(q, 4, 1) = -0.125 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 2)); m_dNxi(q, 5, 1) = -0.125 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 2)); m_dNxi(q, 6, 1) = +0.125 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 2)); m_dNxi(q, 7, 1) = +0.125 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 2)); // - dN / dxi_2 m_dNxi(q, 0, 2) = -0.125 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_dNxi(q, 1, 2) = -0.125 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_dNxi(q, 2, 2) = -0.125 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_dNxi(q, 3, 2) = -0.125 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_dNxi(q, 4, 2) = +0.125 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_dNxi(q, 5, 2) = +0.125 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_dNxi(q, 6, 2) = +0.125 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_dNxi(q, 7, 2) = +0.125 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); } // compute the shape function gradients, based on "x" compute_dN(); } inline size_t Quadrature::nelem() const { return m_nelem; } inline size_t Quadrature::nne() const { return m_nne; } inline size_t Quadrature::ndim() const { return m_ndim; } inline size_t Quadrature::nip() const { return m_nip; } inline xt::xtensor Quadrature::GradN() const { return m_dNx; } template inline void Quadrature::asTensor(const xt::xtensor& arg, xt::xtensor& ret) const { GOOSEFEM_ASSERT(xt::has_shape(arg, {m_nelem, m_nne})); GooseFEM::asTensor<2, rank>(arg, ret); } inline xt::xtensor Quadrature::dV() const { return m_vol; } inline void Quadrature::update_x(const xt::xtensor& x) { GOOSEFEM_ASSERT(x.shape() == m_x.shape()); xt::noalias(m_x) = x; compute_dN(); } inline void Quadrature::compute_dN() { #pragma omp parallel { xt::xtensor J = xt::empty({3, 3}); xt::xtensor Jinv = xt::empty({3, 3}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); J.fill(0.0); for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { for (size_t j = 0; j < m_ndim; ++j) { J(i, j) += dNxi(m, i) * x(m, j); } } } double Jdet = inv(J, Jinv); // dNx(m,i) += Jinv(i,j) * dNxi(m,j); for (size_t m = 0; m < m_nne; ++m) { dNx(m, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1) + Jinv(0, 2) * dNxi(m, 2); dNx(m, 1) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1) + Jinv(1, 2) * dNxi(m, 2); dNx(m, 2) = Jinv(2, 0) * dNxi(m, 0) + Jinv(2, 1) * dNxi(m, 1) + Jinv(2, 2) * dNxi(m, 2); } m_vol(e, q) = m_w(q) * Jdet; } } } } inline void Quadrature::gradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { for (size_t j = 0; j < m_ndim; ++j) { gradu(i, j) += dNx(m, i) * u(m, j); } } } } } } inline void Quadrature::gradN_vector_T( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { for (size_t j = 0; j < m_ndim; ++j) { gradu(j, i) += dNx(m, i) * u(m, j); } } } } } } inline void Quadrature::symGradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); for (size_t m = 0; m < m_nne; ++m) { for (size_t i = 0; i < m_ndim; ++i) { for (size_t j = 0; j < m_ndim; ++j) { eps(i, j) += 0.5 * dNx(m, i) * u(m, j); eps(j, i) += 0.5 * dNx(m, i) * u(m, j); } } } } } } inline void Quadrature::int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, {m_nelem, m_nip})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m * ndim + i, n * ndim + i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { M(m * m_ndim + 0, n * m_ndim + 0) += N(m) * rho * N(n) * vol; M(m * m_ndim + 1, n * m_ndim + 1) += N(m) * rho * N(n) * vol; M(m * m_ndim + 2, n * m_ndim + 2) += N(m) * rho * N(n) * vol; } } } } } inline void Quadrature::int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); for (size_t m = 0; m < m_nne; ++m) { f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0) + dNx(m, 2) * sig(2, 0)) * vol; f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1) + dNx(m, 2) * sig(2, 1)) * vol; f(m, 2) += (dNx(m, 0) * sig(0, 2) + dNx(m, 1) * sig(1, 2) + dNx(m, 2) * sig(2, 2)) * vol; } } } } inline void Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim, m_ndim, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto 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); 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; } } } } } } } } } template inline xt::xtensor Quadrature::AsTensor(const xt::xtensor& qscalar) const { return GooseFEM::AsTensor<2, rank>(qscalar, m_ndim); } inline xt::xarray Quadrature::AsTensor(size_t rank, const xt::xtensor& qscalar) const { return GooseFEM::AsTensor(rank, qscalar, m_ndim); } 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; } +template +inline xt::xtensor Quadrature::AllocateQtensor() const +{ + std::array shape; + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_ndim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xtensor ret = xt::empty(shape); + return ret; +} + +template +inline xt::xtensor Quadrature::AllocateQtensor(double val) const +{ + xt::xtensor ret = this->AllocateQtensor(); + ret.fill(val); + return ret; +} + +inline xt::xarray Quadrature::AllocateQtensor(size_t rank) const +{ + std::vector shape(rank + 2); + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_ndim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xarray ret = xt::empty(shape); + return ret; +} + +inline xt::xarray Quadrature::AllocateQtensor(size_t rank, double val) const +{ + xt::xarray ret = this->AllocateQtensor(rank); + ret.fill(val); + return ret; +} + +inline xt::xtensor Quadrature::AllocateQscalar() const +{ + return this->AllocateQtensor<0>(); +} + +inline xt::xtensor Quadrature::AllocateQscalar(double val) const +{ + return this->AllocateQtensor<0>(val); +} + } // namespace Hex8 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/ElementQuad4.h b/include/GooseFEM/ElementQuad4.h index fa03040..63e566c 100644 --- a/include/GooseFEM/ElementQuad4.h +++ b/include/GooseFEM/ElementQuad4.h @@ -1,142 +1,157 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTQUAD4_H #define GOOSEFEM_ELEMENTQUAD4_H #include "config.h" namespace GooseFEM { namespace Element { namespace Quad4 { template inline double inv(const T& A, T& Ainv); namespace Gauss { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // namespace Gauss 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 } // namespace Nodal namespace MidPoint { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // namespace MidPoint class Quadrature { public: // Fixed dimensions: // ndim = 2 - number of dimensions // nne = 4 - number of nodes per element // // Naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - [nelem, nip] // Constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature() = default; Quadrature(const xt::xtensor& x); Quadrature( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w); // Update the nodal positions (shape of "x" should match the earlier definition) void update_x(const xt::xtensor& x); // Return 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 shape function gradients xt::xtensor GradN() const; // Convert "qscalar" to "qtensor" of certain rank template void asTensor(const xt::xtensor& qscalar, xt::xtensor& qtensor) const; // Return integration volume xt::xtensor dV() const; // Dyadic product (and its transpose and symmetric part) // qtensor(i,j) += dNdx(m,i) * elemvec(m,j) void gradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void gradN_vector_T(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void symGradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; // Integral of the scalar product // elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV void int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const; // Integral of the dot product // elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV void int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const; // Integral of the dot product // elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV void int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const; // Auto-allocation of the functions above xt::xtensor GradN_vector(const xt::xtensor& elemvec) const; xt::xtensor GradN_vector_T(const xt::xtensor& elemvec) const; xt::xtensor SymGradN_vector(const xt::xtensor& elemvec) const; xt::xtensor Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const; xt::xtensor Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const; xt::xtensor Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const; // Convert "qscalar" to "qtensor" of certain rank template xt::xtensor AsTensor(const xt::xtensor& qscalar) const; xt::xarray AsTensor(size_t rank, const xt::xtensor& qscalar) const; + // Return allocated integration point tensor of a certain rank, e.g.: + // - rank == 0 -> qscalar + // - rank == 2 -> qtensor + template + xt::xtensor AllocateQtensor() const; + + template + xt::xtensor AllocateQtensor(double val) const; + + xt::xarray AllocateQtensor(size_t rank) const; + xt::xarray AllocateQtensor(size_t rank, double val) const; + + xt::xtensor AllocateQscalar() const; + xt::xtensor AllocateQscalar(double val) const; + private: // Compute "vol" and "dNdx" based on current "x" void compute_dN(); private: // Dimensions (flexible) size_t m_nelem; // number of elements size_t m_nip; // number of integration points // Dimensions (fixed for this element type) static const size_t m_nne = 4; // number of nodes per element static const size_t m_ndim = 2; // number of dimensions // Data arrays xt::xtensor m_x; // nodal positions stored per element [nelem, nne, ndim] xt::xtensor m_w; // weight of each integration point [nip] xt::xtensor m_xi; // local coordinate of each integration point [nip, ndim] xt::xtensor m_N; // shape functions [nip, nne] xt::xtensor m_dNxi; // shape function grad. wrt local coor. [nip, nne, ndim] xt::xtensor m_dNx; // shape function grad. wrt global coor. [nelem, nip, nne, ndim] xt::xtensor m_vol; // integration point volume [nelem, nip] }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #include "ElementQuad4.hpp" #endif diff --git a/include/GooseFEM/ElementQuad4.hpp b/include/GooseFEM/ElementQuad4.hpp index d46edc0..8b92f2a 100644 --- a/include/GooseFEM/ElementQuad4.hpp +++ b/include/GooseFEM/ElementQuad4.hpp @@ -1,530 +1,578 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTQUAD4_HPP #define GOOSEFEM_ELEMENTQUAD4_HPP #include "ElementQuad4.h" namespace GooseFEM { namespace Element { namespace Quad4 { template inline double inv(const T& A, T& Ainv) { double det = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0); Ainv(0, 0) = A(1, 1) / det; Ainv(0, 1) = -1.0 * A(0, 1) / det; Ainv(1, 0) = -1.0 * A(1, 0) / det; Ainv(1, 1) = A(0, 0) / det; return det; } namespace Gauss { inline size_t nip() { return 4; } inline xt::xtensor xi() { size_t nip = 4; size_t ndim = 2; xt::xtensor xi = xt::empty({nip, ndim}); xi(0, 0) = -1.0 / std::sqrt(3.0); xi(0, 1) = -1.0 / std::sqrt(3.0); xi(1, 0) = +1.0 / std::sqrt(3.0); xi(1, 1) = -1.0 / std::sqrt(3.0); xi(2, 0) = +1.0 / std::sqrt(3.0); xi(2, 1) = +1.0 / std::sqrt(3.0); xi(3, 0) = -1.0 / std::sqrt(3.0); xi(3, 1) = +1.0 / std::sqrt(3.0); return xi; } inline xt::xtensor w() { size_t nip = 4; xt::xtensor w = xt::empty({nip}); w(0) = 1.0; w(1) = 1.0; w(2) = 1.0; w(3) = 1.0; return w; } } // namespace Gauss namespace Nodal { inline size_t nip() { return 4; } inline xt::xtensor xi() { size_t nip = 4; size_t ndim = 2; xt::xtensor xi = xt::empty({nip, ndim}); xi(0, 0) = -1.0; xi(0, 1) = -1.0; xi(1, 0) = +1.0; xi(1, 1) = -1.0; xi(2, 0) = +1.0; xi(2, 1) = +1.0; xi(3, 0) = -1.0; xi(3, 1) = +1.0; return xi; } inline xt::xtensor w() { size_t nip = 4; xt::xtensor w = xt::empty({nip}); w(0) = 1.0; w(1) = 1.0; w(2) = 1.0; w(3) = 1.0; return w; } } // namespace Nodal namespace MidPoint { inline size_t nip() { return 1; } inline xt::xtensor xi() { size_t nip = 1; size_t ndim = 2; xt::xtensor xi = xt::empty({nip, ndim}); xi(0, 0) = 0.0; xi(0, 1) = 0.0; return xi; } inline xt::xtensor w() { size_t nip = 1; xt::xtensor w = xt::empty({nip}); w(0) = 1.0; return w; } } // namespace MidPoint inline Quadrature::Quadrature(const xt::xtensor& x) : Quadrature(x, Gauss::xi(), Gauss::w()) { } inline Quadrature::Quadrature( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w) : m_x(x), m_w(w), m_xi(xi) { GOOSEFEM_ASSERT(m_x.shape(1) == m_nne); GOOSEFEM_ASSERT(m_x.shape(2) == m_ndim); m_nelem = m_x.shape(0); m_nip = m_w.size(); GOOSEFEM_ASSERT(m_xi.shape(0) == m_nip); GOOSEFEM_ASSERT(m_xi.shape(1) == m_ndim); GOOSEFEM_ASSERT(m_w.size() == m_nip); 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}); for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 1) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 2) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_N(q, 3) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 1, 0) = +0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 2, 0) = +0.25 * (1.0 + m_xi(q, 1)); m_dNxi(q, 3, 0) = -0.25 * (1.0 + m_xi(q, 1)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.25 * (1.0 - m_xi(q, 0)); m_dNxi(q, 1, 1) = -0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 2, 1) = +0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 3, 1) = +0.25 * (1.0 - m_xi(q, 0)); } compute_dN(); } inline size_t Quadrature::nelem() const { return m_nelem; } inline size_t Quadrature::nne() const { return m_nne; } inline size_t Quadrature::ndim() const { return m_ndim; } inline size_t Quadrature::nip() const { return m_nip; } inline xt::xtensor Quadrature::GradN() const { return m_dNx; } template inline void Quadrature::asTensor(const xt::xtensor& arg, xt::xtensor& ret) const { GOOSEFEM_ASSERT(xt::has_shape(arg, {m_nelem, m_nne})); GooseFEM::asTensor<2, rank>(arg, ret); } inline xt::xtensor Quadrature::dV() const { return m_vol; } inline void Quadrature::update_x(const xt::xtensor& x) { GOOSEFEM_ASSERT(x.shape() == m_x.shape()); xt::noalias(m_x) = x; compute_dN(); } inline void Quadrature::compute_dN() { #pragma omp parallel { xt::xtensor J = xt::empty({2, 2}); xt::xtensor Jinv = xt::empty({2, 2}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); // J(i,j) += dNxi(m,i) * x(m,j); J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) + dNxi(3, 0) * x(3, 0); J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) + dNxi(3, 0) * x(3, 1); J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) + dNxi(3, 1) * x(3, 0); J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) + dNxi(3, 1) * x(3, 1); double Jdet = inv(J, Jinv); // 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); } m_vol(e, q) = m_w(q) * Jdet; } } } } inline void Quadrature::gradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += dNx(m,i) * u(m,j) gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); gradu(0, 1) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1); gradu(1, 0) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0); gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); } } } inline void Quadrature::gradN_vector_T( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += dNx(m,i) * u(m,j) gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); gradu(1, 0) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1); gradu(0, 1) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0); gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); } } } inline void Quadrature::symGradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += dNx(m,i) * u(m,j) // eps(j,i) = 0.5 * (gradu(i,j) + gradu(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) = 0.5 * (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)); eps(1, 0) = eps(0, 1); } } } inline void Quadrature::int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, {m_nelem, m_nip})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { M(m * m_ndim + 0, n * m_ndim + 0) += N(m) * rho * N(n) * vol; M(m * m_ndim + 1, n * m_ndim + 1) += N(m) * rho * N(n) * vol; } } } } } inline void Quadrature::int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); for (size_t m = 0; m < m_nne; ++m) { f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0)) * vol; f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1)) * vol; } } } } inline void Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_ndim, m_ndim, m_ndim, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto 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); 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; } } } } } } } } } template inline xt::xtensor Quadrature::AsTensor(const xt::xtensor& qscalar) const { return GooseFEM::AsTensor<2, rank>(qscalar, m_ndim); } inline xt::xarray Quadrature::AsTensor(size_t rank, const xt::xtensor& qscalar) const { return GooseFEM::AsTensor(rank, qscalar, m_ndim); } 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; } +template +inline xt::xtensor Quadrature::AllocateQtensor() const +{ + std::array shape; + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_ndim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xtensor ret = xt::empty(shape); + return ret; +} + +template +inline xt::xtensor Quadrature::AllocateQtensor(double val) const +{ + xt::xtensor ret = this->AllocateQtensor(); + ret.fill(val); + return ret; +} + +inline xt::xarray Quadrature::AllocateQtensor(size_t rank) const +{ + std::vector shape(rank + 2); + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_ndim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xarray ret = xt::empty(shape); + return ret; +} + +inline xt::xarray Quadrature::AllocateQtensor(size_t rank, double val) const +{ + xt::xarray ret = this->AllocateQtensor(rank); + ret.fill(val); + return ret; +} + +inline xt::xtensor Quadrature::AllocateQscalar() const +{ + return this->AllocateQtensor<0>(); +} + +inline xt::xtensor Quadrature::AllocateQscalar(double val) const +{ + return this->AllocateQtensor<0>(val); +} + } // namespace Quad4 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/ElementQuad4Axisymmetric.h b/include/GooseFEM/ElementQuad4Axisymmetric.h index 7e9e6a9..deff6a9 100644 --- a/include/GooseFEM/ElementQuad4Axisymmetric.h +++ b/include/GooseFEM/ElementQuad4Axisymmetric.h @@ -1,122 +1,134 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTQUAD4AXISYMMETRIC_H #define GOOSEFEM_ELEMENTQUAD4AXISYMMETRIC_H #include "config.h" namespace GooseFEM { namespace Element { namespace Quad4 { class QuadratureAxisymmetric { public: // Fixed dimensions: // ndim = 2 - number of dimensions // nne = 4 - number of nodes per element // tdim = 3 - number of dimensions of tensors // // Naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, tdim, tdim] // "qscalar" - integration point scalar - [nelem, nip] // Constructor: integration point coordinates and weights are optional (default: Gauss) QuadratureAxisymmetric() = default; QuadratureAxisymmetric(const xt::xtensor& x); QuadratureAxisymmetric( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w); // Update the nodal positions (shape of "x" should match the earlier definition) void update_x(const xt::xtensor& x); // 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 xt::xtensor GradN() const; // Convert "qscalar" to "qtensor" of certain rank template void asTensor(const xt::xtensor& qscalar, xt::xtensor& qtensor) const; // Return integration volume xt::xtensor dV() const; // Dyadic product (and its transpose and symmetric part) // qtensor(i,j) += B(m,i,j,k) * elemvec(m,k) void gradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void gradN_vector_T(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void symGradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; // 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 assembled product // fm = ( Bm^T : qtensor ) dV void int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const; // Integral of the assembled product // Kmn = ( Bm^T : qtensor : Bn ) dV void int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const; // Auto-allocation of the functions above xt::xtensor GradN_vector(const xt::xtensor& elemvec) const; xt::xtensor GradN_vector_T(const xt::xtensor& elemvec) const; xt::xtensor SymGradN_vector(const xt::xtensor& elemvec) const; xt::xtensor Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const; xt::xtensor Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const; xt::xtensor Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const; // Convert "qscalar" to "qtensor" of certain rank template xt::xtensor AsTensor(const xt::xtensor& qscalar) const; xt::xarray AsTensor(size_t rank, const xt::xtensor& qscalar) const; + template + xt::xtensor AllocateQtensor() const; + + template + xt::xtensor AllocateQtensor(double val) const; + + xt::xarray AllocateQtensor(size_t rank) const; + xt::xarray AllocateQtensor(size_t rank, double val) const; + + xt::xtensor AllocateQscalar() const; + xt::xtensor AllocateQscalar(double val) const; + private: // Compute "vol" and "B" based on current "x" void compute_dN(); private: // Dimensions (flexible) size_t m_nelem; // number of elements size_t m_nip; // number of integration points // Dimensions (fixed for this element type) static const size_t m_nne = 4; // number of nodes per element static const size_t m_ndim = 2; // number of dimensions static const size_t m_tdim = 3; // number of dimensions of tensors // 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 grad. wrt local coor. [nip, nne, ndim] xt::xtensor m_B; // B-matrix [nelem, nne, tdim, tdim, tdim] xt::xtensor m_vol; // integration point volume [nelem, nip] }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #include "ElementQuad4Axisymmetric.hpp" #endif diff --git a/include/GooseFEM/ElementQuad4Axisymmetric.hpp b/include/GooseFEM/ElementQuad4Axisymmetric.hpp index ed7c8e7..53df33b 100644 --- a/include/GooseFEM/ElementQuad4Axisymmetric.hpp +++ b/include/GooseFEM/ElementQuad4Axisymmetric.hpp @@ -1,473 +1,521 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTQUADAXISYMMETRIC_HPP #define GOOSEFEM_ELEMENTQUADAXISYMMETRIC_HPP #include "ElementQuad4Axisymmetric.h" namespace GooseFEM { namespace Element { namespace Quad4 { inline QuadratureAxisymmetric::QuadratureAxisymmetric(const xt::xtensor& x) : QuadratureAxisymmetric(x, Gauss::xi(), Gauss::w()) { } inline QuadratureAxisymmetric::QuadratureAxisymmetric( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w) : m_x(x), m_w(w), m_xi(xi) { GOOSEFEM_ASSERT(m_x.shape(1) == m_nne); GOOSEFEM_ASSERT(m_x.shape(2) == m_ndim); m_nelem = m_x.shape(0); m_nip = m_w.size(); GOOSEFEM_ASSERT(m_xi.shape(0) == m_nip); GOOSEFEM_ASSERT(m_xi.shape(1) == m_ndim); GOOSEFEM_ASSERT(m_w.size() == m_nip); m_N = xt::empty({m_nip, m_nne}); m_dNxi = xt::empty({m_nip, m_nne, m_ndim}); m_B = xt::empty({m_nelem, m_nip, m_nne, m_tdim, m_tdim, m_tdim}); m_vol = xt::empty({m_nelem, m_nip}); for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 1) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 2) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_N(q, 3) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 1, 0) = +0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 2, 0) = +0.25 * (1.0 + m_xi(q, 1)); m_dNxi(q, 3, 0) = -0.25 * (1.0 + m_xi(q, 1)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.25 * (1.0 - m_xi(q, 0)); m_dNxi(q, 1, 1) = -0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 2, 1) = +0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 3, 1) = +0.25 * (1.0 - m_xi(q, 0)); } compute_dN(); } inline size_t QuadratureAxisymmetric::nelem() const { return m_nelem; } inline size_t QuadratureAxisymmetric::nne() const { return m_nne; } inline size_t QuadratureAxisymmetric::ndim() const { return m_ndim; } inline size_t QuadratureAxisymmetric::nip() const { return m_nip; } template inline void QuadratureAxisymmetric::asTensor(const xt::xtensor& arg, xt::xtensor& ret) const { GOOSEFEM_ASSERT(xt::has_shape(arg, {m_nelem, m_nne})); GooseFEM::asTensor<2, rank>(arg, ret); } inline xt::xtensor QuadratureAxisymmetric::dV() const { return m_vol; } inline void QuadratureAxisymmetric::update_x(const xt::xtensor& x) { GOOSEFEM_ASSERT(x.shape() == m_x.shape()); xt::noalias(m_x) = x; compute_dN(); } inline void QuadratureAxisymmetric::compute_dN() { // most components remain zero, and are not written m_B.fill(0.0); #pragma omp parallel { xt::xtensor J = xt::empty({2, 2}); xt::xtensor Jinv = xt::empty({2, 2}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto N = xt::adapt(&m_N(q, 0), xt::xshape()); // J(i,j) += dNxi(m,i) * x(m,j); J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) + dNxi(3, 0) * x(3, 0); J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) + dNxi(3, 0) * x(3, 1); J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) + dNxi(3, 1) * x(3, 0); J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) + dNxi(3, 1) * x(3, 1); double Jdet = inv(J, Jinv); // radius for computation of volume double rq = N(0) * x(0, 1) + N(1) * x(1, 1) + N(2) * x(2, 1) + N(3) * x(3, 1); // dNx(m,i) += Jinv(i,j) * dNxi(m,j) for (size_t m = 0; m < m_nne; ++m) { // B(m, r, r, r) = dNdx(m,1) B(m, 0, 0, 0) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); // B(m, r, z, z) = dNdx(m,1) B(m, 0, 2, 2) = Jinv(1, 0) * dNxi(m, 0) + Jinv(1, 1) * dNxi(m, 1); // B(m, t, t, r) B(m, 1, 1, 0) = 1.0 / rq * N(m); // B(m, z, r, r) = dNdx(m,0) B(m, 2, 0, 0) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); // B(m, z, z, z) = dNdx(m,0) B(m, 2, 2, 2) = Jinv(0, 0) * dNxi(m, 0) + Jinv(0, 1) * dNxi(m, 1); } m_vol(e, q) = m_w(q) * Jdet * 2.0 * M_PI * rq; } } } } inline void QuadratureAxisymmetric::gradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += B(m,i,j,k) * u(m,perm(k)) // (where perm(0) = 1, perm(2) = 0) gradu(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); gradu(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); gradu(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); gradu(0, 2) = B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0); gradu(2, 0) = B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1); } } } inline void QuadratureAxisymmetric::gradN_vector_T( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += B(m,i,j,k) * u(m,perm(k)) // (where perm(0) = 1, perm(2) = 0) gradu(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); gradu(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); gradu(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); gradu(2, 0) = B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0); gradu(0, 2) = B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1); } } } inline void QuadratureAxisymmetric::symGradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += B(m,i,j,k) * u(m,perm(k)) // eps(j,i) = 0.5 * (gradu(i,j) + gradu(j,i)) // (where perm(0) = 1, perm(2) = 0) eps(0, 0) = B(0, 0, 0, 0) * u(0, 1) + B(1, 0, 0, 0) * u(1, 1) + B(2, 0, 0, 0) * u(2, 1) + B(3, 0, 0, 0) * u(3, 1); eps(1, 1) = B(0, 1, 1, 0) * u(0, 1) + B(1, 1, 1, 0) * u(1, 1) + B(2, 1, 1, 0) * u(2, 1) + B(3, 1, 1, 0) * u(3, 1); eps(2, 2) = B(0, 2, 2, 2) * u(0, 0) + B(1, 2, 2, 2) * u(1, 0) + B(2, 2, 2, 2) * u(2, 0) + B(3, 2, 2, 2) * u(3, 0); eps(2, 0) = 0.5 * (B(0, 0, 2, 2) * u(0, 0) + B(1, 0, 2, 2) * u(1, 0) + B(2, 0, 2, 2) * u(2, 0) + B(3, 0, 2, 2) * u(3, 0) + B(0, 2, 0, 0) * u(0, 1) + B(1, 2, 0, 0) * u(1, 1) + B(2, 2, 0, 0) * u(2, 1) + B(3, 2, 0, 0) * u(3, 1)); eps(0, 2) = eps(2, 0); } } } inline void QuadratureAxisymmetric::int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, {m_nelem, m_nip})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { M(m * m_ndim + 0, n * m_ndim + 0) += N(m) * rho * N(n) * vol; M(m * m_ndim + 1, n * m_ndim + 1) += N(m) * rho * N(n) * vol; } } } } } inline void QuadratureAxisymmetric::int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); // f(m,i) += B(m,i,j,perm(k)) * sig(i,j) * dV // (where perm(0) = 1, perm(2) = 0) for (size_t m = 0; m < m_nne; ++m) { f(m, 0) += vol * (B(m, 2, 2, 2) * sig(2, 2) + B(m, 0, 2, 2) * sig(0, 2)); f(m, 1) += vol * (B(m, 0, 0, 0) * sig(0, 0) + B(m, 1, 1, 0) * sig(1, 1) + B(m, 2, 0, 0) * sig(2, 0)); } } } } inline void QuadratureAxisymmetric::int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim, m_tdim, m_tdim})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto B = xt::adapt(&m_B(e, q, 0, 0, 0, 0), xt::xshape()); auto C = xt::adapt(&qtensor(e, q, 0, 0, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); // K(m*m_ndim+perm(c), n*m_ndim+perm(f)) = B(m,a,b,c) * C(a,b,d,e) * B(n,e,d,f) * vol; // (where perm(0) = 1, perm(2) = 0) for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 0, 0, 0) * C(0, 0, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 0, 0, 0) * C(0, 0, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 0, 0, 0) * C(0, 0, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 1, 1, 0) * C(1, 1, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 1, 1, 0) * C(1, 1, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 1, 1, 0) * C(1, 1, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 2, 2, 2) * C(2, 2, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 2, 2, 2) * C(2, 2, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 2, 2, 2) * C(2, 2, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 0, 2, 2) * C(0, 2, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 0) += B(m, 0, 2, 2) * C(0, 2, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 0, n * m_ndim + 1) += B(m, 0, 2, 2) * C(0, 2, 0, 2) * B(n, 2, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 0, 0) * B(n, 0, 0, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 1, 1) * B(n, 1, 1, 0) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 2, 0, 0) * C(2, 0, 2, 2) * B(n, 2, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 0) += B(m, 2, 0, 0) * C(2, 0, 2, 0) * B(n, 0, 2, 2) * vol; K(m * m_ndim + 1, n * m_ndim + 1) += B(m, 2, 0, 0) * C(2, 0, 0, 2) * B(n, 2, 0, 0) * vol; } } } } } template inline xt::xtensor QuadratureAxisymmetric::AsTensor(const xt::xtensor& qscalar) const { return GooseFEM::AsTensor<2, rank>(qscalar, m_tdim); } inline xt::xarray QuadratureAxisymmetric::AsTensor(size_t rank, const xt::xtensor& qscalar) const { return GooseFEM::AsTensor(rank, qscalar, m_tdim); } inline xt::xtensor QuadratureAxisymmetric::GradN_vector(const xt::xtensor& elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); this->gradN_vector(elemvec, qtensor); return qtensor; } inline xt::xtensor QuadratureAxisymmetric::GradN_vector_T(const xt::xtensor& elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); this->gradN_vector_T(elemvec, qtensor); return qtensor; } inline xt::xtensor QuadratureAxisymmetric::SymGradN_vector(const xt::xtensor& elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); this->symGradN_vector(elemvec, qtensor); return qtensor; } inline xt::xtensor QuadratureAxisymmetric::Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const { xt::xtensor elemmat = xt::empty({m_nelem, m_nne * m_ndim, m_nne * m_ndim}); this->int_N_scalar_NT_dV(qscalar, elemmat); return elemmat; } inline xt::xtensor QuadratureAxisymmetric::Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const { xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); this->int_gradN_dot_tensor2_dV(qtensor, elemvec); return elemvec; } inline xt::xtensor QuadratureAxisymmetric::Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const { xt::xtensor elemmat = xt::empty({m_nelem, m_ndim * m_nne, m_ndim * m_nne}); this->int_gradN_dot_tensor4_dot_gradNT_dV(qtensor, elemmat); return elemmat; } +template +inline xt::xtensor QuadratureAxisymmetric::AllocateQtensor() const +{ + std::array shape; + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_tdim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xtensor ret = xt::empty(shape); + return ret; +} + +template +inline xt::xtensor QuadratureAxisymmetric::AllocateQtensor(double val) const +{ + xt::xtensor ret = this->AllocateQtensor(); + ret.fill(val); + return ret; +} + +inline xt::xarray QuadratureAxisymmetric::AllocateQtensor(size_t rank) const +{ + std::vector shape(rank + 2); + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_tdim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xarray ret = xt::empty(shape); + return ret; +} + +inline xt::xarray QuadratureAxisymmetric::AllocateQtensor(size_t rank, double val) const +{ + xt::xarray ret = this->AllocateQtensor(rank); + ret.fill(val); + return ret; +} + +inline xt::xtensor QuadratureAxisymmetric::AllocateQscalar() const +{ + return this->AllocateQtensor<0>(); +} + +inline xt::xtensor QuadratureAxisymmetric::AllocateQscalar(double val) const +{ + return this->AllocateQtensor<0>(val); +} + } // namespace Quad4 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/ElementQuad4Planar.h b/include/GooseFEM/ElementQuad4Planar.h index 8eb6652..cea9f59 100644 --- a/include/GooseFEM/ElementQuad4Planar.h +++ b/include/GooseFEM/ElementQuad4Planar.h @@ -1,127 +1,139 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTQUAD4PLANAR_H #define GOOSEFEM_ELEMENTQUAD4PLANAR_H #include "config.h" namespace GooseFEM { namespace Element { namespace Quad4 { class QuadraturePlanar { public: // Fixed dimensions: // ndim = 2 - number of dimensions // nne = 4 - number of nodes per element // tdim = 3 - number of dimensions of tensors // // Naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, tdim, tdim] // "qscalar" - integration point scalar - [nelem, nip] // Constructor: integration point coordinates and weights are optional (default: Gauss) QuadraturePlanar() = default; QuadraturePlanar(const xt::xtensor& x, double thick = 1.0); QuadraturePlanar( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w, double thick = 1.0); // 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 shape function gradients xt::xtensor GradN() const; // Convert "qscalar" to "qtensor" of certain rank template void asTensor(const xt::xtensor& qscalar, xt::xtensor& qtensor) const; // Return integration volume xt::xtensor dV() const; // Dyadic product (and its transpose and symmetric part) // qtensor(i,j) += dNdx(m,i) * elemvec(m,j) void gradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void gradN_vector_T(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; void symGradN_vector(const xt::xtensor& elemvec, xt::xtensor& qtensor) const; // Integral of the scalar product // elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV void int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const; // Integral of the dot product // elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV void int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const; // Integral of the dot product // elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV void int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const; // Auto-allocation of the functions above xt::xtensor GradN_vector(const xt::xtensor& elemvec) const; xt::xtensor GradN_vector_T(const xt::xtensor& elemvec) const; xt::xtensor SymGradN_vector(const xt::xtensor& elemvec) const; xt::xtensor Int_N_scalar_NT_dV(const xt::xtensor& qscalar) const; xt::xtensor Int_gradN_dot_tensor2_dV(const xt::xtensor& qtensor) const; xt::xtensor Int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor& qtensor) const; // Convert "qscalar" to "qtensor" of certain rank template xt::xtensor AsTensor(const xt::xtensor& qscalar) const; xt::xarray AsTensor(size_t rank, const xt::xtensor& qscalar) const; + template + xt::xtensor AllocateQtensor() const; + + template + xt::xtensor AllocateQtensor(double val) const; + + xt::xarray AllocateQtensor(size_t rank) const; + xt::xarray AllocateQtensor(size_t rank, double val) const; + + xt::xtensor AllocateQscalar() const; + xt::xtensor AllocateQscalar(double val) const; + private: // Compute "vol" and "dNdx" based on current "x" void compute_dN(); private: // Dimensions (flexible) size_t m_nelem; // number of elements size_t m_nip; // number of integration points // Dimensions (fixed for this element type) static const size_t m_nne = 4; // number of nodes per element static const size_t m_ndim = 2; // number of dimensions static const size_t m_tdim = 3; // number of dimensions of tensors // 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 grad. wrt local coor. [nip, nne, ndim] xt::xtensor m_dNx; // shape function grad. wrt global coor. [nelem, nip, nne, ndim] xt::xtensor m_vol; // integration point volume [nelem, nip] // Thickness double m_thick; }; } // namespace Quad4 } // namespace Element } // namespace GooseFEM #include "ElementQuad4Planar.hpp" #endif diff --git a/include/GooseFEM/ElementQuad4Planar.hpp b/include/GooseFEM/ElementQuad4Planar.hpp index 7be7b62..913a4c7 100644 --- a/include/GooseFEM/ElementQuad4Planar.hpp +++ b/include/GooseFEM/ElementQuad4Planar.hpp @@ -1,405 +1,453 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_ELEMENTQUAD4PLANAR_HPP #define GOOSEFEM_ELEMENTQUAD4PLANAR_HPP #include "ElementQuad4Planar.h" namespace GooseFEM { namespace Element { namespace Quad4 { inline QuadraturePlanar::QuadraturePlanar(const xt::xtensor& x, double thick) : QuadraturePlanar(x, Gauss::xi(), Gauss::w(), thick) { } inline QuadraturePlanar::QuadraturePlanar( const xt::xtensor& x, const xt::xtensor& xi, const xt::xtensor& w, double thick) : m_x(x), m_w(w), m_xi(xi), m_thick(thick) { GOOSEFEM_ASSERT(m_x.shape(1) == m_nne); GOOSEFEM_ASSERT(m_x.shape(2) == m_ndim); m_nelem = m_x.shape(0); m_nip = m_w.size(); GOOSEFEM_ASSERT(m_xi.shape(0) == m_nip); GOOSEFEM_ASSERT(m_xi.shape(1) == m_ndim); GOOSEFEM_ASSERT(m_w.size() == m_nip); 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}); for (size_t q = 0; q < m_nip; ++q) { m_N(q, 0) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 1) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 - m_xi(q, 1)); m_N(q, 2) = 0.25 * (1.0 + m_xi(q, 0)) * (1.0 + m_xi(q, 1)); m_N(q, 3) = 0.25 * (1.0 - m_xi(q, 0)) * (1.0 + m_xi(q, 1)); } for (size_t q = 0; q < m_nip; ++q) { // - dN / dxi_0 m_dNxi(q, 0, 0) = -0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 1, 0) = +0.25 * (1.0 - m_xi(q, 1)); m_dNxi(q, 2, 0) = +0.25 * (1.0 + m_xi(q, 1)); m_dNxi(q, 3, 0) = -0.25 * (1.0 + m_xi(q, 1)); // - dN / dxi_1 m_dNxi(q, 0, 1) = -0.25 * (1.0 - m_xi(q, 0)); m_dNxi(q, 1, 1) = -0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 2, 1) = +0.25 * (1.0 + m_xi(q, 0)); m_dNxi(q, 3, 1) = +0.25 * (1.0 - m_xi(q, 0)); } compute_dN(); } inline size_t QuadraturePlanar::nelem() const { return m_nelem; } inline size_t QuadraturePlanar::nne() const { return m_nne; } inline size_t QuadraturePlanar::ndim() const { return m_ndim; } inline size_t QuadraturePlanar::nip() const { return m_nip; } inline xt::xtensor QuadraturePlanar::GradN() const { return m_dNx; } template inline void QuadraturePlanar::asTensor(const xt::xtensor& arg, xt::xtensor& ret) const { GOOSEFEM_ASSERT(xt::has_shape(arg, {m_nelem, m_nne})); GooseFEM::asTensor<2, rank>(arg, ret); } inline xt::xtensor QuadraturePlanar::dV() const { return m_vol; } inline void QuadraturePlanar::update_x(const xt::xtensor& x) { GOOSEFEM_ASSERT(x.shape() == m_x.shape()); xt::noalias(m_x) = x; compute_dN(); } inline void QuadraturePlanar::compute_dN() { #pragma omp parallel { xt::xtensor J = xt::empty({2, 2}); xt::xtensor Jinv = xt::empty({2, 2}); #pragma omp for for (size_t e = 0; e < m_nelem; ++e) { auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape()); auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); // J(i,j) += dNxi(m,i) * x(m,j); J(0, 0) = dNxi(0, 0) * x(0, 0) + dNxi(1, 0) * x(1, 0) + dNxi(2, 0) * x(2, 0) + dNxi(3, 0) * x(3, 0); J(0, 1) = dNxi(0, 0) * x(0, 1) + dNxi(1, 0) * x(1, 1) + dNxi(2, 0) * x(2, 1) + dNxi(3, 0) * x(3, 1); J(1, 0) = dNxi(0, 1) * x(0, 0) + dNxi(1, 1) * x(1, 0) + dNxi(2, 1) * x(2, 0) + dNxi(3, 1) * x(3, 0); J(1, 1) = dNxi(0, 1) * x(0, 1) + dNxi(1, 1) * x(1, 1) + dNxi(2, 1) * x(2, 1) + dNxi(3, 1) * x(3, 1); double Jdet = inv(J, Jinv); // 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); } m_vol(e, q) = m_w(q) * Jdet * m_thick; } } } } inline void QuadraturePlanar::gradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += dNx(m,i) * u(m,j) gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); gradu(0, 1) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1); gradu(1, 0) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0); gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); } } } inline void QuadraturePlanar::gradN_vector_T( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(j,i) += dNx(m,i) * u(m,j) gradu(0, 0) = dNx(0, 0) * u(0, 0) + dNx(1, 0) * u(1, 0) + dNx(2, 0) * u(2, 0) + dNx(3, 0) * u(3, 0); gradu(1, 0) = dNx(0, 0) * u(0, 1) + dNx(1, 0) * u(1, 1) + dNx(2, 0) * u(2, 1) + dNx(3, 0) * u(3, 1); gradu(0, 1) = dNx(0, 1) * u(0, 0) + dNx(1, 1) * u(1, 0) + dNx(2, 1) * u(2, 0) + dNx(3, 1) * u(3, 0); gradu(1, 1) = dNx(0, 1) * u(0, 1) + dNx(1, 1) * u(1, 1) + dNx(2, 1) * u(2, 1) + dNx(3, 1) * u(3, 1); } } } inline void QuadraturePlanar::symGradN_vector( const xt::xtensor& elemvec, xt::xtensor& qtensor) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); qtensor.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); // gradu(i,j) += dNx(m,i) * u(m,j) // eps(j,i) = 0.5 * (gradu(i,j) + gradu(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) = 0.5 * (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)); eps(1, 0) = eps(0, 1); } } } inline void QuadraturePlanar::int_N_scalar_NT_dV( const xt::xtensor& qscalar, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qscalar, {m_nelem, m_nip})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto N = xt::adapt(&m_N(q, 0), xt::xshape()); auto& vol = m_vol(e, q); auto& rho = qscalar(e, q); // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV for (size_t m = 0; m < m_nne; ++m) { for (size_t n = 0; n < m_nne; ++n) { M(m * m_ndim + 0, n * m_ndim + 0) += N(m) * rho * N(n) * vol; M(m * m_ndim + 1, n * m_ndim + 1) += N(m) * rho * N(n) * vol; } } } } } inline void QuadraturePlanar::int_gradN_dot_tensor2_dV( const xt::xtensor& qtensor, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); elemvec.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape()); auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape()); auto& vol = m_vol(e, q); for (size_t m = 0; m < m_nne; ++m) { f(m, 0) += (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0)) * vol; f(m, 1) += (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1)) * vol; } } } } inline void QuadraturePlanar::int_gradN_dot_tensor4_dot_gradNT_dV( const xt::xtensor& qtensor, xt::xtensor& elemmat) const { GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim, m_tdim, m_tdim})); GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim})); elemmat.fill(0.0); #pragma omp parallel for for (size_t e = 0; e < m_nelem; ++e) { auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape()); for (size_t q = 0; q < m_nip; ++q) { auto 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); 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; } } } } } } } } } template inline xt::xtensor QuadraturePlanar::AsTensor(const xt::xtensor& qscalar) const { return GooseFEM::AsTensor<2, rank>(qscalar, m_tdim); } inline xt::xarray QuadraturePlanar::AsTensor(size_t rank, const xt::xtensor& qscalar) const { return GooseFEM::AsTensor(rank, qscalar, m_tdim); } inline xt::xtensor QuadraturePlanar::GradN_vector(const xt::xtensor& elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); this->gradN_vector(elemvec, qtensor); return qtensor; } inline xt::xtensor QuadraturePlanar::GradN_vector_T(const xt::xtensor& elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); this->gradN_vector_T(elemvec, qtensor); return qtensor; } inline xt::xtensor QuadraturePlanar::SymGradN_vector(const xt::xtensor& elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); this->symGradN_vector(elemvec, qtensor); return qtensor; } inline xt::xtensor QuadraturePlanar::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 QuadraturePlanar::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 QuadraturePlanar::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; } +template +inline xt::xtensor QuadraturePlanar::AllocateQtensor() const +{ + std::array shape; + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_tdim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xtensor ret = xt::empty(shape); + return ret; +} + +template +inline xt::xtensor QuadraturePlanar::AllocateQtensor(double val) const +{ + xt::xtensor ret = this->AllocateQtensor(); + ret.fill(val); + return ret; +} + +inline xt::xarray QuadraturePlanar::AllocateQtensor(size_t rank) const +{ + std::vector shape(rank + 2); + shape[0] = m_nelem; + shape[1] = m_nip; + size_t n = m_tdim; + std::fill(shape.begin() + 2, shape.end(), n); + xt::xarray ret = xt::empty(shape); + return ret; +} + +inline xt::xarray QuadraturePlanar::AllocateQtensor(size_t rank, double val) const +{ + xt::xarray ret = this->AllocateQtensor(rank); + ret.fill(val); + return ret; +} + +inline xt::xtensor QuadraturePlanar::AllocateQscalar() const +{ + return this->AllocateQtensor<0>(); +} + +inline xt::xtensor QuadraturePlanar::AllocateQscalar(double val) const +{ + return this->AllocateQtensor<0>(val); +} + } // namespace Quad4 } // namespace Element } // namespace GooseFEM #endif diff --git a/include/GooseFEM/Vector.h b/include/GooseFEM/Vector.h index ac34a8c..daf2197 100644 --- a/include/GooseFEM/Vector.h +++ b/include/GooseFEM/Vector.h @@ -1,86 +1,94 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTOR_H #define GOOSEFEM_VECTOR_H #include "config.h" namespace GooseFEM { /* "nodevec" - nodal vectors - [nnode, ndim] "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] "dofval" - DOF values - [ndof] */ class Vector { public: // Constructor Vector() = default; Vector(const xt::xtensor& conn, const xt::xtensor& dofs); // Dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs // DOF lists xt::xtensor dofs() const; // DOFs // Copy nodevec to another nodevec void copy(const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const; // Convert to "dofval" (overwrite entries that occur more than once) -- (auto allocation below) void asDofs(const xt::xtensor& nodevec, xt::xtensor& dofval) const; void asDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const; // Convert to "nodevec" (overwrite entries that occur more than once) -- (auto allocation below) void asNode(const xt::xtensor& dofval, xt::xtensor& nodevec) const; void asNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const; // Convert to "elemvec" (overwrite entries that occur more than once) -- (auto allocation below) void asElement(const xt::xtensor& dofval, xt::xtensor& elemvec) const; void asElement(const xt::xtensor& nodevec, xt::xtensor& elemvec) const; // Assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) void assembleDofs(const xt::xtensor& nodevec, xt::xtensor& dofval) const; void assembleDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const; // Assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) void assembleNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const; // Auto-allocation of the functions above xt::xtensor AsDofs(const xt::xtensor& nodevec) const; xt::xtensor AsDofs(const xt::xtensor& elemvec) const; xt::xtensor AsNode(const xt::xtensor& dofval) const; xt::xtensor AsNode(const xt::xtensor& elemvec) const; xt::xtensor AsElement(const xt::xtensor& dofval) const; xt::xtensor AsElement(const xt::xtensor& nodevec) const; xt::xtensor AssembleDofs(const xt::xtensor& nodevec) const; xt::xtensor AssembleDofs(const xt::xtensor& elemvec) const; xt::xtensor AssembleNode(const xt::xtensor& elemvec) const; + // Get zero-allocated dofval, nodevec, elemvec + xt::xtensor AllocateDofval() const; + xt::xtensor AllocateNodevec() const; + xt::xtensor AllocateElemvec() const; + xt::xtensor AllocateDofval(double val) const; + xt::xtensor AllocateNodevec(double val) const; + xt::xtensor AllocateElemvec(double val) const; + private: // Bookkeeping 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 }; } // namespace GooseFEM #include "Vector.hpp" #endif diff --git a/include/GooseFEM/Vector.hpp b/include/GooseFEM/Vector.hpp index a4813df..a5ed263 100644 --- a/include/GooseFEM/Vector.hpp +++ b/include/GooseFEM/Vector.hpp @@ -1,271 +1,310 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTOR_HPP #define GOOSEFEM_VECTOR_HPP #include "Vector.h" namespace GooseFEM { inline Vector::Vector(const xt::xtensor& conn, const xt::xtensor& dofs) : m_conn(conn), m_dofs(dofs) { m_nelem = m_conn.shape(0); m_nne = m_conn.shape(1); m_nnode = m_dofs.shape(0); m_ndim = m_dofs.shape(1); m_ndof = xt::amax(m_dofs)() + 1; GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); } inline size_t Vector::nelem() const { return m_nelem; } inline size_t Vector::nne() const { return m_nne; } inline size_t Vector::nnode() const { return m_nnode; } inline size_t Vector::ndim() const { return m_ndim; } inline size_t Vector::ndof() const { return m_ndof; } inline xt::xtensor Vector::dofs() const { return m_dofs; } inline void Vector::copy(const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim})); xt::noalias(nodevec_dest) = nodevec_src; } inline void Vector::asDofs(const xt::xtensor& nodevec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m, i)) = nodevec(m, i); } } } inline void Vector::asDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); #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); } } } } inline void Vector::asNode(const xt::xtensor& dofval, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { nodevec(m, i) = dofval(m_dofs(m, i)); } } } inline void Vector::asNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); nodevec.fill(0.0); #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); } } } } inline void Vector::asElement(const xt::xtensor& dofval, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, 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)); } } } } inline void Vector::asElement(const xt::xtensor& nodevec, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, 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); } } } } inline void Vector::assembleDofs(const xt::xtensor& nodevec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m, i)) += nodevec(m, i); } } } inline void Vector::assembleDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_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); } } } } inline void Vector::assembleNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); xt::xtensor dofval = this->AssembleDofs(elemvec); this->asNode(dofval, nodevec); } inline xt::xtensor Vector::AsDofs(const xt::xtensor& nodevec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->asDofs(nodevec, dofval); return dofval; } inline xt::xtensor Vector::AsDofs(const xt::xtensor& elemvec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->asDofs(elemvec, dofval); return dofval; } 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; } 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; } 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; } 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; } inline xt::xtensor Vector::AssembleDofs(const xt::xtensor& nodevec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->assembleDofs(nodevec, dofval); return dofval; } inline xt::xtensor Vector::AssembleDofs(const xt::xtensor& elemvec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->assembleDofs(elemvec, dofval); return dofval; } 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; } +inline xt::xtensor Vector::AllocateDofval() const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + return dofval; +} + +inline xt::xtensor Vector::AllocateNodevec() const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + return nodevec; +} + +inline xt::xtensor Vector::AllocateElemvec() const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + return elemvec; +} + +inline xt::xtensor Vector::AllocateDofval(double val) const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + dofval.fill(val); + return dofval; +} + +inline xt::xtensor Vector::AllocateNodevec(double val) const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + nodevec.fill(val); + return nodevec; +} + +inline xt::xtensor Vector::AllocateElemvec(double val) const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + elemvec.fill(val); + return elemvec; +} + } // namespace GooseFEM #endif diff --git a/include/GooseFEM/VectorPartitioned.h b/include/GooseFEM/VectorPartitioned.h index 3b6ee05..fc277fd 100644 --- a/include/GooseFEM/VectorPartitioned.h +++ b/include/GooseFEM/VectorPartitioned.h @@ -1,169 +1,177 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTORPARTITIONED_H #define GOOSEFEM_VECTORPARTITIONED_H #include "config.h" namespace GooseFEM { /* "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 VectorPartitioned { public: // Constructor VectorPartitioned() = default; 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 // Copy (part of) nodevec/dofval to another nodevec/dofval void copy( const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const; void copy_u( const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const; // "iiu" updated void copy_p( const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const; // "iip" updated // Convert to "dofval" (overwrite entries that occur more than once) 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& dofval, xt::xtensor& dofval_u) 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& dofval, xt::xtensor& dofval_p) 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; // 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& nodevec) const; xt::xtensor AsDofs(const xt::xtensor& elemvec) const; xt::xtensor AsDofs_u(const xt::xtensor& dofval) 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& dofval) const; xt::xtensor AsDofs_p(const xt::xtensor& nodevec) const; xt::xtensor AsDofs_p(const xt::xtensor& elemvec) 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; 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; 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; xt::xtensor AssembleNode(const xt::xtensor& elemvec) const; xt::xtensor Copy( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const; xt::xtensor Copy_u( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const; xt::xtensor Copy_p( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const; + // Get zero-allocated dofval, nodevec, elemvec + xt::xtensor AllocateDofval() const; + xt::xtensor AllocateNodevec() const; + xt::xtensor AllocateElemvec() const; + xt::xtensor AllocateDofval(double val) const; + xt::xtensor AllocateNodevec(double val) const; + xt::xtensor AllocateElemvec(double val) 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 GooseFEM #include "VectorPartitioned.hpp" #endif diff --git a/include/GooseFEM/VectorPartitioned.hpp b/include/GooseFEM/VectorPartitioned.hpp index 9d275be..6aaeda2 100644 --- a/include/GooseFEM/VectorPartitioned.hpp +++ b/include/GooseFEM/VectorPartitioned.hpp @@ -1,705 +1,744 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTORPARTITIONED_HPP #define GOOSEFEM_VECTORPARTITIONED_HPP #include "Mesh.h" #include "VectorPartitioned.h" namespace GooseFEM { inline VectorPartitioned::VectorPartitioned( const xt::xtensor& conn, const xt::xtensor& dofs, const xt::xtensor& iip) : m_conn(conn), m_dofs(dofs), m_iip(iip) { m_nelem = m_conn.shape(0); m_nne = m_conn.shape(1); m_nnode = m_dofs.shape(0); m_ndim = m_dofs.shape(1); m_iiu = xt::setdiff1d(m_dofs, m_iip); m_ndof = xt::amax(m_dofs)() + 1; m_nnp = m_iip.size(); m_nnu = m_iiu.size(); m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs); GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode); GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)()); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); } inline size_t VectorPartitioned::nelem() const { return m_nelem; } inline size_t VectorPartitioned::nne() const { return m_nne; } inline size_t VectorPartitioned::nnode() const { return m_nnode; } inline size_t VectorPartitioned::ndim() const { return m_ndim; } inline size_t VectorPartitioned::ndof() const { return m_ndof; } inline size_t VectorPartitioned::nnu() const { return m_nnu; } inline size_t VectorPartitioned::nnp() const { return m_nnp; } inline xt::xtensor VectorPartitioned::dofs() const { return m_dofs; } inline xt::xtensor VectorPartitioned::iiu() const { return m_iiu; } inline xt::xtensor VectorPartitioned::iip() const { return m_iip; } inline void VectorPartitioned::copy( const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim})); xt::noalias(nodevec_dest) = nodevec_src; } inline void VectorPartitioned::copy_u( const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { nodevec_dest(m, i) = nodevec_src(m, i); } } } } inline void VectorPartitioned::copy_p( const xt::xtensor& nodevec_src, xt::xtensor& nodevec_dest) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) >= m_nnu) { nodevec_dest(m, i) = nodevec_src(m, i); } } } } inline void VectorPartitioned::asDofs( const xt::xtensor& dofval_u, const xt::xtensor& dofval_p, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); #pragma omp parallel for for (size_t d = 0; d < m_nnu; ++d) { dofval(m_iiu(d)) = dofval_u(d); } #pragma omp parallel for for (size_t d = 0; d < m_nnp; ++d) { dofval(m_iip(d)) = dofval_p(d); } } inline void VectorPartitioned::asDofs( const xt::xtensor& nodevec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m, i)) = nodevec(m, i); } } } inline void VectorPartitioned::asDofs_u( const xt::xtensor& dofval, xt::xtensor& dofval_u) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); #pragma omp parallel for for (size_t d = 0; d < m_nnu; ++d) { dofval_u(d) = dofval(m_iiu(d)); } } inline void VectorPartitioned::asDofs_u( const xt::xtensor& nodevec, xt::xtensor& dofval_u) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); dofval_u.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { dofval_u(m_part(m, i)) = nodevec(m, i); } } } } inline void VectorPartitioned::asDofs_p( const xt::xtensor& dofval, xt::xtensor& dofval_p) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); #pragma omp parallel for for (size_t d = 0; d < m_nnp; ++d) { dofval_p(d) = dofval(m_iip(d)); } } inline void VectorPartitioned::asDofs_p( const xt::xtensor& nodevec, xt::xtensor& dofval_p) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); dofval_p.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) >= m_nnu) { dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i); } } } } inline void VectorPartitioned::asDofs( const xt::xtensor& elemvec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); #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); } } } } inline void VectorPartitioned::asDofs_u( const xt::xtensor& elemvec, xt::xtensor& dofval_u) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); dofval_u.fill(0.0); #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); } } } } } inline void VectorPartitioned::asDofs_p( const xt::xtensor& elemvec, xt::xtensor& dofval_p) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); dofval_p.fill(0.0); #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); } } } } } inline void VectorPartitioned::asNode( const xt::xtensor& dofval, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { nodevec(m, i) = dofval(m_dofs(m, i)); } } } inline void VectorPartitioned::asNode( const xt::xtensor& dofval_u, const xt::xtensor& dofval_p, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { nodevec(m, i) = dofval_u(m_part(m, i)); } else { nodevec(m, i) = dofval_p(m_part(m, i) - m_nnu); } } } } inline void VectorPartitioned::asNode( const xt::xtensor& elemvec, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); nodevec.fill(0.0); #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); } } } } inline void VectorPartitioned::asElement( const xt::xtensor& dofval, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, 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)); } } } } inline void VectorPartitioned::asElement( const xt::xtensor& dofval_u, const xt::xtensor& dofval_p, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, 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) < m_nnu) { elemvec(e, m, i) = dofval_u(m_part(m_conn(e, m), i)); } else { elemvec(e, m, i) = dofval_p(m_part(m_conn(e, m), i) - m_nnu); } } } } } inline void VectorPartitioned::asElement( const xt::xtensor& nodevec, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, 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); } } } } inline void VectorPartitioned::assembleDofs( const xt::xtensor& nodevec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval.size() == m_ndof); dofval.fill(0.0); for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { dofval(m_dofs(m, i)) += nodevec(m, i); } } } inline void VectorPartitioned::assembleDofs_u( const xt::xtensor& nodevec, xt::xtensor& dofval_u) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_u.size() == m_nnu); dofval_u.fill(0.0); for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) < m_nnu) { dofval_u(m_part(m, i)) += nodevec(m, i); } } } } inline void VectorPartitioned::assembleDofs_p( const xt::xtensor& nodevec, xt::xtensor& dofval_p) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_p.size() == m_nnp); dofval_p.fill(0.0); for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_part(m, i) >= m_nnu) { dofval_p(m_part(m, i) - m_nnu) += nodevec(m, i); } } } } inline void VectorPartitioned::assembleDofs( const xt::xtensor& elemvec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_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); } } } } inline void VectorPartitioned::assembleDofs_u( const xt::xtensor& elemvec, xt::xtensor& dofval_u) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_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_part(m_conn(e, m), i)) += elemvec(e, m, i); } } } } } inline void VectorPartitioned::assembleDofs_p( const xt::xtensor& elemvec, xt::xtensor& dofval_p) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_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_part(m_conn(e, m), i) - m_nnu) += elemvec(e, m, i); } } } } } inline void VectorPartitioned::assembleNode( const xt::xtensor& elemvec, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); xt::xtensor dofval = this->AssembleDofs(elemvec); this->asNode(dofval, nodevec); } 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; } 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& dofval) const { xt::xtensor dofval_u = xt::empty({m_nnu}); this->asDofs_u(dofval, dofval_u); return dofval_u; } 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& dofval) const { xt::xtensor dofval_p = xt::empty({m_nnp}); this->asDofs_p(dofval, dofval_p); return dofval_p; } 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 VectorPartitioned::AssembleNode(const xt::xtensor& elemvec) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->assembleNode(elemvec, nodevec); return nodevec; } inline xt::xtensor VectorPartitioned::Copy( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const { xt::xtensor ret = nodevec_dest; this->copy(nodevec_src, ret); return ret; } inline xt::xtensor VectorPartitioned::Copy_u( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const { xt::xtensor ret = nodevec_dest; this->copy_u(nodevec_src, ret); return ret; } inline xt::xtensor VectorPartitioned::Copy_p( const xt::xtensor& nodevec_src, const xt::xtensor& nodevec_dest) const { xt::xtensor ret = nodevec_dest; this->copy_p(nodevec_src, ret); return ret; } +inline xt::xtensor VectorPartitioned::AllocateDofval() const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + return dofval; +} + +inline xt::xtensor VectorPartitioned::AllocateNodevec() const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + return nodevec; +} + +inline xt::xtensor VectorPartitioned::AllocateElemvec() const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + return elemvec; +} + +inline xt::xtensor VectorPartitioned::AllocateDofval(double val) const +{ + xt::xtensor dofval = xt::zeros({m_ndof}); + dofval.fill(val); + return dofval; +} + +inline xt::xtensor VectorPartitioned::AllocateNodevec(double val) const +{ + xt::xtensor nodevec = xt::zeros({m_nnode, m_ndim}); + nodevec.fill(val); + return nodevec; +} + +inline xt::xtensor VectorPartitioned::AllocateElemvec(double val) const +{ + xt::xtensor elemvec = xt::zeros({m_nelem, m_nne, m_ndim}); + elemvec.fill(val); + return elemvec; +} + } // namespace GooseFEM #endif diff --git a/include/GooseFEM/VectorPartitionedTyings.h b/include/GooseFEM/VectorPartitionedTyings.h index 72c64e8..80ba8e5 100644 --- a/include/GooseFEM/VectorPartitionedTyings.h +++ b/include/GooseFEM/VectorPartitionedTyings.h @@ -1,120 +1,128 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTORPARTITIONEDTYINGS_H #define GOOSEFEM_VECTORPARTITIONEDTYINGS_H #include "config.h" #include #include namespace GooseFEM { /* "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 VectorPartitionedTyings { public: // Constructor VectorPartitionedTyings() = default; VectorPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp, const Eigen::SparseMatrix& Cdi); // 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 independent, unknown DOFs size_t nnp() const; // number of independent, prescribed DOFs size_t nni() const; // number of independent DOFs size_t nnd() const; // number of dependent DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // independent, unknown DOFs xt::xtensor iip() const; // independent, prescribed DOFs xt::xtensor iii() const; // independent DOFs xt::xtensor iid() const; // dependent DOFs // Copy (part of) nodevec/dofval to another nodevec/dofval void copy_p( const xt::xtensor& dofval_src, xt::xtensor& dofval_dest) const; // "iip" updated // Convert to "dofval" (overwrite entries that occur more than once) void asDofs_i( const xt::xtensor& nodevec, xt::xtensor& dofval_i, bool apply_tyings = true) const; // Convert to "nodevec" (overwrite entries that occur more than once) -- (auto allocation below) void asNode(const xt::xtensor& dofval, xt::xtensor& nodevec) const; // Convert to "elemvec" (overwrite entries that occur more than once) -- (auto allocation below) void asElement(const xt::xtensor& nodevec, xt::xtensor& elemvec) const; // Assemble "dofval" (adds entries that occur more that once) -- (auto allocation below) void assembleDofs(const xt::xtensor& elemvec, xt::xtensor& dofval) const; // Assemble "nodevec" (adds entries that occur more that once) -- (auto allocation below) void assembleNode(const xt::xtensor& elemvec, xt::xtensor& nodevec) const; // Auto-allocation of the functions above xt::xtensor AsDofs_i(const xt::xtensor& nodevec) const; xt::xtensor AsNode(const xt::xtensor& dofval) const; xt::xtensor AsElement(const xt::xtensor& nodevec) const; xt::xtensor AssembleDofs(const xt::xtensor& elemvec) const; xt::xtensor AssembleNode(const xt::xtensor& elemvec) const; + // Get zero-allocated dofval, nodevec, elemvec + xt::xtensor AllocateDofval() const; + xt::xtensor AllocateNodevec() const; + xt::xtensor AllocateElemvec() const; + xt::xtensor AllocateDofval(double val) const; + xt::xtensor AllocateNodevec(double val) const; + xt::xtensor AllocateElemvec(double val) const; + private: // Bookkeeping xt::xtensor m_conn; // connectivity [nelem, nne] xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] xt::xtensor m_iiu; // unknown DOFs [nnu] xt::xtensor m_iip; // prescribed DOFs [nnp] xt::xtensor m_iid; // dependent DOFs [nnd] // 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 independent, unknown DOFs size_t m_nnp; // number of independent, prescribed DOFs size_t m_nni; // number of independent DOFs size_t m_nnd; // number of dependent DOFs // Tyings Eigen::SparseMatrix m_Cdu; Eigen::SparseMatrix m_Cdp; Eigen::SparseMatrix m_Cdi; Eigen::SparseMatrix m_Cud; Eigen::SparseMatrix m_Cpd; Eigen::SparseMatrix m_Cid; // equivalent Eigen functions Eigen::VectorXd Eigen_asDofs_d(const xt::xtensor& nodevec) const; }; } // namespace GooseFEM #include "VectorPartitionedTyings.hpp" #endif diff --git a/include/GooseFEM/VectorPartitionedTyings.hpp b/include/GooseFEM/VectorPartitionedTyings.hpp index dde0412..2df978e 100644 --- a/include/GooseFEM/VectorPartitionedTyings.hpp +++ b/include/GooseFEM/VectorPartitionedTyings.hpp @@ -1,278 +1,317 @@ /* (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM */ #ifndef GOOSEFEM_VECTORPARTITIONEDTYINGS_HPP #define GOOSEFEM_VECTORPARTITIONEDTYINGS_HPP #include "VectorPartitionedTyings.h" namespace GooseFEM { inline VectorPartitionedTyings::VectorPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp, const Eigen::SparseMatrix& Cdi) : m_conn(conn), m_dofs(dofs), m_Cdu(Cdu), m_Cdp(Cdp), m_Cdi(Cdi) { GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows()); GOOSEFEM_ASSERT(Cdi.rows() == Cdp.rows()); m_nnu = static_cast(m_Cdu.cols()); m_nnp = static_cast(m_Cdp.cols()); m_nnd = static_cast(m_Cdp.rows()); m_nni = m_nnu + m_nnp; m_ndof = m_nni + m_nnd; m_iiu = xt::arange(m_nnu); m_iip = xt::arange(m_nnu, m_nnu + m_nnp); m_iid = xt::arange(m_nni, m_nni + m_nnd); m_nelem = m_conn.shape(0); m_nne = m_conn.shape(1); m_nnode = m_dofs.shape(0); m_ndim = m_dofs.shape(1); m_Cud = m_Cdu.transpose(); m_Cpd = m_Cdp.transpose(); m_Cid = m_Cdi.transpose(); GOOSEFEM_ASSERT(static_cast(m_Cdi.cols()) == m_nni); GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim); GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)() + 1); } inline size_t VectorPartitionedTyings::nelem() const { return m_nelem; } inline size_t VectorPartitionedTyings::nne() const { return m_nne; } inline size_t VectorPartitionedTyings::nnode() const { return m_nnode; } inline size_t VectorPartitionedTyings::ndim() const { return m_ndim; } inline size_t VectorPartitionedTyings::ndof() const { return m_ndof; } inline size_t VectorPartitionedTyings::nnu() const { return m_nnu; } inline size_t VectorPartitionedTyings::nnp() const { return m_nnp; } inline size_t VectorPartitionedTyings::nni() const { return m_nni; } inline size_t VectorPartitionedTyings::nnd() const { return m_nnd; } inline xt::xtensor VectorPartitionedTyings::dofs() const { return m_dofs; } inline xt::xtensor VectorPartitionedTyings::iiu() const { return m_iiu; } inline xt::xtensor VectorPartitionedTyings::iip() const { return m_iip; } inline xt::xtensor VectorPartitionedTyings::iii() const { return xt::arange(m_nni); } inline xt::xtensor VectorPartitionedTyings::iid() const { return m_iid; } inline void VectorPartitionedTyings::copy_p( const xt::xtensor& dofval_src, xt::xtensor& dofval_dest) const { GOOSEFEM_ASSERT(dofval_src.size() == m_ndof || dofval_src.size() == m_nni); GOOSEFEM_ASSERT(dofval_dest.size() == m_ndof || dofval_dest.size() == m_nni); #pragma omp parallel for for (size_t i = m_nnu; i < m_nni; ++i) { dofval_dest(i) = dofval_src(i); } } inline void VectorPartitionedTyings::asDofs_i( const xt::xtensor& nodevec, xt::xtensor& dofval_i, bool apply_tyings) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(dofval_i.size() == m_nni); dofval_i.fill(0.0); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_dofs(m, i) < m_nni) { dofval_i(m_dofs(m, i)) = nodevec(m, i); } } } if (!apply_tyings) { return; } Eigen::VectorXd Dofval_d = this->Eigen_asDofs_d(nodevec); Eigen::VectorXd Dofval_i = m_Cid * Dofval_d; #pragma omp parallel for for (size_t i = 0; i < m_nni; ++i) { dofval_i(i) += Dofval_i(i); } } inline void VectorPartitionedTyings::asNode( const xt::xtensor& dofval, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(dofval.size() == m_ndof); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { nodevec(m, i) = dofval(m_dofs(m, i)); } } } inline void VectorPartitionedTyings::asElement( const xt::xtensor& nodevec, xt::xtensor& elemvec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, 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); } } } } inline void VectorPartitionedTyings::assembleDofs( const xt::xtensor& elemvec, xt::xtensor& dofval) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_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); } } } } inline void VectorPartitionedTyings::assembleNode( const xt::xtensor& elemvec, xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); xt::xtensor dofval = this->AssembleDofs(elemvec); this->asNode(dofval, nodevec); } inline xt::xtensor VectorPartitionedTyings::AsDofs_i(const xt::xtensor& nodevec) const { xt::xtensor dofval = xt::empty({m_nni}); this->asDofs_i(nodevec, dofval); return dofval; } inline xt::xtensor VectorPartitionedTyings::AsNode(const xt::xtensor& dofval) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->asNode(dofval, nodevec); return nodevec; } inline xt::xtensor VectorPartitionedTyings::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 VectorPartitionedTyings::AssembleDofs(const xt::xtensor& elemvec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->assembleDofs(elemvec, dofval); return dofval; } inline xt::xtensor VectorPartitionedTyings::AssembleNode(const xt::xtensor& elemvec) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->assembleNode(elemvec, nodevec); return nodevec; } inline Eigen::VectorXd VectorPartitionedTyings::Eigen_asDofs_d(const xt::xtensor& nodevec) const { GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim})); Eigen::VectorXd dofval_d(m_nnd, 1); #pragma omp parallel for for (size_t m = 0; m < m_nnode; ++m) { for (size_t i = 0; i < m_ndim; ++i) { if (m_dofs(m, i) >= m_nni) { dofval_d(m_dofs(m, i) - m_nni) = nodevec(m, i); } } } return dofval_d; } +inline xt::xtensor VectorPartitionedTyings::AllocateDofval() const +{ + xt::xtensor dofval = xt::empty({m_ndof}); + return dofval; +} + +inline xt::xtensor VectorPartitionedTyings::AllocateNodevec() const +{ + xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); + return nodevec; +} + +inline xt::xtensor VectorPartitionedTyings::AllocateElemvec() const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + return elemvec; +} + +inline xt::xtensor VectorPartitionedTyings::AllocateDofval(double val) const +{ + xt::xtensor dofval = xt::zeros({m_ndof}); + dofval.fill(val); + return dofval; +} + +inline xt::xtensor VectorPartitionedTyings::AllocateNodevec(double val) const +{ + xt::xtensor nodevec = xt::zeros({m_nnode, m_ndim}); + nodevec.fill(val); + return nodevec; +} + +inline xt::xtensor VectorPartitionedTyings::AllocateElemvec(double val) const +{ + xt::xtensor elemvec = xt::zeros({m_nelem, m_nne, m_ndim}); + elemvec.fill(val); + return elemvec; +} + } // namespace GooseFEM #endif diff --git a/python/ElementHex8.hpp b/python/ElementHex8.hpp index 9452414..fa5c43c 100644 --- a/python/ElementHex8.hpp +++ b/python/ElementHex8.hpp @@ -1,117 +1,148 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_ElementHex8(py::module& m) { py::class_(m, "Quadrature") .def(py::init&>(), "Quadrature", py::arg("x")) .def( py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(), "Quadrature", py::arg("x"), py::arg("xi"), py::arg("w")) .def( "update_x", &GooseFEM::Element::Hex8::Quadrature::update_x, "Update the nodal positions") .def("nelem", &GooseFEM::Element::Hex8::Quadrature::nelem, "Number of elements") .def("nne", &GooseFEM::Element::Hex8::Quadrature::nne, "Number of nodes per element") .def("ndim", &GooseFEM::Element::Hex8::Quadrature::ndim, "Number of dimensions") .def("nip", &GooseFEM::Element::Hex8::Quadrature::nip, "Number of integration points") .def("dV", &GooseFEM::Element::Hex8::Quadrature::dV, "Integration point volume (qscalar)") .def( "GradN_vector", py::overload_cast&>( &GooseFEM::Element::Hex8::Quadrature::GradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "GradN_vector_T", py::overload_cast&>( &GooseFEM::Element::Hex8::Quadrature::GradN_vector_T, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "SymGradN_vector", py::overload_cast&>( &GooseFEM::Element::Hex8::Quadrature::SymGradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "Int_N_scalar_NT_dV", py::overload_cast&>( &GooseFEM::Element::Hex8::Quadrature::Int_N_scalar_NT_dV, py::const_), "Integration, returns 'elemmat'", py::arg("qscalar")) .def( "Int_gradN_dot_tensor2_dV", py::overload_cast&>( &GooseFEM::Element::Hex8::Quadrature::Int_gradN_dot_tensor2_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "Int_gradN_dot_tensor4_dot_gradNT_dV", py::overload_cast&>( &GooseFEM::Element::Hex8::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "AsTensor", (xt::xarray(GooseFEM::Element::Hex8::Quadrature::*)( size_t, const xt::xtensor&) const) & GooseFEM::Element::Hex8::Quadrature::AsTensor, "Convert 'qscalar' to 'qtensor' of certain rank") + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Hex8::Quadrature::*)( + size_t) const) & + GooseFEM::Element::Hex8::Quadrature::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank")) + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Hex8::Quadrature::*)( + size_t, double) const) & + GooseFEM::Element::Hex8::Quadrature::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank"), + py::arg("val")) + + .def( + "AllocateQscalar", + py::overload_cast<>( + &GooseFEM::Element::Hex8::Quadrature::AllocateQscalar, py::const_), + "Allocate 'qscalar'") + + .def( + "AllocateQscalar", + py::overload_cast( + &GooseFEM::Element::Hex8::Quadrature::AllocateQscalar, py::const_), + "Allocate 'qscalar'", + py::arg("val")) + .def("__repr__", [](const GooseFEM::Element::Hex8::Quadrature&) { return ""; }); } void init_ElementHex8Gauss(py::module& m) { m.def("nip", &GooseFEM::Element::Hex8::Gauss::nip, "Return number of integration point"); m.def("xi", &GooseFEM::Element::Hex8::Gauss::xi, "Return integration point coordinates"); m.def("w", &GooseFEM::Element::Hex8::Gauss::w, "Return integration point weights"); } void init_ElementHex8Nodal(py::module& m) { m.def("nip", &GooseFEM::Element::Hex8::Nodal::nip, "Return number of integration point"); m.def("xi", &GooseFEM::Element::Hex8::Nodal::xi, "Return integration point coordinates"); m.def("w", &GooseFEM::Element::Hex8::Nodal::w, "Return integration point weights"); } diff --git a/python/ElementQuad4.hpp b/python/ElementQuad4.hpp index c068e76..44d35e8 100644 --- a/python/ElementQuad4.hpp +++ b/python/ElementQuad4.hpp @@ -1,117 +1,148 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_ElementQuad4(py::module& m) { py::class_(m, "Quadrature") .def(py::init&>(), "Quadrature", py::arg("x")) .def( py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(), "Quadrature", py::arg("x"), py::arg("xi"), py::arg("w")) .def( "update_x", &GooseFEM::Element::Quad4::Quadrature::update_x, "Update the nodal positions") .def("nelem", &GooseFEM::Element::Quad4::Quadrature::nelem, "Number of elements") .def("nne", &GooseFEM::Element::Quad4::Quadrature::nne, "Number of nodes per element") .def("ndim", &GooseFEM::Element::Quad4::Quadrature::ndim, "Number of dimensions") .def("nip", &GooseFEM::Element::Quad4::Quadrature::nip, "Number of integration points") .def("dV", &GooseFEM::Element::Quad4::Quadrature::dV, "Integration point volume (qscalar)") .def( "GradN_vector", py::overload_cast&>( &GooseFEM::Element::Quad4::Quadrature::GradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "GradN_vector_T", py::overload_cast&>( &GooseFEM::Element::Quad4::Quadrature::GradN_vector_T, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "SymGradN_vector", py::overload_cast&>( &GooseFEM::Element::Quad4::Quadrature::SymGradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "Int_N_scalar_NT_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::Quadrature::Int_N_scalar_NT_dV, py::const_), "Integration, returns 'elemmat'", py::arg("qscalar")) .def( "Int_gradN_dot_tensor2_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor2_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "Int_gradN_dot_tensor4_dot_gradNT_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "AsTensor", (xt::xarray(GooseFEM::Element::Quad4::Quadrature::*)( size_t, const xt::xtensor&) const) & GooseFEM::Element::Quad4::Quadrature::AsTensor, "Convert 'qscalar' to 'qtensor' of certain rank") + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Quad4::Quadrature::*)( + size_t) const) & + GooseFEM::Element::Quad4::Quadrature::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank")) + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Quad4::Quadrature::*)( + size_t, double) const) & + GooseFEM::Element::Quad4::Quadrature::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank"), + py::arg("val")) + + .def( + "AllocateQscalar", + py::overload_cast<>( + &GooseFEM::Element::Quad4::Quadrature::AllocateQscalar, py::const_), + "Allocate 'qscalar'") + + .def( + "AllocateQscalar", + py::overload_cast( + &GooseFEM::Element::Quad4::Quadrature::AllocateQscalar, py::const_), + "Allocate 'qscalar'", + py::arg("val")) + .def("__repr__", [](const GooseFEM::Element::Quad4::Quadrature&) { return ""; }); } void init_ElementQuad4Gauss(py::module& m) { m.def("nip", &GooseFEM::Element::Quad4::Gauss::nip, "Return number of integration point"); m.def("xi", &GooseFEM::Element::Quad4::Gauss::xi, "Return integration point coordinates"); m.def("w", &GooseFEM::Element::Quad4::Gauss::w, "Return integration point weights"); } void init_ElementQuad4Nodal(py::module& m) { m.def("nip", &GooseFEM::Element::Quad4::Nodal::nip, "Return number of integration point"); m.def("xi", &GooseFEM::Element::Quad4::Nodal::xi, "Return integration point coordinates"); m.def("w", &GooseFEM::Element::Quad4::Nodal::w, "Return integration point weights"); } diff --git a/python/ElementQuad4Axisymmetric.hpp b/python/ElementQuad4Axisymmetric.hpp index 81f18ac..f53137e 100644 --- a/python/ElementQuad4Axisymmetric.hpp +++ b/python/ElementQuad4Axisymmetric.hpp @@ -1,110 +1,141 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_ElementQuad4Axisymmetric(py::module& m) { py::class_(m, "QuadratureAxisymmetric") .def(py::init&>(), "QuadratureAxisymmetric", py::arg("x")) .def( py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(), "QuadratureAxisymmetric", py::arg("x"), py::arg("xi"), py::arg("w")) .def( "update_x", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::update_x, "Update the nodal positions") .def( "nelem", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nelem, "Number of elements") .def( "nne", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nne, "Number of nodes per element") .def( "ndim", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::ndim, "Number of dimensions") .def( "nip", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nip, "Number of integration points") .def( "dV", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::dV, "Integration point volume (qscalar)") .def( "GradN_vector", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadratureAxisymmetric::GradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "GradN_vector_T", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadratureAxisymmetric::GradN_vector_T, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "SymGradN_vector", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadratureAxisymmetric::SymGradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "Int_N_scalar_NT_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_N_scalar_NT_dV, py::const_), "Integration, returns 'elemmat'", py::arg("qscalar")) .def( "Int_gradN_dot_tensor2_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_gradN_dot_tensor2_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "Int_gradN_dot_tensor4_dot_gradNT_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadratureAxisymmetric:: Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "AsTensor", (xt::xarray(GooseFEM::Element::Quad4::QuadratureAxisymmetric::*)( size_t, const xt::xtensor&) const) & GooseFEM::Element::Quad4::QuadratureAxisymmetric::AsTensor, "Convert 'qscalar' to 'qtensor' of certain rank") + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Quad4::QuadratureAxisymmetric::*)( + size_t) const) & + GooseFEM::Element::Quad4::QuadratureAxisymmetric::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank")) + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Quad4::QuadratureAxisymmetric::*)( + size_t, double) const) & + GooseFEM::Element::Quad4::QuadratureAxisymmetric::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank"), + py::arg("val")) + + .def( + "AllocateQscalar", + py::overload_cast<>( + &GooseFEM::Element::Quad4::QuadratureAxisymmetric::AllocateQscalar, py::const_), + "Allocate 'qscalar'") + + .def( + "AllocateQscalar", + py::overload_cast( + &GooseFEM::Element::Quad4::QuadratureAxisymmetric::AllocateQscalar, py::const_), + "Allocate 'qscalar'", + py::arg("val")) + .def("__repr__", [](const GooseFEM::Element::Quad4::QuadratureAxisymmetric&) { return ""; }); } diff --git a/python/ElementQuad4Planar.hpp b/python/ElementQuad4Planar.hpp index e7f557c..5ed9a9e 100644 --- a/python/ElementQuad4Planar.hpp +++ b/python/ElementQuad4Planar.hpp @@ -1,101 +1,132 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_ElementQuad4Planar(py::module& m) { py::class_(m, "QuadraturePlanar") .def(py::init&>(), "QuadraturePlanar", py::arg("x")) .def( py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(), "QuadraturePlanar", py::arg("x"), py::arg("xi"), py::arg("w")) .def( "update_x", &GooseFEM::Element::Quad4::QuadraturePlanar::update_x, "Update the nodal positions") .def("nelem", &GooseFEM::Element::Quad4::QuadraturePlanar::nelem, "Number of elements") .def("nne", &GooseFEM::Element::Quad4::QuadraturePlanar::nne, "Number of nodes per element") .def("ndim", &GooseFEM::Element::Quad4::QuadraturePlanar::ndim, "Number of dimensions") .def( "nip", &GooseFEM::Element::Quad4::QuadraturePlanar::nip, "Number of integration points") .def( "dV", &GooseFEM::Element::Quad4::QuadraturePlanar::dV, "Integration point volume (qscalar)") .def( "GradN_vector", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadraturePlanar::GradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "GradN_vector_T", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadraturePlanar::GradN_vector_T, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "SymGradN_vector", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadraturePlanar::SymGradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec")) .def( "Int_N_scalar_NT_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadraturePlanar::Int_N_scalar_NT_dV, py::const_), "Integration, returns 'elemmat'", py::arg("qscalar")) .def( "Int_gradN_dot_tensor2_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadraturePlanar::Int_gradN_dot_tensor2_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "Int_gradN_dot_tensor4_dot_gradNT_dV", py::overload_cast&>( &GooseFEM::Element::Quad4::QuadraturePlanar::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor")) .def( "AsTensor", (xt::xarray(GooseFEM::Element::Quad4::QuadraturePlanar::*)( size_t, const xt::xtensor&) const) & GooseFEM::Element::Quad4::QuadraturePlanar::AsTensor, "Convert 'qscalar' to 'qtensor' of certain rank") + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Quad4::QuadraturePlanar::*)( + size_t) const) & + GooseFEM::Element::Quad4::QuadraturePlanar::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank")) + + .def( + "AllocateQtensor", + (xt::xarray(GooseFEM::Element::Quad4::QuadraturePlanar::*)( + size_t, double) const) & + GooseFEM::Element::Quad4::QuadraturePlanar::AllocateQtensor, + "Allocate 'qtensor'", + py::arg("rank"), + py::arg("val")) + + .def( + "AllocateQscalar", + py::overload_cast<>( + &GooseFEM::Element::Quad4::QuadraturePlanar::AllocateQscalar, py::const_), + "Allocate 'qscalar'") + + .def( + "AllocateQscalar", + py::overload_cast( + &GooseFEM::Element::Quad4::QuadraturePlanar::AllocateQscalar, py::const_), + "Allocate 'qscalar'", + py::arg("val")) + .def("__repr__", [](const GooseFEM::Element::Quad4::QuadraturePlanar&) { return ""; }); } diff --git a/python/Vector.hpp b/python/Vector.hpp index a2e09fc..5bc28f1 100644 --- a/python/Vector.hpp +++ b/python/Vector.hpp @@ -1,96 +1,120 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_Vector(py::module& m) { py::class_(m, "Vector") .def( py::init&, const xt::xtensor&>(), "Switch between dofval/nodevec/elemvec", py::arg("conn"), py::arg("dofs")) .def("nelem", &GooseFEM::Vector::nelem, "Number of element") .def("nne", &GooseFEM::Vector::nne, "Number of nodes per element") .def("nnode", &GooseFEM::Vector::nnode, "Number of nodes") .def("ndim", &GooseFEM::Vector::ndim, "Number of dimensions") .def("ndof", &GooseFEM::Vector::ndof, "Number of degrees-of-freedom") .def("dofs", &GooseFEM::Vector::dofs, "Return degrees-of-freedom") .def( "AsDofs", py::overload_cast&>(&GooseFEM::Vector::AsDofs, py::const_), "Set 'dofval", py::arg("nodevec")) .def( "AsDofs", py::overload_cast&>(&GooseFEM::Vector::AsDofs, py::const_), "Set 'dofval", py::arg("elemvec")) .def( "AsNode", py::overload_cast&>(&GooseFEM::Vector::AsNode, py::const_), "Set 'nodevec", py::arg("dofval")) .def( "AsNode", py::overload_cast&>(&GooseFEM::Vector::AsNode, py::const_), "Set 'nodevec", py::arg("elemvec")) .def( "AsElement", py::overload_cast&>( &GooseFEM::Vector::AsElement, py::const_), "Set 'elemvec", py::arg("dofval")) .def( "AsElement", py::overload_cast&>( &GooseFEM::Vector::AsElement, py::const_), "Set 'elemvec", py::arg("nodevec")) .def( "AssembleDofs", py::overload_cast&>( &GooseFEM::Vector::AssembleDofs, py::const_), "Assemble 'dofval'", py::arg("nodevec")) .def( "AssembleDofs", py::overload_cast&>( &GooseFEM::Vector::AssembleDofs, py::const_), "Assemble 'dofval'", py::arg("elemvec")) .def( "AssembleNode", py::overload_cast&>( &GooseFEM::Vector::AssembleNode, py::const_), "Assemble 'nodevec'", py::arg("elemvec")) + .def( + "AllocateDofval", + py::overload_cast<>(&GooseFEM::Vector::AllocateDofval, py::const_)) + + .def( + "AllocateDofval", + py::overload_cast(&GooseFEM::Vector::AllocateDofval, py::const_)) + + .def( + "AllocateNodevec", + py::overload_cast<>(&GooseFEM::Vector::AllocateNodevec, py::const_)) + + .def( + "AllocateNodevec", + py::overload_cast(&GooseFEM::Vector::AllocateNodevec, py::const_)) + + .def( + "AllocateElemvec", + py::overload_cast<>(&GooseFEM::Vector::AllocateElemvec, py::const_)) + + .def( + "AllocateElemvec", + py::overload_cast(&GooseFEM::Vector::AllocateElemvec, py::const_)) + .def("__repr__", [](const GooseFEM::Vector&) { return ""; }); } diff --git a/python/VectorPartitioned.hpp b/python/VectorPartitioned.hpp index c94c288..9056929 100644 --- a/python/VectorPartitioned.hpp +++ b/python/VectorPartitioned.hpp @@ -1,200 +1,224 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_VectorPartitioned(py::module& m) { py::class_(m, "VectorPartitioned") .def( py::init< const xt::xtensor&, const xt::xtensor&, const xt::xtensor&>(), "Switch between dofval/nodevec/elemvec", py::arg("conn"), py::arg("dofs"), py::arg("iip")) .def("nelem", &GooseFEM::VectorPartitioned::nelem, "Number of element") .def("nne", &GooseFEM::VectorPartitioned::nne, "Number of nodes per element") .def("nnode", &GooseFEM::VectorPartitioned::nnode, "Number of nodes") .def("ndim", &GooseFEM::VectorPartitioned::ndim, "Number of dimensions") .def("ndof", &GooseFEM::VectorPartitioned::ndof, "Number of degrees-of-freedom") .def("nnu", &GooseFEM::VectorPartitioned::nnu, "Number of unknown degrees-of-freedom") .def("nnp", &GooseFEM::VectorPartitioned::nnp, "Number of prescribed degrees-of-freedom") .def("dofs", &GooseFEM::VectorPartitioned::dofs, "Return degrees-of-freedom") .def("iiu", &GooseFEM::VectorPartitioned::iiu, "Return unknown degrees-of-freedom") .def("iip", &GooseFEM::VectorPartitioned::iip, "Return prescribed degrees-of-freedom") .def( "AsDofs", py::overload_cast&, const xt::xtensor&>( &GooseFEM::VectorPartitioned::AsDofs, py::const_), "Set 'dofval", py::arg("dofval_u"), py::arg("dofval_p")) .def( "AsDofs", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsDofs, py::const_), "Set 'dofval", py::arg("nodevec")) .def( "AsDofs", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsDofs, py::const_), "Set 'dofval", py::arg("elemvec")) .def( "AsDofs_u", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsDofs_u, py::const_), "Set 'dofval", py::arg("nodevec")) .def( "AsDofs_u", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsDofs_u, py::const_), "Set 'dofval", py::arg("elemvec")) .def( "AsDofs_p", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsDofs_p, py::const_), "Set 'dofval", py::arg("nodevec")) .def( "AsDofs_p", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsDofs_p, py::const_), "Set 'dofval", py::arg("elemvec")) .def( "AsNode", py::overload_cast&, const xt::xtensor&>( &GooseFEM::VectorPartitioned::AsNode, py::const_), "Set 'nodevec", py::arg("dofval_u"), py::arg("dofval_p")) .def( "AsNode", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsNode, py::const_), "Set 'nodevec", py::arg("dofval")) .def( "AsNode", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsNode, py::const_), "Set 'nodevec", py::arg("elemvec")) .def( "AsElement", py::overload_cast&, const xt::xtensor&>( &GooseFEM::VectorPartitioned::AsElement, py::const_), "Set 'elemvec", py::arg("dofval_u"), py::arg("dofval_p")) .def( "AsElement", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsElement, py::const_), "Set 'elemvec", py::arg("dofval")) .def( "AsElement", py::overload_cast&>( &GooseFEM::VectorPartitioned::AsElement, py::const_), "Set 'elemvec", py::arg("nodevec")) .def( "AssembleDofs", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleDofs, py::const_), "Assemble 'dofval'", py::arg("nodevec")) .def( "AssembleDofs", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleDofs, py::const_), "Assemble 'dofval'", py::arg("elemvec")) .def( "AssembleDofs_u", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleDofs_u, py::const_), "Assemble 'dofval'", py::arg("nodevec")) .def( "AssembleDofs_u", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleDofs_u, py::const_), "Assemble 'dofval'", py::arg("elemvec")) .def( "AssembleDofs_p", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleDofs_p, py::const_), "Assemble 'dofval'", py::arg("nodevec")) .def( "AssembleDofs_p", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleDofs_p, py::const_), "Assemble 'dofval'", py::arg("elemvec")) .def( "AssembleNode", py::overload_cast&>( &GooseFEM::VectorPartitioned::AssembleNode, py::const_), "Assemble 'nodevec'", py::arg("elemvec")) .def("Copy", &GooseFEM::VectorPartitioned::Copy, "Copy") .def("Copy_u", &GooseFEM::VectorPartitioned::Copy_u, "Copy iiu") .def("Copy_p", &GooseFEM::VectorPartitioned::Copy_p, "Copy iip") + .def( + "AllocateDofval", + py::overload_cast<>(&GooseFEM::VectorPartitioned::AllocateDofval, py::const_)) + + .def( + "AllocateDofval", + py::overload_cast(&GooseFEM::VectorPartitioned::AllocateDofval, py::const_)) + + .def( + "AllocateNodevec", + py::overload_cast<>(&GooseFEM::VectorPartitioned::AllocateNodevec, py::const_)) + + .def( + "AllocateNodevec", + py::overload_cast(&GooseFEM::VectorPartitioned::AllocateNodevec, py::const_)) + + .def( + "AllocateElemvec", + py::overload_cast<>(&GooseFEM::VectorPartitioned::AllocateElemvec, py::const_)) + + .def( + "AllocateElemvec", + py::overload_cast(&GooseFEM::VectorPartitioned::AllocateElemvec, py::const_)) + .def("__repr__", [](const GooseFEM::VectorPartitioned&) { return ""; }); } diff --git a/python/VectorPartitionedTyings.hpp b/python/VectorPartitionedTyings.hpp index 6f9a336..c8e955d 100644 --- a/python/VectorPartitionedTyings.hpp +++ b/python/VectorPartitionedTyings.hpp @@ -1,90 +1,114 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include namespace py = pybind11; void init_VectorPartitionedTyings(py::module& m) { py::class_(m, "VectorPartitionedTyings") .def( py::init< const xt::xtensor&, const xt::xtensor&, const Eigen::SparseMatrix&, const Eigen::SparseMatrix&, const Eigen::SparseMatrix&>(), "Switch between dofval/nodevec/elemvec", py::arg("conn"), py::arg("dofs"), py::arg("Cdu"), py::arg("Cdp"), py::arg("Cdi")) .def("nelem", &GooseFEM::VectorPartitionedTyings::nelem, "Number of element") .def("nne", &GooseFEM::VectorPartitionedTyings::nne, "Number of nodes per element") .def("nnode", &GooseFEM::VectorPartitionedTyings::nnode, "Number of nodes") .def("ndim", &GooseFEM::VectorPartitionedTyings::ndim, "Number of dimensions") .def("ndof", &GooseFEM::VectorPartitionedTyings::ndof, "Number of DOFs") .def("nnu", &GooseFEM::VectorPartitionedTyings::nnu, "Number of unknown DOFs") .def("nnp", &GooseFEM::VectorPartitionedTyings::nnp, "Number of prescribed DOFs") .def("nni", &GooseFEM::VectorPartitionedTyings::nni, "Number of independent DOFs") .def("nnd", &GooseFEM::VectorPartitionedTyings::nnd, "Number of dependent DOFs") .def("dofs", &GooseFEM::VectorPartitionedTyings::dofs, "DOFs") .def("iiu", &GooseFEM::VectorPartitionedTyings::iiu, "Unknown DOFs") .def("iip", &GooseFEM::VectorPartitionedTyings::iip, "Prescribed DOFs") .def("iii", &GooseFEM::VectorPartitionedTyings::iii, "Independent DOFs") .def("iid", &GooseFEM::VectorPartitionedTyings::iid, "Dependent DOFs") .def( "AsDofs_i", &GooseFEM::VectorPartitionedTyings::AsDofs_i, "Set 'dofval", py::arg("nodevec")) .def( "AsNode", &GooseFEM::VectorPartitionedTyings::AsNode, "Set 'nodevec", py::arg("dofval")) .def( "AsElement", &GooseFEM::VectorPartitionedTyings::AsElement, "Set 'elemvec", py::arg("nodevec")) .def( "AssembleDofs", &GooseFEM::VectorPartitionedTyings::AssembleDofs, "Assemble 'dofval'", py::arg("elemvec")) .def( "AssembleNode", &GooseFEM::VectorPartitionedTyings::AssembleNode, "Assemble 'nodevec'", py::arg("elemvec")) + .def( + "AllocateDofval", + py::overload_cast<>(&GooseFEM::VectorPartitionedTyings::AllocateDofval, py::const_)) + + .def( + "AllocateDofval", + py::overload_cast(&GooseFEM::VectorPartitionedTyings::AllocateDofval, py::const_)) + + .def( + "AllocateNodevec", + py::overload_cast<>(&GooseFEM::VectorPartitionedTyings::AllocateNodevec, py::const_)) + + .def( + "AllocateNodevec", + py::overload_cast(&GooseFEM::VectorPartitionedTyings::AllocateNodevec, py::const_)) + + .def( + "AllocateElemvec", + py::overload_cast<>(&GooseFEM::VectorPartitionedTyings::AllocateElemvec, py::const_)) + + .def( + "AllocateElemvec", + py::overload_cast(&GooseFEM::VectorPartitionedTyings::AllocateElemvec, py::const_)) + .def("__repr__", [](const GooseFEM::VectorPartitionedTyings&) { return ""; }); } diff --git a/test/basic/ElementHex8.cpp b/test/basic/ElementHex8.cpp index 35c478b..4376c3a 100644 --- a/test/basic/ElementHex8.cpp +++ b/test/basic/ElementHex8.cpp @@ -1,173 +1,180 @@ #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::ElementHex8", "ElementHex8.h") { SECTION("dV - Gauss") { GooseFEM::Mesh::Hex8::Regular mesh(2, 2, 2); auto coor = mesh.coor(); auto top = mesh.nodesTop(); auto right = mesh.nodesRight(); auto back = mesh.nodesBack(); xt::view(coor, xt::keep(back), xt::keep(2)) += 1.0; xt::view(coor, xt::keep(top), xt::keep(1)) += 1.0; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.0; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Hex8::Quadrature quad(vec.AsElement(coor)); auto dV = quad.dV(); auto dV_tensor = quad.AsTensor<2>(dV); + REQUIRE(xt::has_shape(quad.AllocateQscalar(), dV.shape())); + REQUIRE(xt::has_shape(quad.AllocateQscalar(0.0), dV.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor<2>(), dV_tensor.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor<2>(0.0), dV_tensor.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor(2), dV_tensor.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor(2, 0.0), dV_tensor.shape())); + REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(4)), 2.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(5)), 4.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(6)), 4.0 / 8.0)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(7)), 8.0 / 8.0)); for (size_t e = 0; e < mesh.nelem(); ++e) { for (size_t q = 0; q < quad.nip(); ++q) { REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } } } SECTION("dV - Nodal") { GooseFEM::Mesh::Hex8::Regular mesh(2, 2, 2); auto coor = mesh.coor(); auto top = mesh.nodesTop(); auto right = mesh.nodesRight(); auto back = mesh.nodesBack(); xt::view(coor, xt::keep(back), xt::keep(2)) += 1.0; xt::view(coor, xt::keep(top), xt::keep(1)) += 1.0; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.0; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Hex8::Quadrature quad( vec.AsElement(coor), GooseFEM::Element::Hex8::Nodal::xi(), GooseFEM::Element::Hex8::Nodal::w()); auto dV = quad.dV(); auto dV_tensor = quad.AsTensor<2>(dV); REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(4)), 2. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(5)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(6)), 4. / 8.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(7)), 8. / 8.)); for (size_t e = 0; e < mesh.nelem(); ++e) { for (size_t q = 0; q < quad.nip(); ++q) { REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } } } SECTION("int_N_scalar_NT_dV") { GooseFEM::Mesh::Hex8::Regular mesh(3, 3, 3); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Hex8::Quadrature quad( vec.AsElement(mesh.coor()), GooseFEM::Element::Hex8::Nodal::xi(), GooseFEM::Element::Hex8::Nodal::w()); xt::xtensor rho = xt::ones({mesh.nelem(), quad.nip()}); mat.assemble(quad.Int_N_scalar_NT_dV(rho)); auto M = mat.AsDiagonal(); REQUIRE(M.size() == vec.ndof()); REQUIRE(xt::allclose(M, 1.)); } SECTION("symGradN_vector") { GooseFEM::Mesh::Hex8::FineLayer mesh(27, 27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Hex8::Quadrature quad(vec.AsElement(mesh.coor())); xt::xtensor F = xt::zeros({3, 3}); xt::xtensor EPS = xt::zeros({3, 3}); F(0, 1) = 0.1; EPS(0, 1) = 0.05; EPS(1, 0) = 0.05; auto coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) { for (size_t i = 0; i < F.shape()[0]; ++i) { for (size_t j = 0; j < F.shape()[1]; ++j) { disp(n, i) += F(i, j) * coor(n, j); } } } auto eps = quad.SymGradN_vector(vec.AsElement(disp)); auto dV = quad.AsTensor<2>(quad.dV()); auto epsbar = xt::average(eps, dV, {0, 1}); REQUIRE(eps.shape()[0] == mesh.nelem()); REQUIRE(eps.shape()[1] == quad.nip()); REQUIRE(eps.shape()[2] == mesh.ndim()); REQUIRE(eps.shape()[3] == mesh.ndim()); for (size_t e = 0; e < mesh.nelem(); ++e) { for (size_t q = 0; q < quad.nip(); ++q) { REQUIRE(xt::allclose(xt::view(eps, e, q), EPS)); } } REQUIRE(xt::allclose(epsbar, EPS)); } SECTION("symGradN_vector, int_gradN_dot_tensor2s_dV") { GooseFEM::Mesh::Hex8::FineLayer mesh(27, 27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Hex8::Quadrature quad(vec.AsElement(mesh.coor())); xt::xtensor F = xt::zeros({3, 3}); F(0, 1) = 0.1; auto coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) { for (size_t i = 0; i < F.shape()[0]; ++i) { for (size_t j = 0; j < F.shape()[1]; ++j) { disp(n, i) += F(i, j) * coor(n, j); } } } auto eps = quad.SymGradN_vector(vec.AsElement(disp)); auto Fi = vec.AssembleDofs(quad.Int_gradN_dot_tensor2_dV(eps)); REQUIRE(Fi.size() == vec.ndof()); REQUIRE(xt::allclose(Fi, 0.)); } } diff --git a/test/basic/ElementQuad4.cpp b/test/basic/ElementQuad4.cpp index 6087b3b..406487d 100644 --- a/test/basic/ElementQuad4.cpp +++ b/test/basic/ElementQuad4.cpp @@ -1,161 +1,168 @@ #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::ElementQuad4", "ElementQuad4.h") { SECTION("dV - Gauss") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); auto coor = mesh.coor(); auto top = mesh.nodesTopEdge(); auto right = mesh.nodesRightEdge(); xt::view(coor, xt::keep(top), xt::keep(1)) += 1.0; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.0; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(coor)); auto dV = quad.dV(); auto dV_tensor = quad.AsTensor<2>(dV); + REQUIRE(xt::has_shape(quad.AllocateQscalar(), dV.shape())); + REQUIRE(xt::has_shape(quad.AllocateQscalar(0.0), dV.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor<2>(), dV_tensor.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor<2>(0.0), dV_tensor.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor(2), dV_tensor.shape())); + REQUIRE(xt::has_shape(quad.AllocateQtensor(2, 0.0), dV_tensor.shape())); + REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 4.)); for (size_t e = 0; e < mesh.nelem(); ++e) { for (size_t q = 0; q < quad.nip(); ++q) { REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } } } SECTION("dV - Nodal") { GooseFEM::Mesh::Quad4::Regular mesh(2, 2); auto coor = mesh.coor(); auto top = mesh.nodesTopEdge(); auto right = mesh.nodesRightEdge(); xt::view(coor, xt::keep(top), xt::keep(1)) += 1.0; xt::view(coor, xt::keep(right), xt::keep(0)) += 1.0; GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Quad4::Quadrature quad( vec.AsElement(coor), GooseFEM::Element::Quad4::Nodal::xi(), GooseFEM::Element::Quad4::Nodal::w()); auto dV = quad.dV(); auto dV_tensor = quad.AsTensor<2>(dV); REQUIRE(xt::allclose(xt::view(dV, xt::keep(0)), 1. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(1)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(2)), 2. / 4.)); REQUIRE(xt::allclose(xt::view(dV, xt::keep(3)), 4. / 4.)); for (size_t e = 0; e < mesh.nelem(); ++e) { for (size_t q = 0; q < quad.nip(); ++q) { REQUIRE(xt::allclose(xt::view(dV_tensor, xt::keep(e), xt::keep(q)), dV(e, q))); } } } SECTION("int_N_scalar_NT_dV") { GooseFEM::Mesh::Quad4::Regular mesh(3, 3); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Quad4::Quadrature quad( vec.AsElement(mesh.coor()), GooseFEM::Element::Quad4::Nodal::xi(), GooseFEM::Element::Quad4::Nodal::w()); xt::xtensor rho = xt::ones({mesh.nelem(), quad.nip()}); mat.assemble(quad.Int_N_scalar_NT_dV(rho)); auto M = mat.AsDiagonal(); REQUIRE(M.size() == vec.ndof()); REQUIRE(xt::allclose(M, 1.)); } SECTION("symGradN_vector") { GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofs()); GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor())); xt::xtensor F = xt::zeros({2, 2}); xt::xtensor EPS = xt::zeros({2, 2}); F(0, 1) = 0.1; EPS(0, 1) = 0.05; EPS(1, 0) = 0.05; xt::xtensor coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) { for (size_t i = 0; i < F.shape()[0]; ++i) { for (size_t j = 0; j < F.shape()[1]; ++j) { disp(n, i) += F(i, j) * coor(n, j); } } } auto eps = quad.SymGradN_vector(vec.AsElement(disp)); auto dV = quad.AsTensor<2>(quad.dV()); auto epsbar = xt::average(eps, dV, {0, 1}); REQUIRE(eps.shape()[0] == mesh.nelem()); REQUIRE(eps.shape()[1] == quad.nip()); REQUIRE(eps.shape()[2] == mesh.ndim()); REQUIRE(eps.shape()[3] == mesh.ndim()); for (size_t e = 0; e < mesh.nelem(); ++e) { for (size_t q = 0; q < quad.nip(); ++q) { REQUIRE(xt::allclose(xt::view(eps, e, q), EPS)); } } REQUIRE(xt::allclose(epsbar, EPS)); } SECTION("symGradN_vector, int_gradN_dot_tensor2s_dV") { GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27); GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic()); GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor())); xt::xtensor F = xt::zeros({2, 2}); F(0, 1) = 0.1; auto coor = mesh.coor(); xt::xtensor disp = xt::zeros(coor.shape()); for (size_t n = 0; n < mesh.nnode(); ++n) { for (size_t i = 0; i < F.shape()[0]; ++i) { for (size_t j = 0; j < F.shape()[1]; ++j) { disp(n, i) += F(i, j) * coor(n, j); } } } auto eps = quad.SymGradN_vector(vec.AsElement(disp)); auto Fi = vec.AssembleDofs(quad.Int_gradN_dot_tensor2_dV(eps)); REQUIRE(Fi.size() == vec.ndof()); REQUIRE(xt::allclose(Fi, 0.)); } } diff --git a/test/basic/Vector.cpp b/test/basic/Vector.cpp index 9d5d7e7..ff4e99f 100644 --- a/test/basic/Vector.cpp +++ b/test/basic/Vector.cpp @@ -1,203 +1,169 @@ #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Vector", "Vector.h") { SECTION("asDofs - nodevec") { - // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); - - // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); - // velocity field - // - allocate xt::xtensor v = xt::empty({mesh.nnode(), std::size_t(2)}); - // - set periodic v(0, 0) = 1.0; v(0, 1) = 0.0; v(1, 0) = 1.0; v(1, 1) = 0.0; v(2, 0) = 1.0; v(2, 1) = 0.0; v(3, 0) = 1.5; v(3, 1) = 0.0; v(4, 0) = 1.5; v(4, 1) = 0.0; v(5, 0) = 1.5; v(5, 1) = 0.0; v(6, 0) = 1.0; v(6, 1) = 0.0; v(7, 0) = 1.0; v(7, 1) = 0.0; v(8, 0) = 1.0; v(8, 1) = 0.0; - // convert to DOFs - xt::xtensor V = vector.AsDofs(v); + auto V = vector.AsDofs(v); - // check - // - size + REQUIRE(xt::has_shape(vector.AllocateDofval(), V.shape())); + REQUIRE(xt::has_shape(vector.AllocateNodevec(), v.shape())); + REQUIRE(xt::has_shape(vector.AllocateDofval(), vector.AsDofs(v).shape())); + REQUIRE(xt::has_shape(vector.AllocateElemvec(), vector.AsElement(v).shape())); + REQUIRE(xt::has_shape(vector.AllocateNodevec(), vector.AsNode(V).shape())); + REQUIRE(xt::has_shape(vector.AllocateElemvec(), vector.AsElement(V).shape())); REQUIRE(V.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); - // - individual entries ISCLOSE(V(0), v(0, 0)); ISCLOSE(V(1), v(0, 1)); ISCLOSE(V(2), v(1, 0)); ISCLOSE(V(3), v(1, 1)); ISCLOSE(V(4), v(3, 0)); ISCLOSE(V(5), v(3, 1)); ISCLOSE(V(6), v(4, 0)); ISCLOSE(V(7), v(4, 1)); } SECTION("asDofs - elemvec") { - // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); - - // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); - // velocity field - // - allocate xt::xtensor v = xt::empty({mesh.nnode(), std::size_t(2)}); - // - set periodic v(0, 0) = 1.0; v(0, 1) = 0.0; v(1, 0) = 1.0; v(1, 1) = 0.0; v(2, 0) = 1.0; v(2, 1) = 0.0; v(3, 0) = 1.5; v(3, 1) = 0.0; v(4, 0) = 1.5; v(4, 1) = 0.0; v(5, 0) = 1.5; v(5, 1) = 0.0; v(6, 0) = 1.0; v(6, 1) = 0.0; v(7, 0) = 1.0; v(7, 1) = 0.0; v(8, 0) = 1.0; v(8, 1) = 0.0; - // convert to DOFs - element - DOFs - xt::xtensor V = vector.AsDofs(vector.AsElement(vector.AsDofs(v))); + auto V = vector.AsDofs(vector.AsElement(vector.AsDofs(v))); - // check - // - size REQUIRE(V.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); - // - individual entries ISCLOSE(V(0), v(0, 0)); ISCLOSE(V(1), v(0, 1)); ISCLOSE(V(2), v(1, 0)); ISCLOSE(V(3), v(1, 1)); ISCLOSE(V(4), v(3, 0)); ISCLOSE(V(5), v(3, 1)); ISCLOSE(V(6), v(4, 0)); ISCLOSE(V(7), v(4, 1)); } SECTION("asDofs - assembleDofs") { - // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); - - // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); - // force field - // - allocate xt::xtensor f = xt::empty({mesh.nnode(), std::size_t(2)}); - // - set periodic f(0, 0) = -1.0; f(0, 1) = -1.0; f(1, 0) = 0.0; f(1, 1) = -1.0; f(2, 0) = 1.0; f(2, 1) = -1.0; f(3, 0) = -1.0; f(3, 1) = 0.0; f(4, 0) = 0.0; f(4, 1) = 0.0; f(5, 0) = 1.0; f(5, 1) = 0.0; f(6, 0) = -1.0; f(6, 1) = 1.0; f(7, 0) = 0.0; f(7, 1) = 1.0; f(8, 0) = 1.0; f(8, 1) = 1.0; - // assemble as DOFs xt::xtensor F = vector.AssembleDofs(f); - // check - // - size REQUIRE(F.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); - // - 'analytical' result ISCLOSE(F(0), 0); ISCLOSE(F(1), 0); ISCLOSE(F(2), 0); ISCLOSE(F(3), 0); ISCLOSE(F(4), 0); ISCLOSE(F(5), 0); ISCLOSE(F(6), 0); ISCLOSE(F(7), 0); } SECTION("asDofs - assembleNode") { - // mesh GooseFEM::Mesh::Quad4::Regular mesh(2, 2); - - // vector-definition GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic()); - // force field - // - allocate xt::xtensor f = xt::empty({mesh.nnode(), std::size_t(2)}); - // - set periodic f(0, 0) = -1.0; f(0, 1) = -1.0; f(1, 0) = 0.0; f(1, 1) = -1.0; f(2, 0) = 1.0; f(2, 1) = -1.0; f(3, 0) = -1.0; f(3, 1) = 0.0; f(4, 0) = 0.0; f(4, 1) = 0.0; f(5, 0) = 1.0; f(5, 1) = 0.0; f(6, 0) = -1.0; f(6, 1) = 1.0; f(7, 0) = 0.0; f(7, 1) = 1.0; f(8, 0) = 1.0; f(8, 1) = 1.0; - // convert to element, assemble as DOFs xt::xtensor F = vector.AssembleDofs(vector.AsElement(f)); - // check - // - size REQUIRE(F.size() == (mesh.nnode() - mesh.nodesPeriodic().shape(0)) * mesh.ndim()); - // - 'analytical' result ISCLOSE(F(0), 0); ISCLOSE(F(1), 0); ISCLOSE(F(2), 0); ISCLOSE(F(3), 0); ISCLOSE(F(4), 0); ISCLOSE(F(5), 0); ISCLOSE(F(6), 0); ISCLOSE(F(7), 0); } } diff --git a/test/basic/VectorPartitioned.cpp b/test/basic/VectorPartitioned.cpp index a6faa47..9d06944 100644 --- a/test/basic/VectorPartitioned.cpp +++ b/test/basic/VectorPartitioned.cpp @@ -1,131 +1,68 @@ #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::VectorPartitioned", "VectorPartitioned.h") { SECTION("asDofs") { - // mesh - GooseFEM::Mesh::Quad4::Regular mesh(9, 9); - // prescribed DOFs - - xt::xtensor dofs = mesh.dofs(); - - xt::xtensor nodesLeft = mesh.nodesLeftOpenEdge(); - xt::xtensor nodesRight = mesh.nodesRightOpenEdge(); - xt::xtensor nodesTop = mesh.nodesTopEdge(); - xt::xtensor nodesBottom = mesh.nodesBottomEdge(); + auto dofs = mesh.dofs(); + auto nodesLeft = mesh.nodesLeftOpenEdge(); + auto nodesRight = mesh.nodesRightOpenEdge(); + auto nodesTop = mesh.nodesTopEdge(); + auto nodesBottom = mesh.nodesBottomEdge(); - for (size_t i = 0; i < nodesLeft.size(); ++i) - for (size_t j = 0; j < mesh.ndim(); ++j) - dofs(nodesRight(i), j) = dofs(nodesLeft(i), j); + xt::view(dofs, xt::keep(nodesRight)) = xt::view(dofs, xt::keep(nodesLeft)); - xt::xtensor iip = xt::empty({4 * nodesBottom.size()}); - - size_t i = 0; - for (auto& n : nodesTop) { - iip(i) = dofs(n, 0); - ++i; - } - for (auto& n : nodesTop) { - iip(i) = dofs(n, 1); - ++i; - } - for (auto& n : nodesBottom) { - iip(i) = dofs(n, 0); - ++i; - } - for (auto& n : nodesBottom) { - iip(i) = dofs(n, 1); - ++i; - } - - // random displacement + size_t ni = nodesBottom.size(); + xt::xtensor iip = xt::empty({4 * ni}); + xt::view(iip, xt::range(0 * ni, 1 * ni)) = xt::view(dofs, xt::keep(nodesTop), 0); + xt::view(iip, xt::range(1 * ni, 2 * ni)) = xt::view(dofs, xt::keep(nodesTop), 1); + xt::view(iip, xt::range(2 * ni, 3 * ni)) = xt::view(dofs, xt::keep(nodesBottom), 0); + xt::view(iip, xt::range(3 * ni, 4 * ni)) = xt::view(dofs, xt::keep(nodesBottom), 1); xt::xtensor u = xt::random::randn({mesh.nnode(), mesh.ndim()}); - - for (size_t i = 0; i < nodesLeft.size(); ++i) - for (size_t j = 0; j < mesh.ndim(); ++j) - u(nodesRight(i), j) = u(nodesLeft(i), j); - - // vector-definition + xt::view(u, xt::keep(nodesRight)) = xt::view(u, xt::keep(nodesLeft)); GooseFEM::VectorPartitioned vector(mesh.conn(), dofs, iip); - - // convert - - xt::xtensor u_u = vector.AsDofs_u(u); - xt::xtensor u_p = vector.AsDofs_p(u); - + auto u_u = vector.AsDofs_u(u); + auto u_p = vector.AsDofs_p(u); REQUIRE(xt::allclose(u, vector.AsNode(u_u, u_p))); } SECTION("copy_u, copy_p") { - // mesh - GooseFEM::Mesh::Quad4::Regular mesh(9, 9); - // prescribed DOFs + auto dofs = mesh.dofs(); + auto nodesLeft = mesh.nodesLeftOpenEdge(); + auto nodesRight = mesh.nodesRightOpenEdge(); + auto nodesTop = mesh.nodesTopEdge(); + auto nodesBottom = mesh.nodesBottomEdge(); - xt::xtensor dofs = mesh.dofs(); + xt::view(dofs, xt::keep(nodesRight)) = xt::view(dofs, xt::keep(nodesLeft)); - xt::xtensor nodesLeft = mesh.nodesLeftOpenEdge(); - xt::xtensor nodesRight = mesh.nodesRightOpenEdge(); - xt::xtensor nodesTop = mesh.nodesTopEdge(); - xt::xtensor nodesBottom = mesh.nodesBottomEdge(); - - for (size_t i = 0; i < nodesLeft.size(); ++i) - for (size_t j = 0; j < mesh.ndim(); ++j) - dofs(nodesRight(i), j) = dofs(nodesLeft(i), j); - - xt::xtensor iip = xt::empty({4 * nodesBottom.size()}); - - size_t i = 0; - for (auto& n : nodesTop) { - iip(i) = dofs(n, 0); - ++i; - } - for (auto& n : nodesTop) { - iip(i) = dofs(n, 1); - ++i; - } - for (auto& n : nodesBottom) { - iip(i) = dofs(n, 0); - ++i; - } - for (auto& n : nodesBottom) { - iip(i) = dofs(n, 1); - ++i; - } - - // random displacement + size_t ni = nodesBottom.size(); + xt::xtensor iip = xt::empty({4 * ni}); + xt::view(iip, xt::range(0 * ni, 1 * ni)) = xt::view(dofs, xt::keep(nodesTop), 0); + xt::view(iip, xt::range(1 * ni, 2 * ni)) = xt::view(dofs, xt::keep(nodesTop), 1); + xt::view(iip, xt::range(2 * ni, 3 * ni)) = xt::view(dofs, xt::keep(nodesBottom), 0); + xt::view(iip, xt::range(3 * ni, 4 * ni)) = xt::view(dofs, xt::keep(nodesBottom), 1); xt::xtensor u = xt::random::randn({mesh.nnode(), mesh.ndim()}); - - for (size_t i = 0; i < nodesLeft.size(); ++i) - for (size_t j = 0; j < mesh.ndim(); ++j) - u(nodesRight(i), j) = u(nodesLeft(i), j); - - // vector-definition + xt::view(u, xt::keep(nodesRight)) = xt::view(u, xt::keep(nodesLeft)); GooseFEM::VectorPartitioned vector(mesh.conn(), dofs, iip); - - // convert - xt::xtensor v = xt::empty({mesh.nnode(), mesh.ndim()}); - vector.copy_u(u, v); vector.copy_p(u, v); - REQUIRE(xt::allclose(u, v)); } }