File Element.h#

Convenience methods for integration point data.

Copyright: Copyright 2017. Tom de Geus. All rights reserved.
License: This project is released under the GNU Public License (GPLv3).

namespace GooseFEM

Toolbox to perform finite element computations.

namespace Element

Element quadrature and interpolation.

Functions

inline array_type::tensor<double, 3> asElementVector(const array_type::tensor<size_t, 2> &conn, const array_type::tensor<double, 2> &nodevec)#

Convert nodal vector with (“nodevec”, shape:[nnode, ndim]) to nodal vector stored per element (“elemvec”, shape: [nelem, nne, ndim]).

Parameters:

conn – Connectivity.
nodevec – “nodevec”.

Returns:

“elemvec”.

inline array_type::tensor<double, 2> assembleNodeVector(const array_type::tensor<size_t, 2> &conn, const array_type::tensor<double, 3> &elemvec)#

Assemble nodal vector stored per element (“elemvec”, shape [nelem, nne, ndim]) to nodal vector (“nodevec”, shape [nnode, ndim]).

Parameters:

conn – Connectivity.
elemvec – “elemvec”.

Returns:

“nodevec”.

template<class E> inline bool isSequential(const E &dofs)#

Check that DOFs leave no holes.

Parameters:: dofs – DOFs (“nodevec”)
Returns:: true if there are no holds.

bool isDiagonal(const array_type::tensor<double, 3> &elemmat)#

Check that all of the matrices stored per elemmat (shape: [nelem, nne * ndim, nne * ndim]) are diagonal.

Parameters:: elemmat – Element-vectors (“elemmat”)
Returns:: true if all element matrices are diagonal.

template<class D> class QuadratureBase

#include <GooseFEM/Element.h>

CRTP base class for quadrature.

Subclassed by GooseFEM::Element::QuadratureBaseCartesian< Quadrature >, GooseFEM::Element::QuadratureBaseCartesian< QuadratureAxisymmetric >, GooseFEM::Element::QuadratureBaseCartesian< QuadraturePlanar >, GooseFEM::Element::QuadratureBaseCartesian< D >

Public Types

using derived_type = D: Underlying type.

Public Functions

inline auto nelem() const

Number of elements.

Returns:: Scalar.

inline auto nne() const

Number of nodes per element.

Returns:: Scalar.

inline auto ndim() const

Number of dimensions for node vectors.

Returns:: Scalar.

inline auto tdim() const

Number of dimensions for integration point tensors.

Returns:: Scalar.

inline auto nip() const

Number of integration points.

Returns:: Scalar.

template<class T, class R> inline void asTensor(const T &arg, R &ret) const

Convert “qscalar” to “qtensor” of certain rank.

Fully allocated output passed as reference, use AsTensor to allocate and return data.

Parameters:

arg – A “qscalar”.
ret – A “qtensor”.

template<size_t rank, class T> inline auto AsTensor(const T &arg) const

Convert “qscalar” to “qtensor” of certain rank.

Parameters:: arg – A “qscalar”.
Returns:: “qtensor”.

template<class T> inline auto AsTensor(size_t rank, const T &arg) const

Convert “qscalar” to “qtensor” of certain rank.

Parameters:

rank – Tensor rank.
arg – A “qscalar”.

Returns:

“qtensor”.

inline auto shape_elemvec() const -> std::array<size_t, 3>

Get the shape of an “elemvec”.

Returns:: [nelem, nne, ndim].

inline auto shape_elemvec(size_t arg) const -> std::array<size_t, 3>

Get the shape of an “elemvec”.

Parameters:: arg – The vector dimension.
Returns:: [nelem, nne, tdim].

inline auto shape_elemmat() const -> std::array<size_t, 3>

Get the shape of an “elemmat”.

Returns:: [nelem, nne * ndim, nne * ndim].

template<size_t rank = 0> inline auto shape_qtensor() const -> std::array<size_t, 2 + rank>

Get the shape of a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Template Parameters:: rank – Rank of the tensor. Output is fixed-size: std::array<size_t, rank>.
Returns:: [nelem, nip, tdim, …].

inline auto shape_qtensor(size_t rank) const -> std::vector<size_t>

Get the shape of a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Parameters:: rank – The tensor rank.
Returns:: [nelem, nip, tdim, …].

template<size_t trank> inline auto shape_qtensor(size_t rank, size_t arg) const -> std::array<size_t, 2 + trank>

Get the shape of a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Template Parameters:

rank – Rank of the tensor. Output is fixed-size: std::array<size_t, rank>.

Parameters:

rank – The tensor rank. Effectively useless: to distinguish from the dynamic-sized.
arg – The tensor dimension.

Returns:

[nelem, nip, tdim, …].

inline auto shape_qtensor(size_t rank, size_t arg) const -> std::vector<size_t>

Get the shape of a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Parameters:

rank – The tensor rank.
arg – The tensor dimension.

Returns:

[nelem, nip, tdim, …].

inline auto shape_qscalar() const -> std::array<size_t, 2>

Get the shape of a “qscalar” (a “qtensor” of rank 0)

Returns:: [nelem, nip].

inline auto shape_qvector() const -> std::array<size_t, 3>

Get the shape of a “qvector” (a “qtensor” of rank 1)

Returns:: [nelem, nip, tdim].

inline auto shape_qvector(size_t arg) const -> std::array<size_t, 3>

Get the shape of a “qvector” (a “qtensor” of rank 1)

Parameters:: arg – Tensor dimension.
Returns:: [nelem, nip, tdim].

template<class R> inline auto allocate_elemvec() const

Get an allocated array_type::tensor to store a “elemvec”.

Note: the container is not (zero-)initialised.

Template Parameters:: R – value-type of the array, e.g. double.
Returns:: xt::xarray container of the correct shape.

template<class R> inline auto allocate_elemvec(R val) const

Get an allocated and initialised xt::xarray to store a “elemvec”.

Template Parameters:: R – value-type of the array, e.g. double.
Parameters:: val – The value to which to initialise all items.
Returns:: array_type::tensor container of the correct shape.

template<class R> inline auto allocate_elemmat() const

Get an allocated array_type::tensor to store a “elemmat”.

Note: the container is not (zero-)initialised.

Template Parameters:: R – value-type of the array, e.g. double.
Returns:: xt::xarray container of the correct shape.

template<class R> inline auto allocate_elemmat(R val) const

Get an allocated and initialised xt::xarray to store a “elemmat”.

Template Parameters:: R – value-type of the array, e.g. double.
Parameters:: val – The value to which to initialise all items.
Returns:: array_type::tensor container of the correct shape.

template<size_t rank = 0, class R> inline auto allocate_qtensor() const

Get an allocated array_type::tensor to store a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Default: rank = 0, a.k.a. scalar. Note: the container is not (zero-)initialised.

Template Parameters:: R – value-type of the array, e.g. double.
Returns:: [nelem, nip].

template<size_t rank = 0, class R> inline auto allocate_qtensor(R val) const

Get an allocated and initialised array_type::tensor to store a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Default: rank = 0, a.k.a. scalar.

Template Parameters:: R – value-type of the array, e.g. double.
Parameters:: val – The value to which to initialise all items.
Returns:: array_type::tensor container of the correct shape (and rank).

template<class R> inline auto allocate_qtensor(size_t rank) const

Get an allocated xt::xarray to store a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Note: the container is not (zero-)initialised.

Template Parameters:: R – value-type of the array, e.g. double.
Parameters:: rank – The tensor rank.
Returns:: xt::xarray container of the correct shape.

template<class R> inline auto allocate_qtensor(size_t rank, R val) const

Get an allocated and initialised xt::xarray to store a “qtensor” of a certain rank (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).

Template Parameters:

R – value-type of the array, e.g. double.

Parameters:

rank – The tensor rank.
val – The value to which to initialise all items.

Returns:

array_type::tensor container of the correct shape (and rank).

template<class R> inline auto allocate_qscalar() const

Get an allocated array_type::tensor to store a “qscalar” (a “qtensor” of rank 0).

Note: the container is not (zero-)initialised.

Template Parameters:: R – value-type of the array, e.g. double.
Returns:: xt::xarray container of the correct shape.

template<class R> inline auto allocate_qscalar(R val) const

Get an allocated and initialised xt::xarray to store a “qscalar” (a “qtensor” of rank 0).

Template Parameters:: R – value-type of the array, e.g. double.
Parameters:: val – The value to which to initialise all items.
Returns:: array_type::tensor container of the correct shape (and rank).

Private Functions

inline auto derived_cast() -> derived_type&#

inline auto derived_cast() const -> const derived_type&#

template<class D> class QuadratureBaseCartesian : public GooseFEM::Element::QuadratureBase<D>

#include <GooseFEM/Element.h>

CRTP base class for interpolation and quadrature for a generic element in Cartesian coordinates.

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]

Public Types

using derived_type = D: Underlying type.

Public Functions

template<class T> inline void update_x(const T &x)

Update the nodal positions.

This recomputes:

the shape functions,
the shape function gradients (in local and global) coordinates,
the integration points volumes. Under the small deformations assumption this function should not be called.

Parameters:: x – nodal coordinates (elemvec). Shape should match the earlier definition.

inline auto GradN() const -> const array_type::tensor<double, 4>&

Shape function gradients (in global coordinates).

Returns:: gradN stored per element, per integration point [nelem, nip, nne, ndim].

inline auto dV() const -> const array_type::tensor<double, 2>&

Integration volume.

Returns:: volume stored per element, per integration point [nelem, nip].

template<class T> inline auto InterpQuad_vector(const T &elemvec) const -> array_type::tensor<double, 3>

Interpolate element vector and evaluate at each quadrature point.

\( \vec{u}(\vec{x}_q) = N_i^e(\vec{x}) \vec{u}_i^e \)

Parameters:: elemvec – nodal vector stored per element [nelem, nne, ndim].
Returns:: qvector [nelem, nip, ndim].

template<class T, class R> inline void interpQuad_vector(const T &elemvec, R &qvector) const

Same as InterpQuad_vector(), but writing to preallocated return.

Parameters:

elemvec – nodal vector stored per element [nelem, nne, ndim].
qvector – [nelem, nip, ndim].

template<class T> inline auto GradN_vector(const T &elemvec) const -> array_type::tensor<double, 4>

Element-by-element: dyadic product of the shape function gradients and a nodal vector.

Typical input: nodal displacements. Typical output: quadrature point strains. Within one element::

 for e in range(nelem):
     for q in range(nip):
         for m in range(nne):
             qtensor(e, q, i, j) += dNdx(e, q, m, i) * elemvec(e, m, j)

Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x().

Parameters:: elemvec – [nelem, nne, ndim]
Returns:: qtensor [nelem, nip, tdim, tdim]

template<class T, class R> inline void gradN_vector(const T &elemvec, R &qtensor) const

Same as GradN_vector(), but writing to preallocated return.

Parameters:

elemvec – [nelem, nne, ndim]
qtensor – overwritten [nelem, nip, tdim, tdim]

template<class T> inline auto GradN_vector_T(const T &elemvec) const -> array_type::tensor<double, 4>

The transposed output of GradN_vector().

Within one element::

 for e in range(nelem):
     for q in range(nip):
         for m in range(nne):
             qtensor(e, q, j, i) += dNdx(e, q, m, i) * elemvec(e, m, j)

Parameters:: elemvec – [nelem, nne, ndim]
Returns:: qtensor [nelem, nip, tdim, tdim]

template<class T, class R> inline void gradN_vector_T(const T &elemvec, R &qtensor) const

Same as GradN_vector_T(), but writing to preallocated return.

Parameters:

elemvec – [nelem, nne, ndim]
qtensor – overwritten [nelem, nip, tdim, tdim]

template<class T> inline auto SymGradN_vector(const T &elemvec) const -> array_type::tensor<double, 4>

The symmetric output of GradN_vector().

Without one element::

 for e in range(nelem):
     for q in range(nip):
         for m in range(nne):
             qtensor(e, q, i, j) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j)
             qtensor(e, q, j, i) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j)

Parameters:: elemvec – [nelem, nne, ndim]
Returns:: qtensor [nelem, nip, tdim, tdim]

template<class T, class R> inline void symGradN_vector(const T &elemvec, R &qtensor) const

Same as SymGradN_vector(), but writing to preallocated return.

Parameters:

elemvec – [nelem, nne, ndim]
qtensor – overwritten [nelem, nip, tdim, tdim]

template<class T> inline auto Int_N_vector_dV(const T &qvector) const -> array_type::tensor<double, 3>

Element-by-element: integral of a continuous vector-field.

\( \vec{f}_i^e = \int N_i^e(\vec{x}) \vec{f}(\vec{x}) d\Omega_e \)

which is integration numerically as follows

\( \vec{f}_i^e = \sum\limits_q N_i^e(\vec{x}_q) \vec{f}(\vec{x}_q) \)

Parameters:: qvector – [nelem, nip. ndim]
Returns:: elemvec [nelem, nne. ndim]

template<class T, class R> inline void int_N_vector_dV(const T &qvector, R &elemvec) const

Same as Int_N_vector_dV(), but writing to preallocated return.

Parameters:

qvector – [nelem, nip. ndim]
elemvec – overwritten [nelem, nne. ndim]

template<class T> inline auto Int_N_scalar_NT_dV(const T &qscalar) const -> array_type::tensor<double, 3>

Element-by-element: integral of the scalar product of the shape function with a scalar.

Within one one element::

 for e in range(nelem):
     for q in range(nip):
         for m in range(nne):
             for n in range(nne):
                 elemmat(e, m * ndim + i, n * ndim + i) +=
                     N(e, q, m) * qscalar(e, q) * N(e, q, n) * dV(e, q)

with i a tensor dimension. Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x().

Parameters:: qscalar – [nelem, nip]
Returns:: elemmat [nelem, nne * ndim, nne * ndim]

template<class T, class R> inline void int_N_scalar_NT_dV(const T &qscalar, R &elemmat) const

Same as Int_N_scalar_NT_dV(), but writing to preallocated return.

Parameters:

qscalar – [nelem, nip]
elemmat – overwritten [nelem, nne * ndim, nne * ndim]

template<class T> inline auto Int_gradN_dot_tensor2_dV(const T &qtensor) const -> array_type::tensor<double, 3>

Element-by-element: integral of the dot product of the shape function gradients with a second order tensor.

Typical input: stress. Typical output: nodal force. Within one one element::

 for e in range(nelem):
     for q in range(nip):
         for m in range(nne):
             elemvec(e, m, j) += dNdx(e, q, m, i) * qtensor(e, q, i, j) * dV(e, q)

with i and j tensor dimensions. Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x().

Parameters:: qtensor – [nelem, nip, ndim, ndim]
Returns:: elemvec [nelem, nne. ndim]

template<class T, class R> inline void int_gradN_dot_tensor2_dV(const T &qtensor, R &elemvec) const

Same as Int_gradN_dot_tensor2_dV(), but writing to preallocated return.

Parameters:

qtensor – [nelem, nip, ndim, ndim]
elemvec – overwritten [nelem, nne. ndim]

template<class T> inline auto Int_gradN_dot_tensor4_dot_gradNT_dV(const T &qtensor) const -> array_type::tensor<double, 3>

Element-by-element: integral of the dot products of the shape function gradients with a fourth order tensor.

Typical input: stiffness tensor. Typical output: stiffness matrix. Within one one element::

 for e in range(nelem):
     for q in range(nip):
         for m in range(nne):
             for n in range(nne):
                 elemmat(e, m * ndim + j, n * ndim + k) +=
                     dNdx(e,q,m,i) * qtensor(e,q,i,j,k,l) * dNdx(e,q,n,l) * dV(e,q)

with i, j, k, and l tensor dimensions. Note that the functions and their gradients are precomputed upon construction, or updated when calling update_x().

Parameters:: qtensor – [nelem, nip, ndim, ndim, ndim, ndim]
Returns:: elemmat [nelem, nne * ndim, nne * ndim]

template<class T, class R> inline void int_gradN_dot_tensor4_dot_gradNT_dV(const T &qtensor, R &elemmat) const

Same as Int_gradN_dot_tensor4_dot_gradNT_dV(), but writing to preallocated return.

Parameters:

qtensor – [nelem, nip, ndim, ndim, ndim, ndim]
elemmat – overwritten [nelem, nne * ndim, nne * ndim]

Protected Functions

inline void compute_dN()#: Update the shape function gradients (called when the nodal positions are updated).

Private Functions

inline auto derived_cast() -> derived_type&#

inline auto derived_cast() const -> const derived_type&#

template<class T, class R> inline void interpQuad_vector_impl(const T &elemvec, R &qvector) const#

template<class T, class R> inline void gradN_vector_impl(const T &elemvec, R &qtensor) const#

template<class T, class R> inline void gradN_vector_T_impl(const T &elemvec, R &qtensor) const#

template<class T, class R> inline void symGradN_vector_impl(const T &elemvec, R &qtensor) const#

template<class T, class R> inline void int_N_vector_dV_impl(const T &qvector, R &elemvec) const#

template<class T, class R> inline void int_N_scalar_NT_dV_impl(const T &qscalar, R &elemmat) const#

template<class T, class R> inline void int_gradN_dot_tensor2_dV_impl(const T &qtensor, R &elemvec) const#

template<class T, class R> inline void int_gradN_dot_tensor4_dot_gradNT_dV_impl(const T &qtensor, R &elemmat) const#

inline void compute_dN_impl()#

Friends

friend class QuadratureBase< D >