`mgis.fenics`

moduleThese pages are intended to describe the `FEniCS`

interface to `MGIS`

available in the `mgis.fenics`

Python module. This module has been developed in the spirit of providing sa upport to the implementation of generalized nonlinear behaviours i.e.Â including multiple gradients and dual flux variables. It can also be used as a simple interface to standard mechanical behaviours in small or finite strain (see the plasticity demos below). We briefly describe here the general concepts underlying the module implementation.

Repository: a repository containing the demos sources files is available here

The provided documented demos have been designed to progressively illustrate the use of the interface and the versatility of the approach when implementing complex generalized behaviours, both on the `MFront`

and `FEniCS`

sides. We recommend browsing the demos in the following order:

- Stationnary non-linear heat transfer
- Stationnary non-linear heat transfer: 3D problem and performance comparisons
- Transient heat equation with phase change
- Monolithic transient thermoelasticity
- Small-strain von Mises elastoplasticity
- Finite-strain elastoplasticity within the logarithmic strain framework
- Multiphase model for fiber-reinforced materials
- Phase-field approach to brittle fracture

`mgis.fenics`

moduleThis module has been developed based on the concepts exposed in the Elasto-plastic analysis implemented using the `MFront`

code generator demo published on Numerical Tours of Computational Mechanics using FEniCS. We suggest reading first this demo as an introduction to this module implementation concepts. In particular, the module relies on the notion of `Quadrature`

function spaces to express the constitutive relation at the quadrature points level.

Compared to this initial solution, the `mgis.fenics`

module has several advantages and new features:

- The module reduces the size of boiler-plate code. In particular,
`Quadrature`

functionsâ€™ definition will however be handled directly inside the`MFrontNonlinearProblem`

class. - The module also relies on the
`NewtonSolver`

available in`FEniCS`

instead of programming the Newton method manually. - The module provides access to non-linear optimisation solvers (see the
`MFrontOptimisationProblem`

class). - Finite-strain behaviours are supported as well as parallel computations.

`MFrontNonlinearMaterial`

classThis class handles the loading of a `MFront`

behaviour through `MGIS`

. In particular, it will contain the following important attributes:

`behaviour`

: an instance of`MGIS`

`Behaviour`

class which handles all the information about a specific`MFront`

behaviour. It is created by the`load`

function which takes the path to a library, the name of a behaviour and a modelling hypothesis.`data_manager`

: an instance of`MGIS`

`MaterialDataManager`

class which handles a bunch of integration points. It is instantiated using`behaviour`

and the number of integration points

Note about finite strain behavioursBefore loading the behaviour, it is checked if the behaviour is a finite-strain one or not. In the former case, specific finite-strain options are used when calling

`load`

. Such options specify that the stress measure will be post-processed by`MGIS`

from Cauchy to First Piola Kirchhoff (PK1) stress and that the tangent operator will be given by \(\dfrac{\partial \text{PK1}}{\partial F}\) (`DPK1_DF`

).

A set of helper functions enables to retrieve information about `MFront`

object names and sizes of the form

`self.get_{1}_{2}`

where `{1}`

can be any of `material_property`

, `external_state_variable`

, `internal_state_variable`

, `gradient`

, `flux`

or `tangent_block`

and `{2}`

can be any of `names`

or `sizes`

.

`MFrontNonlinearProblem`

classThis class handles the definition and resolution of a nonlinear problem associated with a `MFrontNonlinearMaterial`

. Its main attributes are:

`u`

: the unknown mechanical field \(u\) (a`dolfin.Function`

)`material`

: the associated`MFrontNonlinearMaterial`

`state_variables`

: a dictionary of internal and external state variables. External state variables are represented as`dolfin`

objects whereas internal state variables are represented as Quadrature functions.`gradients`

: a dictionary of`Gradient`

objects`fluxes`

: a dictionary of`Flux`

objects`residual`

: the nonlinear residual \(F(u)\)`tangent_form`

: the tangent bilinear form associated with the residual`solver`

: the non-linear solver (default is`NewtonSolver`

)

Its main methods are:

`register_gradient`

: registers a`MFront`

gradient with a`UFL`

expression (see the registration concept).`register_external_state_variable`

: registers a`MFront`

external state variable with a`UFL`

expression (see the registration concept).`set_loading`

: defines the external forces linear form \(L(v)\) in the default residual expression`initialize`

: initializes the functions associated with gradients, fluxes, external and internal state variables objects and the corresponding tangent blocks. All gradients and external state variables must have been registered first. Automatic registration is performed at the beginning of the method call. This method has to be called once before calling`solve`

.`update_constitutive_law`

: performs the consitutive law update. Currently, this method is called for all quadrature points before assembling the residual or tangent form. Ideally, this should be done locally during the assembly (see the module current limitations).`solve`

: solves the associated non-linear problem using`solver`

`get_flux`

: returns the function associated with a flux`get_state_variable`

: returns the function associated with an internal state variable

By default, the nonlinear residual is assumed to take the following form: Find \(u\in V\) such that:

\[ \sum_{i=1}^p \int_{\Omega}\boldsymbol{\sigma}_i(u)\cdot\delta\mathbf{g}_i(v) \,\text{dx}- L(v) = 0 \quad \forall v\in V \qquad(1)\]

where the \(\boldsymbol{\sigma}_i(u)\) are a set of **fluxes** as defined in the `MFront`

behaviour using `@Flux`

and \(\mathbf{g}_i\) are the corresponding set of **gradients**, declared using `@Gradient`

. \(\delta \mathbf{g}_i(v)\) denotes the directional derivative of \(\mathbf{g}_i\) along direction \(v\), a `TestFunction`

of the function space \(V\), and \(L(v)\) is a linear form which can be expressed using standard `UFL`

operators. From the mechanical point of view, this residual expresses the balance between internal and external forces.

This generic form is suitable for most quasi-static mechanical problems but it does not necessarily encompass all possible situations. In particular, evolution equations such as transient heat transfer cannot be expressed in this form. However, this is not a limitation since the residual can also be redefined explicitly by the user (see Transient non-linear heat equation).

Either for the default or a user-defined one, the tangent operator associated with the residual is computed automatically using either the `UFL`

symbolic derivation `ufl.derivative`

for simple terms (e.g. involving the unknown field \(u\)) or using tangent operator blocks defined in the `MFront`

behaviour for the nonlinear fluxes and internal state variables. More precisely, the variation of each flux is given by:

\[ \dfrac{\partial \boldsymbol{\sigma}_i}{\partial u} = \sum_{j\in \text{blocks}(i)} \mathbb{T}^{\boldsymbol{\sigma}_i}_{\mathbf{g}_j}\cdot \delta \mathbf{g}_j(u) \]

where the flux \(\boldsymbol{\sigma}_i\) is assumed to depend on the gradients \(\mathbf{g}_j\) for \(j\in \text{blocks}(i)\) (which at least contains \(i\) itself in general) and \(\mathbb{T}^{\boldsymbol{\sigma}_i}_{\mathbf{g}_j} = \dfrac{\partial \boldsymbol{\sigma}_i}{\partial \mathbf{g}_j}\) is the tangent operator associated with the corresponding block. Variations of internal state variables are computed in the same manner.

In the default case, the tangent bilinear form therefore reads as:

\[ a_\text{tangent}(u, v) = \sum_{i=1}^p \sum_{j\in\text{blocks}(i)}\int_{\Omega}\delta\mathbf{g}_j(u)\cdot \mathbb{T}^{\boldsymbol{\sigma}_i}_{\mathbf{g}_j} \cdot\delta\mathbf{g}_i(v) \,\text{dx} \]

`Flux`

and `Gradient`

objectsTwo helper classes have been defined to handle flux and gradient objects:

- the
`Gradient`

class provides a representation of`MFront`

gradient objects. Its main purpose is to provide the corresponding`UFL`

expression, linking`MFront`

and`FEniCS`

concepts. It also handles:- the reshaping from
`UFL`

tensorial representation to`MFront`

vectorial conventions. - the symbolic expression of the gradient variation (directional derivative).
- the representation as a Quadrature function. This class is intended for internal use only. Gradient objects must be declared by the user using the registration concept.

- the reshaping from
- the
`Flux`

class provides a representation of`MFront`

flux objects. Its main purpose is to store the flux values in the form of a Quadrature function, as well as the tangent block structure of the flux, each of them also represented by a Quadrature function. - the
`InternalStateVariable`

class is the same as the`Flux`

class but for internal state variables. Both classes derive from the abstract`QuadratureFunction`

class.

Inspecting the default case (1), `MFront`

provides access to the flux names, shapes and values when performing the constitutive update and also to the corresponding gradient. The definition of the tangent operator blocks inside the `MFront`

behaviour also gives access to the block structure \(\text{blocks}(i)\) for each flux.

The remaining information which must be provided from the `FEniCS`

side are the unknown field \(u\) and its discretization space \(V\), the chosen integration measure \(\,\text{dx}\) (through the `quadrature_degree`

keyword) and, finally, the expression of each declared gradients \(\mathbf{g}_i\) in terms of the unknown field \(u\).

This step is what we call *registration* of each gradient which will be discussed in depth in the demos. Let us just mention that the gradients can registered using the `register_gradient`

method of the `MFrontNonlinearProblem`

class or via an automatic procedure if the gradient name matches predefined common gradient objects e.g. `"Strain"`

, `"TemperatureGradient"`

, etc.

`MFrontOptimisationProblem`

classThis class is a sibing to the `MFrontNonlinearProblem`

class. It enables to solve nonlinear optimisation problems, especially bound-constrained problems using the default `PETScTAOSolver`

of the form:

\[ \min_{b_l \leq u \leq b_u} f(u) \] where \(b_l\) (resp. \(b_u\)) denotes a lower (resp. upper) bound on the optimisation variable \(u\in V\).

By default, the objective function \(f(u)\) corresponds to the material total energy (stored + dissipated) computed from the `get_total_energy()`

method. The optimisation problem requires the definition of the gradient \(F(u)\) which, by default, corresponds to (1) and its jacobian which is computed as discussed before.

The module has been developed using `FEniCS`

version 2019.1.0. An important reimplementation will be planned once the `dolfinx`

project will officially release a stable version. In particular, it will aim at fixing the following current limitations:

- the constitutive update is performed before any assembly procedure of the tangent and residual forms. This adds an extra cost of looping over quadrature points and, more importantly, an important memory cost since all tangent blocks at all quadrature points must be saved. The
`dolfinx`

project should offer the possibility of using custom assemblers in which constitutive integration should be possible at the local assembly level - multi-materials are not completely supported yet. More precisely, spatially varying material properties can be defined but it is not possible to define two different constitutive behaviours on two different parts of the mesh. This feature has not been supported since it is not possible to define functions on sub-meshes yet. This should also be available soon in the next developments.
- memory transfers between
`FEniCS`

and`MGIS`

objects have not been optimized.