Symbolic MPM module (Concept & Formulation)

The symbolic MPM module allows to define problems using syntax describing their weak form. The weak form is represented using

• Variables, defining unknown fields and test fields appearing in the weak form. Variables define physical meaning of the field, its interpolation and degrees of freedom.

• Terms, representing mathematical expressions appearing in the weak form. Terms are evaluated for specific variables. Beside the library of pre-defined terms, an arbitrary term can be entered directly in the input file as a symbolic expression (SymbolicTerm, see below).

• Integrals, representing integrals of terms over specific domains. Integral domains are defined using sets.

• Materials, providing the constitutive response. A material receives the generalized state of an integration point once per equilibrium iteration and the terms then query the cached response (see Material interface: state push and declared state layout).

The characteristic discretized equations are assembled from individual integral contributions by looping over individual elements of each integral domain. Each term can evaluate itself for known variables at given time and also evaluate its linearization with respect to unknown variables. How these contributions are assembled is defined by the engineering model, representing the problem type to be solved.

We will illustrate the concept of MPM Symbolic module on an example of formulation of the elasticity problem able to represent the incompressible behavior. This formulation can be written using mixed approach involving pressure and displacement fields

\[\nabla p + 2\mu\nabla\cdot\rm{dev}[\nabla^s u]+f = 0\;\rm{in}\ \Omega\]

\[\frac{1}{K}p-\nabla\cdot u=0\;\rm{in}\ \Omega\]

\[u=\bar u\;\rm{on}\ \Gamma_u\]

\[\sigma n = \bar t\;\rm{on}\ \Gamma_t\]

The weak form of the above equations can be written as

\[\int_\Omega\overbrace{2\mu\nabla^sw:\rm{dev}[\nabla^su]}^{T_1}\ d\Omega-\int_\Omega\overbrace{ \nabla\cdot w\ p}^{T_2}\ d\Omega=\int_{\Gamma_t}\overbrace{w\cdot \bar t}^{T_3}\ d\Gamma\]

\[-\int_\Omega\underbrace{q\ \nabla\cdot u}_{T_4 = T_2^T}\ d\Omega+\int_\Omega\underbrace{\frac{1}{K}q\ p}_{T_5}\ d\Omega = 0\]

where

• \(u, p\) are unknown displacement and pressure fields (Variables),

• \(w, q\) are the corresponding test fields (Variables).

In order to obtain stable results, the approximations should satisfy Babuška-Brezzi condition or some form of stabilization would be required. To satisfy the B-B condition, the approximation order of displacement should be higher than that of hydrostatic pressure. In the following we will consider 2D case.

Variables

In the input deck, the Variables (\(u,p,w,q\)) can be set up using following syntax

Variable name "u"  interpolation "feilin"   type 1 quantity "Displacements" size 2 dofs 2 1 2
Variable name "w"  interpolation "feilin"   type 1 quantity "Displacements" size 2 dofs 2 1 2
Variable name "p"  interpolation "feiconst" type 0 quantity "Pressure"      size 1 dofs 1 11
Variable name "dp" interpolation "feiconst" type 0 quantity "Pressure"      size 1 dofs 1 11

where interpolation determines the interpolation used for specific field. Here we use linear approximation (interpolation “feilin”) for displacement field and related test field and constant approximation (interpolation “feiconst”) for pressure (so called Q1P0 element). The \(u, w\) fields are vector fields (type 1) with physical meaning of displacement (quantity “Displacements”) and two degrees of freedom (size 2 dofs 2 1 2). And \(p, dp\) are scalar fields (type 0) with physical meaning of pressure (quantity “Pressure” size 1 dofs 1 11). The test field of the pressure, denoted \(q\) in the weak form above, is named dp here, following the shipped input decks.

The quantity attribute names the physical field the variable represents. It is a FieldType, written as the enumerator name with the FT_ prefix dropped — "Displacements", "Pressure", "Pressure2", "Temperature", "Concentration1", "Concentration2". This is what the assembly matches against the state layout declared by the material (Material interface: state push and declared state layout), so it has to say what the field is; dofs says only on which degrees of freedom that field is numbered. The two declarations are independent and are cross-checked at initialization: a deck storing, say, temperature on D_w is legitimate and only produces a warning. Numeric values of quantity are still accepted but should not be used — they are an implementation detail, and the meaning of a given number has changed in the past.

A variable that is a test (weighting) function should say so, with dualto naming the unknown it weights

Variable name "dp" interpolation "feiconst" type 0 quantity "Pressure" size 1 dofs 1 11 dualto "p"

Both records are otherwise indistinguishable — they may even carry different interpolations, which is what a non-symmetric (Petrov-Galerkin) formulation looks like — and without it the assembly cannot tell which field to read nodal unknowns from when pushing the state to the material. It is required only where a deck names a test function as some term’s variable, which happens for pure source terms (those with a zero linearization), since they have no unknown to depend on.

Pre-defined terms

The weak form above consists of several terms to be evaluated

• \(T_1\): This is represented by BTSigmaTerm, evaluated for \(w\) test field and \(u\) as unknown field, under plain strain assumptions (mmode 7) and considering deviatoric part only (lhsmatmode 29), with the stress read back for the residual (rhsmatmode 30)

  BTSigmaTerm 1 variable "u"  testvariable "w" mmode 7 lhsmatmode 29 rhsmatmode 30
  

• \(T_2\): represented by BTamNTerm, evaluated for \(w\) test field and \(p\) as unknown field, under plain strain assumptions (mmode 7) and with scalar parameter equal to 1.0 (atype 28)

  BTamNTerm 2 variable "p" testvariable "w" mmode 7 atype 28
  

• Similarly, we set up remaining three terms \(T_3, T_4=T_2^T\) and \(T_5\). The uvmt 1 attribute selects the value mode in which the unknown field is read (total values instead of the default rates)

  NTamTBTerm 3 variable "u" testvariable "dp" mmode 7 atype 28 uvmt 1
  NTcN 4 variable "p" testvariable "dp" mmode 7 ctype 27 uvmt 1
  NTfTerm 5 variable "u" testvariable "w" mmode 6 flux 2 0. 6.25
  

The mmode attribute takes a MaterialMode, which may be given by name (mmode "_PlaneStrain") instead of by number. The library of available terms is documented in the Term documentation chapter. The terms are integrated over specific domains (\(\Omega,\ \Gamma_t\)), defined using corresponding sets (defined bellow in the example):

Integral 1 domain 1 set 1 term 1
Integral 2 domain 1 set 1 term 2 factor -1.0
Integral 3 domain 1 set 1 term 3 factor -1.0
Integral 4 domain 1 set 1 term 4
Integral 5 domain 1 set 2 term 5

User-defined terms: SymbolicTerm

A term does not have to exist as a C++ class. SymbolicTerm takes the mathematical expressions of the term itself, written in a small matrix language, and compiles them into bytecode at initialization; the assembly loops then execute the compiled form:

SymbolicTerm <id> variable "<name>" testvariable "<name>" mmode <mode>
             lexpression "<lhs expression>" rexpression "<rhs expression>"

lexpression defines the linearization (the tangent/stiffness contribution, a matrix) and rexpression the residual/internal force contribution (a column matrix, converted to a vector). The operators of the weak form are available as functors: Grad_s, Grad, Div and N build the interpolation operator matrices of a variable at the integration point, MDer, MVec and MProp read the material response, and ru/rv read the nodal unknowns/rates of a field on the cell.

The five terms of the Cook membrane formulation above read, in symbolic form (compare with the pre-defined terms of the previous section)

SymbolicTerm 1 variable "u" testvariable "w"  mmode 7 lexpression "Grad_s(w,gp).T*MDer(gp,ts, MatResponseMode::DeviatoricStiffness)*Grad_s(u,gp)" rexpression "Grad_s(w,gp).T*Sig_dev(u,gp,ts)"
SymbolicTerm 2 variable "p" testvariable "w"  mmode 7 lexpression "Div(w,gp).T*N(p,gp)" rexpression "Div(w,gp).T*N(p,gp)*ru(p, cell, ts)"
SymbolicTerm 3 variable "u" testvariable "dp" mmode 7 lexpression "N(dp,gp).T*Div(u,gp)" rexpression "N(dp,gp).T*Div(u,gp)*ru(u, cell, ts)"
SymbolicTerm 4 variable "p" testvariable "dp" mmode 7 lexpression "N(dp,gp).T*2.4e-7*N(p,gp)" rexpression "N(dp,gp).T*2.4e-7*N(p,gp)*ru(p, cell, ts)"
NTfTerm     5 variable "u" testvariable "w"  mmode 6 flux 2 0. 6.25

with the deck otherwise unchanged; the complete input file is tests/regression/mpm/mpms_cook2_u2p1.in. Response modes are best given by name (MatResponseMode::DeviatoricStiffness) rather than by their integer value. The full description of the expression language, of the available functors and of the evaluation mechanics is in the SymbolicTerm expression language section of the OOFEM Input manual.

Material interface: state push and declared state layout

The material interface used by the module is split into two halves. Once per equilibrium iteration, before any term is evaluated, the solver hands each integration point’s material its generalized state — strains, field gradients, field values — and the material performs all of its constitutive work there, caching the results:

virtual void updateTempState(const FloatArray &stateVector, GaussPoint *gp, TimeStep *tStep);

The terms then only read that cache, selecting what they need by MatResponseMode:

virtual void   giveCharacteristicVector(FloatArray &answer,  MatResponseMode, GaussPoint*, TimeStep*) const;
virtual void   giveCharacteristicMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint*, TimeStep*) const;
virtual double giveCharacteristicValue (              MatResponseMode, GaussPoint*, TimeStep*) const;

This matters for the formulation because the tangent and the residual are assembled in separate sweeps over the integration points: with the state pushed once per iteration, the response a term reads no longer depends on the order in which the sweeps happen to run, and a term never has to assemble and pass the state itself.

What the state vector contains, and in which order, is declared by the material rather than fixed by convention:

virtual StateVariableLayout giveStateVariableIDs(MaterialMode mmode) const;

Each entry of the returned layout is a {FieldType field, StateOperator op, ValueModeType mode} triple, in packing order: which primary field, what is taken of it, and of which value of the field; the value mode may be omitted (see below). For each entry the assembly finds the unknown variable whose quantity is field, reads its nodal unknowns of the given value mode on the cell, applies the operator at the integration point and appends the resulting block. Four operators are available:

StateOperator	Meaning	Block size
SO_Value	the field itself, \(Nu\)	Variable::size
SO_Gradient	gradient of a scalar field, \(\nabla u\)	number of spatial dimensions
SO_SymmetricGradient	symmetric gradient of a vector field, \(\nabla^s u\) (strain, Voigt form)	material mode dependent (6 in _3dMat and the UP modes, 4 in _PlaneStrain, 3 in _PlaneStress, 1 in _1dMat)
SO_Divergence	divergence of a vector field, \(\nabla\cdot u\)	1

So a poromechanical material consuming strain, pressure gradient and pressure declares

StateVariableLayout giveStateVariableIDs(MaterialMode mmode) const override {
    return { { FT_Displacements, SO_SymmetricGradient },
             { FT_Pressure,      SO_Gradient },
             { FT_Pressure,      SO_Value } };
}

and receives \([\,\varepsilon,\ \nabla p,\ p\,]\) packed in that order.

The value mode selects the nodal unknowns the operator is applied to. It defaults to VM_TotalIntrinsic, so an entry written as { field, op } reads the values at the intrinsic time of the step, as before. With the generalized trapezoidal rule of the nonstationary problem the modes are

ValueModeType	Nodal values
VM_TotalIntrinsic (default)	\(\alpha x_{n+1} + (1-\alpha) x_n\)
VM_Total	\(x_{n+1}\)
VM_Velocity	\((x_{n+1} - x_n)/\Delta t\)
VM_Incremental	\(x_{n+1} - x_n\)

The operators are linear, so a rate-dependent material obtains, for example, the strain rate with { FT_Displacements, SO_SymmetricGradient, VM_Velocity } (see the Kelvin–Voigt option of the test material statedeptangentmat and tests/regression/mpm/mpms_staterate.in).

The field axis is matched against the deck’s quantity keyword rather than against dof ids, so a material never obliges a deck to renumber its dofs. Any FieldType of src/core/field.h may be named; these formulations use FT_Displacements, FT_Pressure, FT_Pressure2, FT_Temperature, FT_Concentration1 and FT_Concentration2. Naming the same field twice with different operators is normal, and is how a material asks for both a quantity and its gradient. Only the four operators above are available; anything else is rejected when the state is pushed, with unsupported state operator.

Three further properties are worth noting. The layout is queried per integration point with that point’s material mode, so a material may advertise one layout in _3dUP and decline in modes it does not implement. A material returning an empty layout — the default — does not participate in this protocol and keeps using its physics-specific entry points. And a cell that cannot supply every declared field (a thermo-mechanical material driven for the thermal sub-problem alone, say) is skipped rather than given a partial state, which would be positionally ambiguous; a later query that does need the missing value reports it as unset rather than returning stale data.

The same interface is available to materials implemented in python, where giveStateVariableIDs(mmode) returns a sequence of (FieldType, StateOperator) pairs or (FieldType, StateOperator, ValueModeType) triples, the enumerations being exported by the oofem module:

import oofem

class MyMaterial:
    def giveStateVariableIDs(self, mmode):
        return [(oofem.FieldType.FT_Displacements,
                 oofem.StateOperator.SO_SymmetricGradient),
                (oofem.FieldType.FT_Pressure, oofem.StateOperator.SO_Value),
                (oofem.FieldType.FT_Pressure2, oofem.StateOperator.SO_Value),
                # strain rate
                (oofem.FieldType.FT_Displacements,
                 oofem.StateOperator.SO_SymmetricGradient,
                 oofem.ValueModeType.VM_Velocity)]

    def updateTempState(self, stateVector, gp, tStep, stateDict, tempStateDict):
        eps = stateVector[0:6]
        pw, pa = stateVector[6], stateVector[7]
        epsRate = stateVector[8:14]
        ...   # all constitutive work here, results into tempStateDict

A python material that defines updateTempState but no giveStateVariableIDs advertises an empty layout and is therefore never pushed to; one that defines neither is reported at initialization, naming what to change.

For a field to be available to be pushed at all, some term on the cell must name it as its unknown variable; test functions are excluded through dualto. The design rationale for the (field, operator) split and for the push/pull contract is described in doc/unified_material_interface.md.

Example: Cook membrane

The Cook’s membrane is a standard benchmark problem. It consists of a tapered plate clamped on one of its sides with a transversal distributed load applied to the opposite side. The plate is in plain strain and its dimensions, as well as the material parameters and boundary conditions, are shown in figure bellow

Fig. 1 Cook membrane geometry and boundary conditions

The complete OOFEM input deck for mesh consisting of 2x2 is shown below

cook2_u1p0_2.out
Demo of symbolic mpm problem; Cook membrane benchmark
#
mpmsymbolicstationaryproblem nsteps 1 nvariables 4 nterms 5 nintegrals 5 lhsterms 4 1 2 3 4 rhsterms 1 5 nmodules 1 profileopt 1 maxiter 10
errorcheck
Variable name "u" interpolation "feilin" type 1 quantity "Displacements" size 2 dofs 2 1 2 # displacement
Variable name "w" interpolation "feilin" type 1 quantity "Displacements" size 2 dofs 2 1 2 # test function
Variable name "p" interpolation "feiconst"  type 0 quantity "Pressure" size 1 dofs 1 11 # pressure
Variable name "dp" interpolation "feiconst" type 0 quantity "Pressure" size 1 dofs 1 11 # test function
BTSigmaTerm 1 variable "u"  testvariable "w" mmode 7 lhsmatmode 29 rhsmatmode 30
BTamNTerm 2 variable "p" testvariable "w" mmode 7 atype 28
NTamTBTerm 3 variable "u" testvariable "dp" mmode 7 atype 28 uvmt 1
NTcN 4 variable "p" testvariable "dp" mmode 7 ctype 27 uvmt 1
NTfTerm 5 variable "u" testvariable "w" mmode 6 flux 2 0. 6.25
Integral 1 domain 1 set 1 term 1
Integral 2 domain 1 set 1 term 2 factor -1.0
Integral 3 domain 1 set 1 term 3 factor -1.0
Integral 4 domain 1 set 1 term 4
Integral 5 domain 1 set 2 term 5
domain HeatTransfer
outputmanager tstep_all dofman_all element_all
ndofman 9 nelem 6 nbc 1 ncrosssect 1 nic 0 nltf 2 nmat 1 nset 3
Node 1 coords 3 0.0 0.0 0.0
Node 2 coords 3 24.0 22.0 0.0
Node 3 coords 3 48.0 44.0 0.0
Node 4 coords 3 0.0 22.0 0.0
Node 5 coords 3 24.0 37.0 0.0
Node 6 coords 3 48.0 52.0 0.0
Node 7 coords 3 0.0 44.0 0.0
Node 8 coords 3 24.0 52.0 0.0
Node 9 coords 3 48.0 60.0 0.0
q1 1 nodes 4 1 2 5 4 mat 1 crosssect 1
q1 2 nodes 4 2 3 6 5 mat 1 crosssect 1
q1 3 nodes 4 4 5 8 7 mat 1 crosssect 1
q1 4 nodes 4 5 6 9 8 mat 1 crosssect 1
l1 5 nodes 2 3 6 mat 1 crosssect 1
l1 6 nodes 2 6 9 mat 1 crosssect 1
simplecs 1 thick 5.0
isole 1 d 1 e 250 n 0.49999 talpha 1.
# clamped-displacement
boundarycondition 1 loadtimefunction 1 set 3 values 2 0 0   dofs 2 1 2
constantfunction 1 f(t) 1
PiecewiseLinFunction 2 nPoints 4 t 4 -10. 0. 1. 5. f(t) 4 0. 0. 1.0 1.0
set 1 elementranges  {(1 4)}
set 2 elementranges  {(5 6)}
set 3 elementedges 4 1 4 3 4

The complete input deck can be found in tests/regression/mpm/cook2_u1p0_2.in file; the variant driven by the lightweight test engineering model is tests/regression/mpm/cook2_u1p0.in.

Switching to high-order interpolation

To switch to high order element Q2P1 with quadratic interpolation of displacement and linear interpolation of pressure we need just to use different interpolations for variables. The updated variable input records are following

Variable name "u" interpolation "feiquad" type 1 quantity "Displacements" size 2 dofs 2 1 2 # displacement
Variable name "w" interpolation "feiquad" type 1 quantity "Displacements" size 2 dofs 2 1 2 # test function
Variable name "p" interpolation "feilin"  type 0 quantity "Pressure" size 1 dofs 1 11 # pressure
Variable name "dp" interpolation "feilin" type 0 quantity "Pressure" size 1 dofs 1 11 # test function

The complete input deck can be found in tests/regression/mpm/cook2_u2p1.in file, and its SymbolicTerm counterpart in tests/regression/mpm/mpms_cook2_u2p1.in.

To illustrate the convergence, sequence of uniform meshes of the plate is considered, starting from a mesh consisting of a single elements and proceeding by uniform refinement.

The figure below shows the vertical displacement of the plane tip plotted against the number of element segments along each side. The solution is compared to reference solution [1].

References

[1]

Ignacio Romero, Manfred Bischoff, Incompatible Bubbles: A non-conforming finite element formulation for linear elasticity, Computer Methods in Applied Mechanics and Engineering, Volume 196, Issues 9–12, 2007, Pages 1662-1672, ISSN 0045-7825, https://doi.org/10.1016/j.cma.2006.09.010.