diff --git a/docs/GooseFEM.svg b/docs/GooseFEM.svg
new file mode 100644
index 0000000..e786884
--- /dev/null
+++ b/docs/GooseFEM.svg
@@ -0,0 +1,293 @@
+
+
+
+
diff --git a/include/xGooseFEM/Dynamics.h b/include/xGooseFEM/Dynamics.h
index 8578a00..66afee2 100644
--- a/include/xGooseFEM/Dynamics.h
+++ b/include/xGooseFEM/Dynamics.h
@@ -1,64 +1,67 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef XGOOSEFEM_DYNAMICS_H
#define XGOOSEFEM_DYNAMICS_H
// -------------------------------------------------------------------------------------------------
#include "GooseFEM.h"
// ======================================= GooseFEM::Dynamics =======================================
namespace xGooseFEM {
namespace Dynamics {
// ------------------------------------------ dummy class ------------------------------------------
class Geometry
{
public:
+ // destructor
+ virtual ~Geometry() {};
+
// solve for DOF-accelerations [ndof]
virtual xt::xtensor solve_A() { return xt::empty({0}); };
virtual xt::xtensor solve_V() { return xt::empty({0}); };
// return nodal vectors [nnode, ndim]
virtual xt::xtensor u() const { return xt::empty({0,0}); };
virtual xt::xtensor v() const { return xt::empty({0,0}); };
virtual xt::xtensor a() const { return xt::empty({0,0}); };
// return DOF values [ndof]
virtual xt::xtensor dofs_u() const { return xt::empty({0}); };
virtual xt::xtensor dofs_v() const { return xt::empty({0}); };
virtual xt::xtensor dofs_a() const { return xt::empty({0}); };
// overwrite nodal vectors [nnode, ndim]
virtual void set_u(const xt::xtensor &nodevec) { UNUSED(nodevec); return; };
// overwrite nodal vectors, reconstructed from DOF values [ndof]
virtual void set_u(const xt::xtensor &dofval) { UNUSED(dofval); return; };
virtual void set_v(const xt::xtensor &dofval) { UNUSED(dofval); return; };
virtual void set_a(const xt::xtensor &dofval) { UNUSED(dofval); return; };
};
// ------------------------------------ evaluate one time step -------------------------------------
inline void Verlet (Geometry &geometry, double dt, size_t nstep=1);
inline void velocityVerlet(Geometry &geometry, double dt, size_t nstep=1);
namespace Overdamped {
inline void forwardEuler (Geometry &geometry, double dt, size_t nstep=1);
}
// -------------------------------------------------------------------------------------------------
}} // namespace ...
// =================================================================================================
#endif
diff --git a/include/xGooseFEM/ElementQuad4.h b/include/xGooseFEM/ElementQuad4.h
index 099dde1..1bac6d1 100644
--- a/include/xGooseFEM/ElementQuad4.h
+++ b/include/xGooseFEM/ElementQuad4.h
@@ -1,132 +1,156 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef XGOOSEFEM_ELEMENTQUAD4_H
#define XGOOSEFEM_ELEMENTQUAD4_H
// -------------------------------------------------------------------------------------------------
#include "GooseFEM.h"
// =================================== GooseFEM::Element::Quad4 ====================================
namespace xGooseFEM {
namespace Element {
namespace Quad4 {
// ======================================== tensor algebra =========================================
using T2 = xt::xtensor_fixed>;
inline double inv(const T2 &A, T2 &Ainv);
// ================================ GooseFEM::Element::Quad4::Gauss ================================
namespace Gauss {
inline size_t nip(); // number of integration points
inline xt::xtensor xi(); // integration point coordinates (local coordinates)
inline xt::xtensor w(); // integration point weights
}
// ================================ GooseFEM::Element::Quad4::Nodal ================================
namespace Nodal {
inline size_t nip(); // number of integration points
inline xt::xtensor xi(); // integration point coordinates (local coordinates)
inline xt::xtensor w(); // integration point weights
}
// =================================================================================================
// ------------------------------------------ quadrature -------------------------------------------
class Quadrature
{
private:
// dimensions (flexible)
size_t m_nelem; // number of elements
size_t m_nip; // number of integration points
// dimensions (fixed for this element type)
static const size_t m_nne=4; // number of nodes per element
static const size_t m_ndim=2; // number of dimensions
// data arrays
xt::xtensor m_x; // nodal positions stored per element [nelem, nne, ndim]
xt::xtensor m_w; // weight of each integration point [nip]
xt::xtensor m_xi; // local coordinate of each integration point [nip, ndim]
xt::xtensor m_N; // shape functions [nip, nne]
xt::xtensor m_dNxi; // shape function gradients w.r.t. local coordinate [nip, nne, ndim]
xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim]
xt::xtensor m_vol; // integration point volume [nelem, nip]
private:
// compute "vol" and "dNdx" based on current "x"
void compute_dN();
public:
// convention:
// "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim]
// "elemvec" - nodal vectors stored per element - [nelem, nne, ndim]
// "qtensor" - integration point tensor - [nelem, nip, ndim, ndim]
// "qscalar" - integration point scalar - [nelem, nip]
// constructor: integration point coordinates and weights are optional (default: Gauss)
+
Quadrature() = default;
+
Quadrature(const xt::xtensor &x);
+
Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w);
// update the nodal positions (shape of "x" should match the earlier definition)
+
void update_x(const xt::xtensor &x);
// return dimensions
+
size_t nelem() const; // number of elements
size_t nne() const; // number of nodes per element
size_t ndim() const; // number of dimension
size_t nip() const; // number of integration points
// return integration volume
- // - in-place
+
void dV(xt::xtensor &qscalar) const;
void dV(xt::xtensor &qtensor) const; // same volume for all tensor components
- // - return qscalar/qtensor
+
+ // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part
+
+ void gradN_vector (const xt::xtensor &elemvec,
+ xt::xtensor &qtensor) const;
+
+ void gradN_vector_T (const xt::xtensor &elemvec,
+ xt::xtensor &qtensor) const;
+
+ void symGradN_vector(const xt::xtensor &elemvec,
+ xt::xtensor &qtensor) const;
+
+ // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV"
+
+ void int_N_scalar_NT_dV(const xt::xtensor &qscalar,
+ xt::xtensor &elemmat) const;
+
+ // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV"
+
+ void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor,
+ xt::xtensor &elemvec) const;
+
+ // integral of the dot product "elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV"
+
+ void int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor,
+ xt::xtensor &elemmat) const;
+
+ // auto allocation of the functions above
+
xt::xtensor dV() const;
+
xt::xtensor dVtensor() const;
- // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part
- // - in-place
- void gradN_vector (const xt::xtensor &elemvec, xt::xtensor &qtensor) const;
- void gradN_vector_T (const xt::xtensor &elemvec, xt::xtensor &qtensor) const;
- void symGradN_vector(const xt::xtensor &elemvec, xt::xtensor &qtensor) const;
- // - return qtensor
- xt::xtensor gradN_vector (const xt::xtensor &elemvec) const;
- xt::xtensor gradN_vector_T (const xt::xtensor &elemvec) const;
+ xt::xtensor gradN_vector(const xt::xtensor &elemvec) const;
+
+ xt::xtensor gradN_vector_T(const xt::xtensor &elemvec) const;
+
xt::xtensor symGradN_vector(const xt::xtensor &elemvec) const;
- // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV"
- // - in-place
- void int_N_scalar_NT_dV(const xt::xtensor &qscalar, xt::xtensor &elemmat) const;
- // - return elemmat
xt::xtensor int_N_scalar_NT_dV(const xt::xtensor &qscalar) const;
- // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV"
- // - in-place
- void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const;
- // - return elemvec
xt::xtensor int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const;
+ xt::xtensor int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor) const;
+
};
// -------------------------------------------------------------------------------------------------
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/include/xGooseFEM/ElementQuad4.hpp b/include/xGooseFEM/ElementQuad4.hpp
index 8bf8611..1550a14 100644
--- a/include/xGooseFEM/ElementQuad4.hpp
+++ b/include/xGooseFEM/ElementQuad4.hpp
@@ -1,642 +1,702 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef XGOOSEFEM_ELEMENTQUAD4_CPP
#define XGOOSEFEM_ELEMENTQUAD4_CPP
// -------------------------------------------------------------------------------------------------
#include "ElementQuad4.h"
-// =================================== GooseFEM::Element::Quad4 ====================================
+// =================================================================================================
namespace xGooseFEM {
namespace Element {
namespace Quad4 {
-// ======================================== tensor algebra =========================================
+// =================================================================================================
inline double inv(const T2 &A, T2 &Ainv)
{
// compute determinant
double det = A[0] * A[3] - A[1] * A[2];
// compute inverse
Ainv[0] = A[3] / det;
Ainv[1] = -1. * A[1] / det;
Ainv[2] = -1. * A[2] / det;
Ainv[3] = A[0] / det;
return det;
}
-// ================================ GooseFEM::Element::Quad4::Gauss ================================
+// =================================================================================================
namespace Gauss {
-// --------------------------------- number of integration points ----------------------------------
+// -------------------------------------------------------------------------------------------------
inline size_t nip()
{
return 4;
}
-// ----------------------- integration point coordinates (local coordinates) -----------------------
+// -------------------------------------------------------------------------------------------------
inline xt::xtensor xi()
{
size_t nip = 4;
size_t ndim = 2;
xt::xtensor xi = xt::empty({nip,ndim});
xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.);
xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.);
xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.);
xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.);
return xi;
}
-// ----------------------------------- integration point weights -----------------------------------
+// -------------------------------------------------------------------------------------------------
inline xt::xtensor w()
{
size_t nip = 4;
xt::xtensor w = xt::empty({nip});
w(0) = 1.;
w(1) = 1.;
w(2) = 1.;
w(3) = 1.;
return w;
}
// -------------------------------------------------------------------------------------------------
}
-// ================================ GooseFEM::Element::Quad4::Nodal ================================
+// =================================================================================================
namespace Nodal {
-// --------------------------------- number of integration points ----------------------------------
+// -------------------------------------------------------------------------------------------------
inline size_t nip()
{
return 4;
}
-// ----------------------- integration point coordinates (local coordinates) -----------------------
+// -------------------------------------------------------------------------------------------------
inline xt::xtensor xi()
{
size_t nip = 4;
size_t ndim = 2;
xt::xtensor xi = xt::empty({nip,ndim});
xi(0,0) = -1.; xi(0,1) = -1.;
xi(1,0) = +1.; xi(1,1) = -1.;
xi(2,0) = +1.; xi(2,1) = +1.;
xi(3,0) = -1.; xi(3,1) = +1.;
return xi;
}
-// ----------------------------------- integration point weights -----------------------------------
+// -------------------------------------------------------------------------------------------------
inline xt::xtensor w()
{
size_t nip = 4;
xt::xtensor w = xt::empty({nip});
w(0) = 1.;
w(1) = 1.;
w(2) = 1.;
w(3) = 1.;
return w;
}
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
-// ------------------------------------------ constructor ------------------------------------------
+// -------------------------------------------------------------------------------------------------
inline Quadrature::Quadrature(const xt::xtensor &x) : m_x(x)
{
assert( m_x.shape()[1] == m_nne );
assert( m_x.shape()[2] == m_ndim );
// integration scheme
m_nip = Gauss::nip();
m_xi = Gauss::xi();
m_w = Gauss::w();
// extract number of elements
m_nelem = m_x.shape()[0];
// allocate arrays
m_N = xt::empty({ m_nip, m_nne });
m_dNxi = xt::empty({ m_nip, m_nne, m_ndim});
m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim});
m_vol = xt::empty({m_nelem, m_nip });
// shape functions
for ( size_t k = 0 ; k < m_nip ; ++k )
{
m_N(k,0) = .25 * (1.-m_xi(k,0)) * (1.-m_xi(k,1));
m_N(k,1) = .25 * (1.+m_xi(k,0)) * (1.-m_xi(k,1));
m_N(k,2) = .25 * (1.+m_xi(k,0)) * (1.+m_xi(k,1));
m_N(k,3) = .25 * (1.-m_xi(k,0)) * (1.+m_xi(k,1));
}
// shape function gradients in local coordinates
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - dN / dxi_0
m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1));
m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1));
m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1));
m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1));
// - dN / dxi_1
m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0));
m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0));
m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0));
m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0));
}
// compute the shape function gradients, based on "x"
compute_dN();
}
-// ------------------------------------------ constructor ------------------------------------------
+// -------------------------------------------------------------------------------------------------
inline Quadrature::Quadrature(const xt::xtensor &x, const xt::xtensor &xi,
const xt::xtensor &w) : m_x(x), m_w(w), m_xi(xi)
{
assert( m_x.shape()[1] == m_nne );
assert( m_x.shape()[2] == m_ndim );
// extract number of elements and number of integration points
m_nelem = m_x.shape()[0];
m_nip = m_w.size();
assert( m_xi.shape()[0] == m_nip );
assert( m_xi.shape()[1] == m_ndim );
assert( m_w .size() == m_nip );
// allocate arrays
m_N = xt::empty({ m_nip, m_nne });
m_dNxi = xt::empty({ m_nip, m_nne, m_ndim});
m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim});
m_vol = xt::empty({m_nelem, m_nip });
// shape functions
for ( size_t k = 0 ; k < m_nip ; ++k )
{
m_N(k,0) = .25 * (1.-m_xi(k,0)) * (1.-m_xi(k,1));
m_N(k,1) = .25 * (1.+m_xi(k,0)) * (1.-m_xi(k,1));
m_N(k,2) = .25 * (1.+m_xi(k,0)) * (1.+m_xi(k,1));
m_N(k,3) = .25 * (1.-m_xi(k,0)) * (1.+m_xi(k,1));
}
// shape function gradients in local coordinates
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - dN / dxi_0
m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1));
m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1));
m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1));
m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1));
// - dN / dxi_1
m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0));
m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0));
m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0));
m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0));
}
// compute the shape function gradients, based on "x"
compute_dN();
}
-// --------------------------- integration volume (per tensor-component) ---------------------------
+// -------------------------------------------------------------------------------------------------
inline void Quadrature::dV(xt::xtensor &qscalar) const
{
assert( qscalar.shape()[0] == m_nelem );
assert( qscalar.shape()[1] == m_nip );
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
for ( size_t k = 0 ; k < m_nip ; ++k )
qscalar(e,k) = m_vol(e,k);
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::dV(xt::xtensor &qtensor) const
{
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
for ( size_t k = 0 ; k < m_nip ; ++k )
for ( size_t i = 0 ; i < m_ndim ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
qtensor(e,k,i,j) = m_vol(e,k);
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor Quadrature::dV() const
-{
- xt::xtensor out = xt::empty({m_nelem, m_nip});
-
- this->dV(out);
-
- return out;
-}
-
-// -------------------------------------------------------------------------------------------------
-
-inline xt::xtensor Quadrature::dVtensor() const
-{
- xt::xtensor out = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
-
- this->dV(out);
-
- return out;
-}
-
-// -------------------------------------- number of elements ---------------------------------------
-
inline size_t Quadrature::nelem() const
{
return m_nelem;
}
-// ---------------------------------- number of nodes per element ----------------------------------
+// -------------------------------------------------------------------------------------------------
inline size_t Quadrature::nne() const
{
return m_nne;
}
-// ------------------------------------- number of dimensions --------------------------------------
+// -------------------------------------------------------------------------------------------------
inline size_t Quadrature::ndim() const
{
return m_ndim;
}
-// --------------------------------- number of integration points ----------------------------------
+// -------------------------------------------------------------------------------------------------
inline size_t Quadrature::nip() const
{
return m_nip;
}
-// --------------------------------------- update positions ----------------------------------------
+// -------------------------------------------------------------------------------------------------
inline void Quadrature::update_x(const xt::xtensor &x)
{
assert( x.shape()[0] == m_nelem );
assert( x.shape()[1] == m_nne );
assert( x.shape()[2] == m_ndim );
assert( x.size() == m_x.size() );
// update positions
xt::noalias(m_x) = x;
// update the shape function gradients for the new "x"
compute_dN();
}
-// ------------------------ shape function gradients in global coordinates -------------------------
+// -------------------------------------------------------------------------------------------------
inline void Quadrature::compute_dN()
{
// loop over all elements (in parallel)
#pragma omp parallel
{
// allocate local variables
T2 J, Jinv;
#pragma omp for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias nodal positions
auto x = xt::adapt(&m_x(e,0,0), xt::xshape());
// loop over integration points
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - alias
auto dNxi = xt::adapt(&m_dNxi( k,0,0), xt::xshape());
auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape());
// - Jacobian (loops unrolled for efficiency)
// J(i,j) += dNxi(m,i) * x(m,j);
J(0,0) = dNxi(0,0)*x(0,0) + dNxi(1,0)*x(1,0) + dNxi(2,0)*x(2,0) + dNxi(3,0)*x(3,0);
J(0,1) = dNxi(0,0)*x(0,1) + dNxi(1,0)*x(1,1) + dNxi(2,0)*x(2,1) + dNxi(3,0)*x(3,1);
J(1,0) = dNxi(0,1)*x(0,0) + dNxi(1,1)*x(1,0) + dNxi(2,1)*x(2,0) + dNxi(3,1)*x(3,0);
J(1,1) = dNxi(0,1)*x(0,1) + dNxi(1,1)*x(1,1) + dNxi(2,1)*x(2,1) + dNxi(3,1)*x(3,1);
// - determinant and inverse of the Jacobian
double Jdet = inv(J, Jinv);
// - shape function gradients wrt global coordinates (loops partly unrolled for efficiency)
// dNx(m,i) += Jinv(i,j) * dNxi(m,j);
for ( size_t m = 0 ; m < m_nne ; ++m )
{
dNx(m,0) = Jinv(0,0) * dNxi(m,0) + Jinv(0,1) * dNxi(m,1);
dNx(m,1) = Jinv(1,0) * dNxi(m,0) + Jinv(1,1) * dNxi(m,1);
}
// - integration point volume
m_vol(e,k) = m_w(k) * Jdet;
}
}
}
}
-// ------------------- dyadic product "qtensor(i,j) = dNdx(m,i) * elemvec(m,j)" --------------------
+// -------------------------------------------------------------------------------------------------
inline void Quadrature::gradN_vector(
const xt::xtensor &elemvec, xt::xtensor &qtensor) const
{
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
// zero-initialize output: matrix of tensors
qtensor.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias element vector (e.g. nodal displacements)
auto u = xt::adapt(&elemvec(e,0,0), xt::xshape());
// loop over all integration points in element "e"
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape());
auto gradu = xt::adapt(&qtensor(e,k,0,0), xt::xshape());
// - evaluate dyadic product (loops unrolled for efficiency)
// gradu(i,j) += dNx(m,i) * u(m,j)
gradu(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0);
gradu(0,1) = dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1);
gradu(1,0) = dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0);
gradu(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1);
}
}
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const
-{
- xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
-
- this->gradN_vector(elemvec, qtensor);
-
- return qtensor;
-}
-
-// ---------------------------------- transpose of "gradN_vector" ----------------------------------
-
inline void Quadrature::gradN_vector_T(
const xt::xtensor &elemvec, xt::xtensor &qtensor) const
{
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
// zero-initialize output: matrix of tensors
qtensor.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias element vector (e.g. nodal displacements)
auto u = xt::adapt(&elemvec(e,0,0), xt::xshape());
// loop over all integration points in element "e"
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape());
auto gradu = xt::adapt(&qtensor(e,k,0,0), xt::xshape());
// - evaluate transpose of dyadic product (loops unrolled for efficiency)
// gradu(j,i) += dNx(m,i) * u(m,j)
gradu(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0);
gradu(1,0) = dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1);
gradu(0,1) = dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0);
gradu(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1);
}
}
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const
-{
- xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
-
- this->gradN_vector_T(elemvec, qtensor);
-
- return qtensor;
-}
-
-// ------------------------------- symmetric part of "gradN_vector" --------------------------------
-
inline void Quadrature::symGradN_vector(
const xt::xtensor &elemvec, xt::xtensor &qtensor) const
{
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
// zero-initialize output: matrix of tensors
qtensor.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias element vector (e.g. nodal displacements)
auto u = xt::adapt(&elemvec(e,0,0), xt::xshape());
// loop over all integration points in element "e"
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape());
auto eps = xt::adapt(&qtensor(e,k,0,0), xt::xshape());
// - evaluate symmetrized dyadic product (loops unrolled for efficiency)
// grad(i,j) += dNx(m,i) * u(m,j)
// eps (j,i) = 0.5 * ( grad(i,j) + grad(j,i) )
eps(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0);
eps(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1);
eps(0,1) = ( dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1) +
dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0) ) / 2.;
eps(1,0) = eps(0,1);
}
}
}
-// -------------------------------------------------------------------------------------------------
-
-inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const
-{
- xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
-
- this->symGradN_vector(elemvec, qtensor);
-
- return qtensor;
-}
-
// ------- scalar product "elemmat(m*ndim+i,n*ndim+i) = N(m) * qscalar * N(n)"; for all "i" --------
inline void Quadrature::int_N_scalar_NT_dV(
const xt::xtensor &qscalar, xt::xtensor &elemmat) const
{
assert( qscalar.shape()[0] == m_nelem );
assert( qscalar.shape()[1] == m_nip );
assert( elemmat.shape()[0] == m_nelem );
assert( elemmat.shape()[1] == m_nne*m_ndim );
assert( elemmat.shape()[2] == m_nne*m_ndim );
// zero-initialize: matrix of matrices
elemmat.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias (e.g. mass matrix)
auto M = xt::adapt(&elemmat(e,0,0), xt::xshape());
// loop over all integration points in element "e"
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - alias
auto N = xt::adapt(&m_N(k,0), xt::xshape());
auto& vol = m_vol (e,k);
auto& rho = qscalar(e,k);
// - evaluate scalar product, for all dimensions, and assemble
// M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV
for ( size_t m = 0 ; m < m_nne ; ++m ) {
for ( size_t n = 0 ; n < m_nne ; ++n ) {
M(m*m_ndim+0, n*m_ndim+0) += N(m) * rho * N(n) * vol;
M(m*m_ndim+1, n*m_ndim+1) += N(m) * rho * N(n) * vol;
}
}
}
}
}
-// -------------------------------------------------------------------------------------------------
-
-inline xt::xtensor Quadrature::int_N_scalar_NT_dV(const xt::xtensor &qscalar) const
-{
- xt::xtensor elemmat = xt::empty({m_nelem, m_nne*m_ndim, m_nne*m_ndim});
-
- this->int_N_scalar_NT_dV(qscalar, elemmat);
-
- return elemmat;
-}
-
// ------------ integral of dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" ------------
inline void Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor,
xt::xtensor &elemvec) const
{
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nip );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
// zero-initialize output: matrix of vectors
elemvec.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias (e.g. nodal force)
auto f = xt::adapt(&elemvec(e,0,0), xt::xshape());
// loop over all integration points in element "e"
for ( size_t k = 0 ; k < m_nip ; ++k )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,k,0,0), xt::xshape());
auto sig = xt::adapt(&qtensor(e,k,0,0), xt::xshape());
auto& vol = m_vol(e,k);
// - evaluate dot product, and assemble
for ( size_t m = 0 ; m < m_nne ; ++m )
{
f(m,0) += ( dNx(m,0) * sig(0,0) + dNx(m,1) * sig(1,0) ) * vol;
f(m,1) += ( dNx(m,0) * sig(0,1) + dNx(m,1) * sig(1,1) ) * vol;
}
}
}
}
// -------------------------------------------------------------------------------------------------
-inline xt::xtensor Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const
+inline void Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor,
+ xt::xtensor &elemmat) const
+{
+ assert( qtensor.shape()[0] == m_nelem );
+ assert( qtensor.shape()[1] == m_nip );
+ assert( qtensor.shape()[2] == m_ndim );
+ assert( qtensor.shape()[3] == m_ndim );
+ assert( qtensor.shape()[4] == m_ndim );
+ assert( qtensor.shape()[5] == m_ndim );
+
+ assert( elemmat.shape()[0] == m_nelem );
+ assert( elemmat.shape()[1] == m_nne*m_ndim );
+ assert( elemmat.shape()[2] == m_nne*m_ndim );
+
+ // zero-initialize output: matrix of vector
+ elemmat.fill(0.0);
+
+ // loop over all elements (in parallel)
+ #pragma omp parallel for
+ for ( size_t e = 0 ; e < m_nelem ; ++e )
+ {
+ // alias (e.g. nodal force)
+ auto K = xt::adapt(&elemmat(e,0,0), xt::xshape());
+
+ // loop over all integration points in element "e"
+ for ( size_t q = 0 ; q < m_nip ; ++q ){
+
+ // - alias
+ auto dNx = xt::adapt(&m_dNx(e,q,0,0), xt::xshape());
+ auto C = xt::adapt(&qtensor(e,q,0,0,0,0), xt::xshape());
+ auto& vol = m_vol(e,q);
+
+ // - evaluate dot product, and assemble
+ for ( size_t m = 0 ; m < m_nne ; ++m )
+ for ( size_t n = 0 ; n < m_nne ; ++n )
+ for ( size_t i = 0 ; i < m_ndim ; ++i )
+ for ( size_t j = 0 ; j < m_ndim ; ++j )
+ for ( size_t k = 0 ; k < m_ndim ; ++k )
+ for ( size_t l = 0 ; l < m_ndim ; ++l )
+ K(m*m_ndim+j, n*m_ndim+k) += dNx(m,i) * C(i,j,k,l) * dNx(n,l) * vol;
+ }
+ }
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::dV() const
+{
+ xt::xtensor out = xt::empty({m_nelem, m_nip});
+
+ this->dV(out);
+
+ return out;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::dVtensor() const
+{
+ xt::xtensor out = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
+
+ this->dV(out);
+
+ return out;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const
+{
+ xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
+
+ this->gradN_vector(elemvec, qtensor);
+
+ return qtensor;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const
+{
+ xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
+
+ this->gradN_vector_T(elemvec, qtensor);
+
+ return qtensor;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const
+{
+ xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim});
+
+ this->symGradN_vector(elemvec, qtensor);
+
+ return qtensor;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::int_N_scalar_NT_dV(
+ const xt::xtensor &qscalar) const
+{
+ xt::xtensor elemmat = xt::empty({m_nelem, m_nne*m_ndim, m_nne*m_ndim});
+
+ this->int_N_scalar_NT_dV(qscalar, elemmat);
+
+ return elemmat;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+inline xt::xtensor Quadrature::int_gradN_dot_tensor2_dV(
+ const xt::xtensor &qtensor) const
{
xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim});
this->int_gradN_dot_tensor2_dV(qtensor, elemvec);
return elemvec;
}
// -------------------------------------------------------------------------------------------------
+inline xt::xtensor Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV(
+ const xt::xtensor &qtensor) const
+ {
+ xt::xtensor elemmat = xt::empty({m_nelem, m_ndim*m_nne, m_ndim*m_nne});
+
+ this->int_gradN_dot_tensor4_dot_gradNT_dV(qtensor, elemmat);
+
+ return elemmat;
+ }
+
+// -------------------------------------------------------------------------------------------------
+
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/include/xGooseFEM/GooseFEM.h b/include/xGooseFEM/GooseFEM.h
index af1a7c8..d885316 100644
--- a/include/xGooseFEM/GooseFEM.h
+++ b/include/xGooseFEM/GooseFEM.h
@@ -1,104 +1,111 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef XGOOSEFEM_H
#define XGOOSEFEM_H
// =================================================================================================
#define _USE_MATH_DEFINES // to use "M_PI" from "math.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace xt::placeholders;
// =================================================================================================
#define GOOSEFEM_WORLD_VERSION 0
#define GOOSEFEM_MAJOR_VERSION 1
#define GOOSEFEM_MINOR_VERSION 0
#define GOOSEFEM_VERSION_AT_LEAST(x,y,z) \
(GOOSEFEM_WORLD_VERSION>x || (GOOSEFEM_WORLD_VERSION>=x && \
(GOOSEFEM_MAJOR_VERSION>y || (GOOSEFEM_MAJOR_VERSION>=y && \
GOOSEFEM_MINOR_VERSION>=z))))
#define GOOSEFEM_VERSION(x,y,z) \
(GOOSEFEM_WORLD_VERSION==x && \
GOOSEFEM_MAJOR_VERSION==y && \
GOOSEFEM_MINOR_VERSION==z)
// =================================================================================================
// dummy operation that can be use to suppress the "unused parameter" warnings
#define UNUSED(p) ( (void)(p) )
// =================================================================================================
// alias Eigen sparse matrices
namespace xGooseFEM
{
typedef Eigen::SparseMatrix SpMatD;
typedef Eigen::SparseMatrix SpMatS;
+ typedef Eigen::Triplet TripD;
}
// =================================================================================================
#include "Mesh.h"
#include "MeshTri3.h"
#include "MeshQuad4.h"
#include "MeshHex8.h"
#include "Element.h"
#include "ElementQuad4.h"
#include "ElementHex8.h"
#include "Vector.h"
+#include "VectorPartitioned.h"
+#include "MatrixPartitioned.h"
#include "MatrixDiagonal.h"
+#include "MatrixDiagonalPartitioned.h"
#include "Iterate.h"
#include "Dynamics.h"
#include "Mesh.hpp"
#include "MeshTri3.hpp"
#include "MeshQuad4.hpp"
#include "MeshHex8.hpp"
#include "Element.hpp"
#include "ElementQuad4.hpp"
#include "ElementHex8.hpp"
#include "Vector.hpp"
+#include "VectorPartitioned.hpp"
+#include "MatrixPartitioned.hpp"
#include "MatrixDiagonal.hpp"
+#include "MatrixDiagonalPartitioned.hpp"
#include "Iterate.hpp"
#include "Dynamics.hpp"
// =================================================================================================
#endif
diff --git a/include/xGooseFEM/MatrixDiagonal.h b/include/xGooseFEM/MatrixDiagonal.h
index 862c400..0390489 100644
--- a/include/xGooseFEM/MatrixDiagonal.h
+++ b/include/xGooseFEM/MatrixDiagonal.h
@@ -1,147 +1,81 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef XGOOSEFEM_MATRIXDIAGONAL_H
#define XGOOSEFEM_MATRIXDIAGONAL_H
// -------------------------------------------------------------------------------------------------
#include "GooseFEM.h"
// =========================================== GooseFEM ============================================
namespace xGooseFEM {
// -------------------------------------------------------------------------------------------------
class MatrixDiagonal
{
public:
// constructor
MatrixDiagonal() = default;
-
MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs);
- MatrixDiagonal(const xt::xtensor &conn, const xt::xtensor &dofs,
- const xt::xtensor &iip);
-
- // index operators: access plain storage
- double& operator[](size_t i);
- const double& operator[](size_t i) const;
-
- // index operators: access using matrix indices
- double& operator()(size_t i);
- const double& operator()(size_t i) const;
- double& operator()(size_t i, size_t j);
- const double& operator()(size_t i, size_t j) const;
-
// dimensions
size_t nelem() const; // number of elements
size_t nne() const; // number of nodes per element
size_t nnode() const; // number of nodes
size_t ndim() const; // number of dimensions
size_t ndof() const; // number of DOFs
- size_t nnu() const; // number of unknown DOFs
- size_t nnp() const; // number of prescribed DOFs
// DOF lists
xt::xtensor dofs() const; // DOFs
- xt::xtensor iiu() const; // unknown DOFs
- xt::xtensor iip() const; // prescribed DOFs
// product: b_i = A_ij * x_j
- // b = A * x
- xt::xtensor dot (const xt::xtensor &x ) const;
- // b_u = A_uu * x_u + A_up * x_p
- xt::xtensor dot_u(const xt::xtensor &x_u, const xt::xtensor &x_p) const;
- // b_p = A_pu * x_u + A_pp * x_p
- xt::xtensor dot_p(const xt::xtensor &x_u, const xt::xtensor &x_p) const;
+ xt::xtensor dot(const xt::xtensor &x) const;
// assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
// WARNING: ignores any off-diagonal terms
void assemble(const xt::xtensor &elemmat);
- // set matrix components from externally assembled object
- void set (const xt::xtensor &A ) const; // diagonal [ndof]
- void set_uu(const xt::xtensor &A_uu) const; // diagonal [nnu]
- void set_pp(const xt::xtensor &A_pp) const; // diagonal [nnp]
-
- // solve
- // x = A \ b
- xt::xtensor solve (const xt::xtensor &b );
- // x = assembly{ A_uu \ ( b_u - A_up * x_p ) ; x_p }
- xt::xtensor solve (const xt::xtensor &b , const xt::xtensor &x_p);
- // x_u = A_uu \ ( b_u - A_up * x_p )
- xt::xtensor solve_u(const xt::xtensor &b_u, const xt::xtensor &x_p);
-
- // return (sub-)matrix as diagonal matrix (column)
- // assembly{ A_uu ; A_pp }
- xt::xtensor asDiagonal () const;
- // A_uu
- xt::xtensor asDiagonal_uu() const;
- // A_pp
- xt::xtensor