Type traits for mechanical problems

[dbeurle]

Fixed size numerical structures (vectors and matrices) are critical to extracting the maximum performance out of the code by providing as much information as possible to the optimiser. Depending on each simulation type, the matrices representing the material tangent, stress and deformation tensors are known. Some generic classes, for example the internal variables class, already take template parameters for the size of the rank two and rank four tensor. The dimension information (one, two or three) and the degrees of freedom can be encoded into these traits.

To deal with the quirks of each simulation type (axisymmetric, plate, shell, solid, beam, plane strain and plane stress) we should introduce a traits class that specifies this information for each simulation type. This generates a little additional work but should encapsulate the size information.

Firstly, we need to introduce a strongly-typed enum to provide a compile-time logic.

#include <type_traits>

enum class type : uint8_t {plane_stress, plane_strain, solid, plate, shell, axisymmetric};

template <type T, bool>
struct traits;

template <type::solid, bool is_symmetric = true>
struct traits
{
    static auto constexpr value_type = type::solid;
    static auto constexpr size = 3;

    using rank_two_tensor = matrix3;
    using rank_four_tensor = std::conditional<is_symmetric, matrix6, matrix9>::type;

    using internal_variables_type = internal_variables<rank_two_tensor, rank_four_tensor>;
    using constitutive_type = constitutive_model<type::solid, internal_variables_type>;
};

This should be able to handle the unsymmetrical 9x9 constitutive models for three-dimensional theory if required through the template parameter. In addition, if analysis specific post-processing steps are required, further template specialisations can be used for the parts which are not common.

A downside of this approach is the additional boiler plate and the fact that most of the implemented code will still be dependent on the type with what code can be reused, will be.

[dbeurle]

Constitutive model

The design for the constitutive model needs to handle the stress update and tangent matrix operator at the same time. The explicit simulations do not require a material tangent but only the stress update.

A constitutive model dictates:

• (Primary) symmetric or unsymmetric material tangent
• (Secondary) form of the gradient operator (B)

The dependency starts with the choice of the constitutive model. This should then have an is_symmetric compile time flag. This in turn dictates the internal variable matrix size and then allows the correct mesh to be selected with the corresponding element assembly kernels.

[dbeurle]

An outline of all the possible cases this design would need to account for are listed here.

Linear theory

Displacement formulation

The assembly and solution of the linear system K*d = f_ext for each common mechanical theory with a finite element implementation is characterised by the following:

Module	Coordinates	Unknowns
Beam	x, y, z	u1, u2, u3, r1, r2, r3 (displacement and rotation)
Plane stress	x, y	u1, u2 (displacement)
Plane strain	x, y	u1, u2 (displacement)
Axisymmetric	r, t (cylindrical coordinates)	u1, u2 (cylindrical coordinates)
Shell	x, y, z	u1, u2 ,u3, m1, m2, m3
Solid	x, y, z	u1, u2, u3

and each theory uses a special form of the B and D (gradient and material tangent operators respectively) and are tabulated below

Module	Gradient operator	Notes
Plane stress	B (3x3)	Identical to plane strain theory
Plane strain	B (3x3)	Identical to plane stress theory
Solid	B (6x3)	Implemented for symmetric operator

Linear eigenvalue buckling

For buckling problems the above displacement formulation is required to be solved along with the assembly pre-stress matrix or geometric stiffness matrix for the current configuration. An eigenvalue problem results for the buckling. The requirement is that the linear material tangent is build alongside the geometric stiffness.

Natural frequency analysis

Natural frequency problems require the solution of K*x = lambda*M*x. This requires that the finite element mesh (fem_mesh) provides methods for computing the element mass and stiffness matrices.

Non-linear theory

Updated and total Lagrange formulations. LATIN for non-incremental method.

[dbeurle]

The enum characterisations could be:

namespace mechanical
{
enum class type {
                 /// C1 beam theory with Hermite interpolation
                 euler_beam,
                 /// C0 beam theory
                 beam,
                 shell,
                 axisymmetric,
                 plane_stress,
                 plane_strain,
                 solid};

enum class material {linear, nonlinear};

enum class deformation {small, finite};

enum class method {perturbation, incremental, latin};

enum class lagrangian {total, updated};
}

where the impossible combinations could be disallowed at compile time.