[GeoMechanicsApplication] Investigate extracting the equation of motion functionality from the elements

[rfaasse]

As a developer, I want to investigate how to separate the equation of motion functionality from the elements, such that the elements become more maintainable and configurable.

In discussions we discussed the following proposal for the equation of motion policy structure:

In this investigation, the feasibility and (if possible) concrete steps to implement this structure are documented

Acceptance criteria

• A follow-up issue is created and refined for the implementation or further investigations
• Investigation is documented

[rfaasse]

Outcome of the investigation

This investigation has the goal of finding out how to implement the different configurations for solving the equation of motion (steady state, transient and dynamic) into a policy structure for the element.

In summary, this investigation concludes that these different scenarios are currently not handled by the elements themselves, but by the different flavors of time integration schemes. Therefore, at this moment, there is no need for extraction of these policies for the elements, but we can keep this functionality in the schemes. However, although this is not part of the element redesign, there are still missing configurations in the schemes (e.g. Backward Euler does not have dynamic/transient versions).

Background

The idea of the equation of motion policy is that it should handle the top line of the following matrix equations, i.e. the displacements. There are three possible scenarios (as also described in the policy structure in the description of this issue): transient, steady state and dynamic, where different terms are taken into account:

• For steady state, only the stiffness matrix (K) is taken into account in the equations, while the other matrices (M and D) are not considered
• For transient, additional to the stiffness matrix, also the damping matrix (D) is considered
• For dynamic, also the inertia term, or the mass matrix (M) is taken into account.

It is seen that for whichever policy is used, the stiffness matrix needs to be added. Therefore the element will always add that contribution to the LHS/RHS.

The conditional matrices (M/D) are where this policy would start adding value to easily configure which displacement terms are taken into account. The elements have functions (which are overrides from the base Element class from the Kratos Core), to calculate the mass and damping matrices (i.e. CalculateMassMatrix and CalculateDampingMatrix). However, the element itself never calls this function (except for the GeoStructural elements, which are out of scope for this investigation on the Geo/UPw/Pw elements). The only place where these functions are called are the time integration schemes. For example the NewmarkDynamicUPwScheme has the following function to add terms to the LHS:

void CalculateLHSContribution(Element&                       rCurrentElement,
                              LocalSystemMatrixType&         LHS_Contribution,
                              Element::EquationIdVectorType& EquationId,
                              const ProcessInfo&             CurrentProcessInfo) override
{
    KRATOS_TRY

    int thread = OpenMPUtils::ThisThread();

    rCurrentElement.CalculateLeftHandSide(LHS_Contribution, CurrentProcessInfo);
    rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
    rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
    this->AddDynamicsToLHS(LHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], CurrentProcessInfo);
    rCurrentElement.EquationIdVector(EquationId, CurrentProcessInfo);

    KRATOS_CATCH("")
}

So only the dynamic scheme adds the mass and damping matrices to the LHS, while the quasistatic_with_damping scheme only adds damping and omits the mass matrix. Therefore it is seen the configuration for the different ways of solving the equation of motion (either dynamic, transient or steady state) are taken care of by the time integration schemes. Therefore, it does not need to be part of the element re-design to add an equation of motion policy that adds mass/damping conditionally.

One note is that the schemes are not complete (e.g. the backward_euler schemes have no dynamic/transient versions). This means that although the equation of motion policy does not need to be part of the element redesign, we still need to take steps to create the dynamic/transient/steady-state flavors for all available time integration schemes.

[WPK4FEM]

I completely support the above analysis by Richard.

[avdg81]

Very clear and well-written analysis. Thank you Richard!

[mnabideltares]

It is very nicely written, Thank you Richard.

[rfaasse]

Closing this issue after review from team members. The acceptance criteria are met, as the investigation is documented and no follow-up items for creating the policy within the elements is needed (the investigation of the missing scheme configuration can be found in a separate issue, which was already created before.