model_tutorial.md

• Introduction
• Description of the mechanical problem
• Resolution of the linear system
  • Spatial discretization
  • Time discretization
• Creating the new project
  • If you're not using XCode
  • With XCode
    • Create the target and the scheme
• Skeleton of the new model
  • Main file
  • Model files
  • Variational formulation files
  • InputData file
  • Fix header guards
  • First compilation
• Writing the model
  • InputData
  • The demo.lua file
    • First run of the executable
    • Filling the Lua file
      • Transient block
      • NumberingSubset1
      • Unknown1
      • Mesh10
      • Domains
        • Domain1
        • Domain2
        • Domain3
        • Domain10
      • EssentialBoundaryCondition1.
      • FiniteElementSpaces
        • FiniteElementSpace1
        • FiniteElementSpace2
      • Solid
        • Volumic mass \rho_0
        • Young modulus E
        • Poisson ratio \nu
        • Transient source.
        • Result
  • Running again... and commit!
  • Add variational_formulation object in the Model
  • Variational formulation: defining the Parameters
  • VariationalFormulation: add the numbering subset as a data attribute
  • Add method to run the static case
  • VariationalFormulation: allocating the system matrices and vectors
  • VariationalFormulation: define both operators (\underline{\underline{\mathbb{K}}} and \underline{\mathbb{F}} ) required for the static case
    • Source operator
    • Stiffness operator
  • Variational formulation: define a global matrix for stiffness
  • Assembling the stiffness
  • Static case (at last!)
  • Defining the mass operator
  • Additional global linear algebra
  • PrepareDynamicRun()
  • Model::Forward
  • FinalizeStep
• Add CMake file
  • Registering the new model in CMake
  • Building the libraries and executables
    • Define the new library:
    • Link to the MoReFEM library (or libraries)
    • Create the executable
    • Install the sources
    • EnsightOutput and tests

Introduction

The goal of this tutorial is to build step by step a rather simple model: elastic deformation of a bar upon which a force is applied; time step is constant and volumic mass is the same throughout the bar.

This tutorial assumes you're working in macOS and XCode; however the adaptation to make it work in Linux is very short (the most important change is to set up as soon as possible the CMake build, which is presented in the last step here).

The prerequisites are:

• Install MoReFEM.
• Initialize stuff such as git commit hook or installation of XCode templates:

python Scripts/init_morefem.py

• Create a branch named 0_demo_elasticity.

git checkout -b 0_demo_elasticity

This is a toy branch which is not aimed at being reintegrated in the master branch; but we will use it to showcase how git might be helpful while writing a new Model.

It should be stressed that new models are better created outside of the main MoReFEM library; however it adds (few) complications I would rather not address in a first approach (another tutorial will be written soon on this topic).

Description of the mechanical problem

In this tutorial we are solving a simple 2D problem: the deformation of a bar within the theory of linear elasticity using a finite element method. Here is a quick description of the equations depicting the problem and the resolution method that we will implement step by step; our model will in fact be generic enough to work with 3D cases.

The strong formulation of the equilibrium equations of a dynamic mechanical problem on a fixed domain \Omega_t reads:

\begin{cases}
\displaystyle  \rho(\underline{f}(\underline{x}) - \underline{\ddot{y}}(\underline{x})) + \underline{\nabla}_{\underline{x}} \cdot \underline{\underline{\sigma}}(\underline{x}) = 0 \quad \text{in $\Omega_t$} \\
\displaystyle \underline{\underline{\sigma}}(\underline{x}) \cdot \underline{n} = \underline{g}(\underline{x}) \quad \text{on } \Gamma^N \\
\displaystyle \underline{y}(\underline{x}) = \underline{0} \quad \text{on } \, \Gamma^D \\
\displaystyle \underline{\underline{\sigma}}^T(\underline{x}) = \underline{\underline{\sigma}}(\underline{x})
\end{cases}

The domain \Omega_t over which the solid is defined has a fixed boundary \Gamma^D where Dirichlet conditions are applied (here we impose a zero displacement condition) and another boundary \Gamma^N where Neumann conditions are applied (for instance external surfacic loadings in 3D). \underline{\underline{\sigma}} is the Cauchy stress tensor, \rho the volumic mass in the current configuration, \underline{f} represents the applied external forces, \underline{\ddot{y}} is the acceleration, \underline{n} is the outward-pointing normal of the boundary \Gamma^N and \underline{g} the external surfacic load.

The weak formulation of the fundamental law of dynamics in the deformed configuration reads:

\forall \underline{y}^* \in \mathcal{V}, \quad \int_{\Omega_t} \rho(\underline{f}-\underline{\ddot{y})} \cdot \underline{y}^* \textrm{d}\Omega_t + \int_{\Omega_t} \underline{\underline{\sigma}}:\underline{\nabla}_{\underline{x}}\underline{y}^* \textrm{d}\Omega_t  = \int_{\Gamma^N} \underline{g}\cdot \underline{y}^* \text{d}S 

In order to solve this equation without having to deal with the changes of positions due to the fact that each variable is expressed in the reference configuration, we can rewrite this equation within a Total Lagrangian formalism. Here we will only consider a surfacic loading \underline{g}(\underline{x}) without a volumic contribution (namely \underline{f}(\underline{x}) = \underline{0}). Without going into the details of the chain rules involved, the variational formulation of the fundamental law of dynamics in total Lagrangian formalism within elasticity theory reads:

\forall \underline{y}^* \in \mathcal{V}, \quad \displaystyle \int_{\Omega _{0}}^{} \rho_0 \underline{\ddot{y}} \cdot \underline{y}^{*}\textrm{d}\Omega_0 + \int_{\Omega_{0}}^{} \left( \underline{\underline{\underline{\underline{\text{A}}}}} : \underline{\underline{\varepsilon}}(\underline{y}) \right) : \underline{\underline{\varepsilon}}(\underline{y}^*) \textrm{d}\Omega_0
= \int_{\Gamma_0^N} \underline{g}_0 \cdot \underline{y}^* \text{d}S_0 

Where \underline{\underline{\varepsilon}}(\underline{y}) = \underline{\underline{\nabla}} \, \underline{y} + (\underline{\underline{\nabla}} \, \underline{y})^T is the linearized strain tensor and the consitutive behaviour law \underline{\underline{\sigma}} = \underline{\underline{\underline{\underline{\text{A}}}}} : \underline{\underline{\varepsilon}}(\underline{y}).

Here we are only solving for one unkown: the displacement \underline{y}. In order to do so, we will place ourselves within the plane strain hypothesis, namely \varepsilon_{33} = 0 which implies:

\underline{\underline{\sigma}} =  \begin{pmatrix} 
\sigma_{xx} & \sigma_{xy} & 0 & \\
\sigma_{xy} & \sigma_{yy} & 0 & \\
0 & 0 & \sigma_{zz} & 
\end{pmatrix} \quad \text{with} \quad \displaystyle \sigma_{zz} =  \frac{\lambda}{2(\mu+\lambda)}(\sigma_{xx}+\sigma_{yy})

where \lambda and \mu are the Lamé coefficients of the elastic solid.

Thus the constitutive behaviour law in 2D reads, using engineering notation:

\begin{pmatrix}
\sigma_{xx}(\underline{y}) & \\
\sigma_{yy}(\underline{y}) & \\
\sigma_{xy}(\underline{y}) & 
\end{pmatrix} =  \begin{pmatrix}
\lambda + 2\mu & \lambda & 0 & \\
\lambda  & \lambda + 2\mu  & 0 & \\
0 & 0 & \mu & 
\end{pmatrix}\begin{pmatrix}  
\varepsilon_{xx}({\underline{y}}) &\\
\varepsilon_{yy}({\underline{y}}) &\\
2\varepsilon_{xy}({\underline{y}}) &

\end{pmatrix} = \underline{\underline{\hat{\text{A}}}} \cdot \underline{\hat{\varepsilon}} (\underline{y})

Finally, we are left with the following set of equations and boundary conditions:

\begin{cases}
\displaystyle \forall \underline{y}^* \in \mathcal{V}, \quad \displaystyle \int_{\Omega _{0}}^{} \rho_0 \, \underline{y}^{*} \cdot  \underline{\ddot{y}} \, \textrm{d}\Omega_0 + \int_{\Omega_{0}}^{}   \underline{\hat{\varepsilon}}(\underline{y}^*)^T \cdot \underline{\underline{\hat{\text{A}}}} \cdot \underline{\hat{\varepsilon}}(\underline{y}) \, \textrm{d}\Omega_0
= \int_{\Gamma_0^N}  \underline{y}^* \cdot \underline{g}_0 \, \text{d}S_0 \\
\displaystyle \underline{y}(\underline{x}) = \underline{0} \quad \text{on } \, \Gamma^D 
\end{cases}

This bilinear system (with respect to the virtual displacement field \underline{y}^* and \underline{y}) will be solved using the finite element method.

Resolution of the linear system

Spatial discretization

In order to solve our previoulsy defined bilinear system (with respect to \underline{y} and \underline{y}^*), we will be using the standard Galerkin method. It consists of approximating the function of interest (the displacement field in our case) by a finite sum of known shape functions (polynomials usually) \phi_k(\underline{\xi}) weighted by unkown coefficients y_{jk} where k \in [1, \, N + 1] , N being the order of the shape functions used. In 2D, the discretization of the displacement field gives:

y_j = \sum_{k=1}^{N+1} y_{jk} \phi_k (\underline{\xi}) \quad j \in [x,y]  

\underline{y}(\underline{x})_{2D} = \begin{pmatrix} 
y_{x} & \\
y_{y} & 
\end{pmatrix} = \begin{pmatrix}
\phi_1 & \dots &  \phi_{N+1} & 0 & \dots & 0 &  \\
0 & \dots & 0 & \phi_1 & \dots & \phi_{N+1} & 
\end{pmatrix} \begin{pmatrix} y_{x1} & \\
\vdots & \\
y_{xN+1} & \\
y_{y1} & \\
\vdots & \\
y_{yN+1} & 
\end{pmatrix} = \underline{\underline{\mathbb{N}}} \cdot \underline{\mathbb{U}}_h 

To discretize the linearized strain tensor \underline{\underline{\varepsilon}} we can note that:

\hat{\underline{\varepsilon}}(\underline{y}) = \begin{pmatrix} \varepsilon_{xx} \\
\varepsilon_{yy} & \\
2\varepsilon_{xy} & 
\end{pmatrix}  = \begin{pmatrix} 
\partial_x y_x & \\
\partial_y y_y & \\
\partial_y y_x + \partial_x y_y &
\end{pmatrix} =
\begin{pmatrix} 
1 & 0 & 0 & 0 & \\
0 & 0 & 0 & 1 & \\
0 & 1 & 1 & 0 & 
\end{pmatrix} \begin{pmatrix} 
\partial_x y_x & \\
\partial_y y_x & \\
\partial_x y_y & \\
\partial_y y_y & \\
\end{pmatrix}

We can then apply the spatial discretization in a similar fashion using the derivatives of the shape functions:

\hat{\underline{\varepsilon}}(\underline{y}) = \begin{pmatrix} \varepsilon_{xx} \\
\varepsilon_{yy} & \\
2\varepsilon_{xy} & 
\end{pmatrix}  = 
\begin{pmatrix} 
1 & 0 & 0 & 0 & \\
0 & 0 & 0 & 1 & \\
0 & 1 & 1 & 0 & 
\end{pmatrix} \begin{pmatrix}
\partial_x \phi_1 & \dots &  \partial_x \phi_{N+1} & 0 & \dots & 0 &  \\
\partial_y \phi_1 & \dots &  \partial_y \phi_{N+1} & 0 & \dots & 0 &  \\
0 & \dots & 0 & \partial_x \phi_1 & \dots & \partial_x \phi_{N+1} & \\
0 & \dots & 0 & \partial_y \phi_1 & \dots & \partial_y \phi_{N+1} & 
\end{pmatrix} \begin{pmatrix} y_{x1} & \\
\vdots & \\
y_{xN+1} & \\
y_{y1} & \\
\vdots & \\
y_{yN+1} & 
\end{pmatrix} = \underline{\underline{\mathbb{B}}} \cdot \underline{\mathbb{U}}_h

Plugging these discretized forms into our equilibrium equation gives:

\forall \, \underline{\mathbb{U}}^*_h \in \mathcal{V}_h, \quad \displaystyle \int_{\Omega_0} \rho_0  \underline{\mathbb{U}}^{*T}_h \underline{\underline{\mathbb{N}}} ^T \cdot \underline{\underline{\mathbb{N}}} \cdot \underline{\dot{\mathbb{V}}}_h \, \text{d}\Omega_0 + \int_{\Omega_{0}}^{}  \underline{\mathbb{U}}^{*T}_h \cdot \underline{\underline{\mathbb{B}}}^T \cdot  \underline{\underline{\hat{\text{A}}}} \cdot \underline{\underline{\mathbb{B}}} \cdot \underline{\mathbb{U}}_h \, \textrm{d}\Omega_0
= \int_{\Gamma_0^N} \underline{\mathbb{U}}^{*T} \cdot \underline{\underline{\mathbb{N}}}^T \cdot \underline{g}_0 \, \text{d}S_0

where \underline{\dot{\mathbb{V}}}_h is time derivative of the unkown coefficients \dot{y}_{jk} relative to the velocity field (which is itself the time derivative of the unknown weighting coefficients of the displacement field).

This equation can be factorized and simplified as follows:

\forall \, \underline{\mathbb{U}}^*_h \in \mathcal{V}_h, \quad \displaystyle  \underline{\mathbb{U}}^{*T}_h \left[ \int_{\Omega_0} \rho_0  \underline{\underline{\mathbb{N}}} ^T \cdot \underline{\underline{\mathbb{N}}} \, \text{d}\Omega_0 \right] \underline{\dot{\mathbb{V}}}_h + \underline{\mathbb{U}}^{*T}_h \left[ \int_{\Omega_{0}}^{}  \underline{\underline{\mathbb{B}}}^T \cdot  \underline{\underline{\hat{\text{A}}}} \cdot \underline{\underline{\mathbb{B}}} \, \textrm{d}\Omega_0 \right] \underline{\mathbb{U}}_h
= \underline{\mathbb{U}}^{*T}_h \left[ \int_{\Gamma_0^N} \underline{\underline{\mathbb{N}}}^T \cdot \underline{g}_0 \, \text{d}S_0 \right]

\underline{\mathbb{U}}^{*T}_h  \cdot \underline{\underline{\mathbb{M}}} \cdot \underline{\dot{\mathbb{V}}}_h +  \underline{\mathbb{U}}^{*T}_h \cdot \underline{\underline{\mathbb{K}}} \cdot \underline{\mathbb{U}}_h = \underline{\mathbb{U}}^{*T}_h \cdot \underline{\mathbb{F}}

Finally, we have:

\underline{\underline{\mathbb{M}}} \cdot \underline{\dot{\mathbb{V}}}_h + \underline{\underline{\mathbb{K}}} \cdot \underline{\mathbb{U}}_h  = \underline{\mathbb{F}} 

where \underline{\underline{\mathbb{M}}} corresponds to the mass matrix, \underline{\underline{\mathbb{K}}} to the sitffness matrix and \underline{\mathbb{F}} is the discretized right-hand side (corresponding to a surfacic load). Each of these 3 quantities correspond to different operators that will be assembled in the variational formulation later in this tutorial.

Time discretization

Here the only time dependency for our elastic problem is the term associated to the inertia, involving the acceleration field \ddot{\underline{y}}:

\int_{\Omega _{0}}^{} \rho_0 \, \underline{y}^{*} \cdot  \underline{\ddot{y}} \, \textrm{d}\Omega_0

This means that in order to solve our system, we just need to update the acceleration values (no need to solve a linear system) with a selected time scheme , once we have the static solution. In this demo, we will implement the Newmark time scheme, which reads:

\frac{\dot{\underline{y}}^{n+1} + \dot{\underline{y}}^n}{2} = \frac{\underline{y}^{n+1} - \underline{y}^n}{\Delta t}

Creating the new project

We will create a new project in Sources/ModelInstances called DemoElasticity.

If you're not using XCode

You can still refer to what is indicated in the XCode section, with the following adaptations:

• Create the directory manually in Sources/ModelInstances.
• Set up the CMake presented later in this tutorial now.
• When a XCodeTemplate file is mentioned:

• Copy it (or them in case of a class) from XCodeTemplates directory into your working one, e.g. for Model class if you're in DemoElasticity directory:

cp ../../../XCodeTemplates/Model.xctemplate/FILEBASENAME.cpp Model.cpp cp ../../../XCodeTemplates/Model.xctemplate/FILEBASENAME.hpp Model.hpp cp ../../../XCodeTemplates/Model.xctemplate/FILEBASENAME.hxx Model.hxx

• Edit the newly copied files, and replace expressions such as ___VARIABLE_groupName:identifier___ by the value that was suggested for groupName (always 'DemoElasticity' for this tutorial).

With XCode

Create the target and the scheme

• Open the MoReFEM project in XCode.

• In XCode, double click on the MoReFEM at the top of the middle panel (besides the XCode icon). In Targets, click on one Model (typically Elasticity) and choose to duplicate it. Rename the duplicata DemoElasticity (default name is Elasticity copy).

• Click on the new target, and remove in the Build Phases tab all that is in the Compilation Sources block.

• Edit the scheme:

• Either press Command, shift and < simultaneously.
• Or through the menu: Product > Scheme
• Choose Manage schemes and click once on Elasticity copy; you will hence be able to rename it (DemoElasticity for instance).
• In Run > Arguments, add a line:

-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/DemoElasticity/demo.lua    

in Arguments Passed on Launch. This is really what would be provided to the command line in a terminal; it indicates the path of the Lua file to use (that should not exist yet in our case). This Lua file will be used to define all the parameters required by our model.

• Make sure DemoElasticity is the current scheme (on the top left of XCode)

• In the left panel right click on ModelInstances and choose New group; name it DemoElasticity.

Skeleton of the new model

Main file

• Right click on the newly created group, and choose New file:
• Choose Main in MoReFEM section (that should appear due to the init_morefem.py script)
• Choose DemoElasticity as both GroupName and ProblemName. It just provides default names for namespaces and relative paths; you might edit manually the source files if your choices were not right.
• Let main in Save As.
• In Targets, untick the default (probably 'Utilities') and tick the target you've just created (DemoElasticity).

Model files

• Repeat mostly what was done for the main file, except now you choose Model in MoReFEM section. Three files are created here with three extensions:

• cpp: code source compiled; definition of functions and methods.
• hpp: declaration of classes and functions with documentation; no implementation there.
• hxx: definition of functions and methods that are not suitable for cpp: template and inline.

Variational formulation files

• Repeat mostly what was done for the main file, except you now choose VariationalFormulation in MoReFEM section.

InputData file

• Repeat mostly what was done for the main file, except now you choose InputData in MoReFEM section.

Fix header guards

XCode templates aren't rich enough to handle properly header guards (we want a unique header guards in each hpp and hxx file); a dedicated Python script is present to do this task. This script should be used whenever you rename or move a header file so that the new header guard is properly generated (it uses its location and filename to generate the unique identifier).

In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:

python Scripts/header_guards.py  

You should see:

File ModelInstances/DemoElasticity/VariationalFormulation.hxx was modified for header guards.
File ModelInstances/DemoElasticity/Model.hpp was modified for header guards.
File ModelInstances/DemoElasticity/InputData.hpp was modified for header guards.
File ModelInstances/DemoElasticity/Model.hxx was modified for header guards.
File ModelInstances/DemoElasticity/VariationalFormulation.hpp was modified for header guards.

First compilation

Now you should be able to compile the code by:

• By pressing Command + B
• Or through the menu: Product > Build
• If you're using cmake, call to make or ninja (depending on the chosen generator) in a terminal.

We can therefore issue our first commit (with a dummy ticket number: we won't integrate these changes onto master):

In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 Skeleton of new model 'DemoElasticity' added." 

Writing the model

Now we can really start dealing with the model we intend to write. First step is to define the data the end user might provide in a Lua file; to do so we have to edit the InputData.hpp file and replace the default generated context by the one required by our model:

InputData

The point here is just to provide names more expressive in the code; for instance in the library we may use:

domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::neumann))

rather than:

domain_manager.GetDomain(14)

The former is longer but much more expressive, especially for complex models with many objects to manipulate. These objects will be designated by their index in the Lua file so we need to match the moniker used in the code to the actual index used in Lua file. The only rule in the index chosen is that the same can't be chosen twice for the same type of object (e.g. you can't define two meshes with the same index); the indexes will be reflected a bit later in the generated Lua file.

//! \copydoc doxygen_hide_mesh_enum
enum class MeshIndex
{
mesh = 10 // only one mesh considered in current model; we arbitrarily index it by a 10.
};


//! \copydoc doxygen_hide_domain_enum
enum class DomainIndex
{
highest_dimension = 1,
neumann = 2,
dirichlet = 3,
full_mesh = 10
};


//! \copydoc doxygen_hide_felt_space_enum
enum class FEltSpaceIndex
{
highest_dimension = 1, // doesn't matter the name is the same or not as in DomainIndex
neumann = 2
};


//! \copydoc doxygen_hide_unknown_enum
enum class UnknownIndex
{
displacement = 1
};


//! \copydoc doxygen_hide_solver_enum
enum class SolverIndex
{
solver = 1
};


//! \copydoc doxygen_hide_numbering_subset_enum
enum class NumberingSubsetIndex
{
monolithic = 1
};


//! \copydoc doxygen_hide_source_enum
enum class SourceIndex
{
surfacic = 1
};


//! \copydoc doxygen_hide_boundary_condition_enum
enum class BoundaryConditionIndex
{
sole = 1
};

These enum classes basically just define compilation constants; we will now use them to define the actual objects we will use in our model:

//! \copydoc doxygen_hide_input_parameter_tuple
using InputDataTuple = std::tuple
<
InputDataNS::TimeManager,

InputDataNS::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::monolithic)>,

InputDataNS::Unknown<EnumUnderlyingType(UnknownIndex::displacement)>,

InputDataNS::Mesh<EnumUnderlyingType(MeshIndex::mesh)>,

InputDataNS::Domain<EnumUnderlyingType(DomainIndex::highest_dimension)>,
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::neumann)>,
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::dirichlet)>,
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::full_mesh)>,

InputDataNS::DirichletBoundaryCondition<EnumUnderlyingType(BoundaryConditionIndex::sole)>,

InputDataNS::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::highest_dimension)>,
InputDataNS::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::neumann)>,

InputDataNS::Petsc<EnumUnderlyingType(SolverIndex::solver)>,

InputDataNS::Solid::VolumicMass,
InputDataNS::Solid::YoungModulus,
InputDataNS::Solid::PoissonRatio,
InputDataNS::Solid::PlaneStressStrain,

InputDataNS::VectorialTransientSource<EnumUnderlyingType(SourceIndex::surfacic)>,

InputDataNS::Result
>;

The code should compile as well, so:

In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:

```` git add Sources/ModelInstances/DemoElasticity git commit -m "#0 DemoElasticity: input data filled with the relevant data for the intended model." ````

The demo.lua file

First run of the executable

Let's run the code:

• By pressing Command + R
• Or through the menu: Product > Run
• If you're using cmake, use the generated demo_elasticity executable:

demo_elasticity -i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/DemoElasticity/demo.lua 

You should get the message:

${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/DemoElasticity/demo.lua wasn't existing and has just been created on processor 0; please edit it and then copy it onto all machines intended to run the code in parallel.
Exception caught from MoReFEMData<InputData>: Exception found at Sources/Core/MoReFEMData/MoReFEMData.hxx, line 150: Input parameter file to edit and complete!

This is normal: the path for the Lua file was valid but the file didn't exist yet; the choice was made in this case to create a Lua file with the blocks mentioned in the InputData.hpp file and default value written when possible.

This is not the case all the time, so you really need to edit this file and specify your input data.

In XCode, got to DemoElasticity group, right click on it and choose Add files to MoReFEM; demo.lua should appear there and you only need to double click on it.

Filling the Lua file

WARNING: This file must respect Lua syntax; you must put a , at the end of each line you fill (save for the last line in a block for which it is not mandatory - but there is no harm putting it nonetheless).

For each field, check the Expected format in the comment, which provides the way the entry should be provided.

If not specified otherwise, I indicate here only the fields that are to be changed; if a field appears in the Lua file that is not mentioned in the following, do not modify or remove it.

Transient block

This block sets the time interval over which the model is run. Choose timeMax = 0.5 for instance.

NumberingSubset1

There is only one "NumberingSubset" provided in the model; you just have to choose a name for it . For instance, you can call it name ="monolithic".

Unknown1

In this Elastic model only one unknown is considered: solid displacement \underline{y} which is a vectorial quantity.

So choose "displacement" (or whatever name you want: the only constraint is that two different unknowns cannot share the same name) and "vectorial".

name = "displacement",
nature = "vectorial"

Mesh10

Let's choose a toy mesh present in the library:

mesh = "${MOREFEM_ROOT}/Data/Mesh/elasticity_Nx50_Ny20_force_label.mesh",
dimension = 2

Domains

The domains we define are just a geometrical description of our mesh, they decribe how the mesh is partionned. This is required so that we can apply the different operators on the relevant domains.

Domain1

First domain is the one upon which the finite element space containing the 2D elements will be built. This corresponds to \displaystyle \Omega_0 \setminus \{ \Gamma^D \cup \Gamma^N \}.

mesh_index = { 10 }, -- related to Mesh1
dimension_list = { 2 }, -- only elements with dimension = 2
mesh_label_list = { }, -- no constraint upon MeshLabel
geometric_element_type_list = { } -- no constraint upon geometric element type

Domain2

Second Domain is the area upon which the Neumann boundary condition is applied. This corresponds to \Gamma^N.

mesh_index = { 10 }, -- related to Mesh1
dimension_list = { 1 }, -- only elements with dimension = 2
mesh_label_list = { 2 }, -- consider only edges with this mesh label
geometric_element_type_list = { } -- no constraint upon geometric element type

Domain3

Third Domain is the area upon which Dirichlet boundary condition is applied. This corresponds to \Gamma^D.

mesh_index = { 10 }, -- related to Mesh1
dimension_list = { }, -- no constraint upon dimension
mesh_label_list = { 1 }, -- consider only geometric elements with this mesh label
geometric_element_type_list = { } -- no constraint upon geometric element type

Domain10

Fourth domain is the whole mesh with nothing left aside. This corresponds to \Omega_0.

mesh_index = { 10 }, -- related to Mesh1
dimension_list = { }, -- no constraint upon dimension
mesh_label_list = { }, -- no constraint upon MeshLabel
geometric_element_type_list = { } -- no constraint upon geometric element type

FiniteElementSpaces

FiniteElementSpace1

The finite element space upon which the elastic operator is defined.

god_of_dof_index = 10, -- equivalent to the mesh_index
domain_index = 1, -- see Domain1
unknown_list = { "displacement" }, -- or whatever you have named Unknown1 earlier
shape_function_list = { "P1b" }, -- the moniker for the type of shape function you want. Putting a stupid value will result in an error message that will provides you all valid values. 
numbering_subset_list = { 1 }  -- the numbering subset upon which the unknown is numbered.

FiniteElementSpace2

The finite element space upon which the Neumann condition is defined. The Neumann condition itself is a specific variationnal term added to the right-hand side of our system, which is defined below as a Transient Source.

god_of_dof_index = 10, -- equivalent to the mesh_index
domain_index = 2, -- see Domain2
unknown_list = { "displacement" }, -- or whatever you have named Unknown1 earlier
shape_function_list = { "P1" }, -- not 'P1b': we're talking about segments here.
numbering_subset_list = { 1 } 

EssentialBoundaryCondition1.

This corresponds to \underline{y}(\underline{x}) = \underline{0} \quad \text{on} \quad \Gamma^D. To apply this boundary condition we modify the stiffness matrix of our system (either by pseudo-elimination or penalization) on the relevant degrees of freedom (those encompassed by Domain3 here).

name = "sole", 
component = "Comp12", -- means a value provided for X and another for Y
unknown = "displacement", -- or whatever you have named Unknown1 earlier
value = { 0., 0. }, -- the values for X and Y constrained by the Dirichlet boundary condition
domain_index = 3, -- The domain for the condition is Domain3

Solid

Volumic mass \rho_0

nature = "constant",
value = 1.3

Young modulus E

nature = "lua_function", -- for the demo: we'll just provide a constant function...
value = [[
        function (x, y, z)
        return 8307692.
        end
        ]]

Poisson ratio \nu

nature = "constant",
value = 0.04

Transient source.

This corresponds to \displaystyle \underline{g}_0(\underline{x}) \quad \text{in} \quad \int_{\Gamma_0^N} \underline{g}_0 \cdot \underline{y}^* \, \text{d}S_0

nature = { "constant", "constant", "constant" },
value = { 0., 5.e-3, 0. }

Result

output_directory = "${MOREFEM_RESULT_DIR}/DemoElasticity",

Running again... and commit!

With this file, the code should run smoothly and exit with exit code 0.

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 DemoElasticity: filling the demo.lua file."

Add variational_formulation object in the Model

In Model.hpp:

// At the top of Model.hpp with others include:
# include "ModelInstances/DemoElasticity/VariationalFormulation.hpp"

// At the bottom of the class definition
private:

//! Non constant access to the underlying VariationalFormulation object.
VariationalFormulation& GetNonCstVariationalFormulation() noexcept;

//! Access to the underlying VariationalFormulation object.
const VariationalFormulation& GetVariationalFormulation() const noexcept;


private:

//! Underlying variational formulation.
VariationalFormulation::unique_ptr variational_formulation_ = nullptr;

In Model.hxx:

inline const VariationalFormulation& Model::GetVariationalFormulation() const noexcept
{
assert(!(!variational_formulation_));
return *variational_formulation_;
}


inline VariationalFormulation& Model::GetNonCstVariationalFormulation() noexcept
{
return const_cast<VariationalFormulation&>(GetVariationalFormulation());
}

In Model.cpp, fill SupplInitialize():

void Model::SupplInitialize()
{
decltype(auto) god_of_dof = this->GetGodOfDof(EnumUnderlyingType(MeshIndex::mesh));
decltype(auto) morefem_data = parent::GetMoReFEMData();
decltype(auto) time_manager = parent::GetNonCstTimeManager();

{
decltype(auto) bc_manager = DirichletBoundaryConditionManager::GetInstance(__FILE__, __LINE__);
auto&& bc_list = { bc_manager.GetDirichletBoundaryConditionPtr("sole") };

variational_formulation_ = std::make_unique<VariationalFormulation>(morefem_data,
god_of_dof,
std::move(bc_list),
time_manager);
}

decltype(auto) variational_formulation = GetNonCstVariationalFormulation();
variational_formulation.Init(morefem_data.GetInputData());
}

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 DemoElasticity: variational formulation added in the Model."

Variational formulation: defining the Parameters

The MoReFEM library automatically loads the Lua file in the initialization process, but it doesn't mean all the content is instantly interpreted as directly usable MoReFEM objects. Some are (for instance the Mesh is fully built without further ado) but others like Parameter objects need to be built explicitly.

To begin with, by definition a Parameter is an object which purpose is to be evaluated at geometric coordinates; there are two flavours of such coordinates:

• Coordinates on the mesh.
• Coordinates in the reference element.

The data in the Lua file should explicitly specify which type of coords is handled in its comment.

In VariationalFormulation.hpp, in the class declaration:

At the top of the class, with using self, using parent, etc...:

//! Alias to the type of the source Parameter.
using source_parameter_type = Parameter<ParameterNS::Type::vector, LocalCoords>;

At the bottom of the class:

private:

//! Volumic mass.
const ScalarParameter<>& GetVolumicMass() const noexcept;

//! Young modulus.
const ScalarParameter<>& GetYoungModulus() const noexcept;

//! Poisson ratio.
const ScalarParameter<>& GetPoissonRatio() const noexcept;

//! Source parameter.
const source_parameter_type& GetSourceParameter() const noexcept;


private:

//! Volumic mass.
ScalarParameter<>::unique_ptr volumic_mass_ = nullptr;

//! Young modulus.
ScalarParameter<>::unique_ptr young_modulus_ = nullptr;

//! Poisson ratio.
ScalarParameter<>::unique_ptr poisson_ratio_ = nullptr;

//! Source parameter.
typename source_parameter_type::unique_ptr source_parameter_ = nullptr;

In VariationalFormulation.hxx:

inline const ScalarParameter<>& VariationalFormulation
::GetVolumicMass() const noexcept
{
assert(!(!volumic_mass_));
return *volumic_mass_;
}


inline const ScalarParameter<>& VariationalFormulation
::GetYoungModulus() const noexcept
{
assert(!(!young_modulus_));
return *young_modulus_;
}


inline const ScalarParameter<>& VariationalFormulation
::GetPoissonRatio() const noexcept
{
assert(!(!poisson_ratio_));
return *poisson_ratio_;
}


inline const VariationalFormulation::source_parameter_type& VariationalFormulation
::GetSourceParameter() const noexcept
{
assert(!(!source_parameter_));
return *source_parameter_;
}

In VariationalFormulation.cpp, fill SupplInit():

#include "Parameters/InitParameterFromInputData/InitParameterFromInputData.hpp"

void VariationalFormulation::SupplInit(const input_data_type& input_data)
{
decltype(auto) domain_manager =
DomainManager::GetInstance(__FILE__, __LINE__);

decltype(auto) full_mesh_domain = 
domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::full_mesh), __FILE__, __LINE__);

volumic_mass_ = InitScalarParameterFromInputData<InputDataNS::Solid::VolumicMass>("Volumic mass",
full_mesh_domain,
input_data);

if (!GetVolumicMass().IsConstant())
throw Exception("Current elastic model is restricted to a constant volumic mass!", __FILE__, __LINE__);

young_modulus_ = InitScalarParameterFromInputData<InputDataNS::Solid::YoungModulus>("Young modulus",
full_mesh_domain,
input_data);

poisson_ratio_ = InitScalarParameterFromInputData<InputDataNS::Solid::PoissonRatio>("Poisson ratio",
full_mesh_domain,
input_data);

decltype(auto) neumann_domain = 
domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::neumann), __FILE__, __LINE__);

source_parameter_ = InitThreeDimensionalParameterFromInputData
<
InputDataNS::VectorialTransientSource<EnumUnderlyingType(SourceIndex::surfacic)>
>("Surfacic source",
neumann_domain,
input_data);

}

In a terminal at the root of MoReFEM:

(after checking it compiles of course...):

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 DemoElasticity: parameters properly added to the variational formulation."

VariationalFormulation: add the numbering subset as a data attribute

We will have to use extensively the numbering subset in the variational formulation; it is obviously possible to get it from the GodOfDof each time but it is actually more practical to provide a direct access in the VariationalFormulation:

In VariationalFormulation.hpp:

Modify the constructor to add numbering subset as an argument: the most efficient way to store it is as a reference, and to do so we need to define it at the cobject construction:

/*!
* \copydoc doxygen_hide_varf_constructor
* \param[in] numbering_subset The only \a NumberingSubset considered in this model.
*/
explicit VariationalFormulation(const morefem_data_type& morefem_data,
const GodOfDof& god_of_dof,
DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
const NumberingSubset& numbering_subset,
TimeManager& time_manager);

and add the data attribute and a public accessor:

public:

//! Accessor to the numbering subset.
const NumberingSubset& GetNumberingSubset() const noexcept;

private:

//! Only numbering subset considered in this model.
const NumberingSubset& numbering_subset_;

In VariationalFormulation.hxx:

inline const NumberingSubset& VariationalFormulation::GetNumberingSubset() const noexcept
{
return numbering_subset_;
}

In VariationalFormulation.cpp, modify the constructor:

VariationalFormulation::VariationalFormulation(const morefem_data_type& morefem_data,
const GodOfDof& god_of_dof,
DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
const NumberingSubset& numbering_subset,
TimeManager& time_manager)
: parent(morefem_data,
time_manager,
god_of_dof,
std::move(boundary_condition_list)),
numbering_subset_(numbering_subset)
{ }

In Model.cpp, modify the construction of the variational formulation in SupplInitialize():

decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::monolithic));

variational_formulation_ = std::make_unique<VariationalFormulation>(morefem_data,
god_of_dof,
std::move(bc_list),
numbering_subset,
time_manager);

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 DemoElasticity: add numbering subset to the variational formulation."

Add method to run the static case

Let's provide a method in variational formulation to run the static case in the initialization; at the moment we will let it empty. This method will be used to solve the initial static problem \displaystyle \underline{\underline{\mathbb{K}}} \cdot \underline{\mathbb{U}}_h = \underline{\mathbb{F}} .

In VariationalFormulation.hpp, in the class:

public:

//! Run the static case.
void RunStaticCase();

In VariationalFormulation.cpp, in the class:

void VariationalFormulation::RunStaticCase()
{ }

In Model.cpp:

void Model::SupplInitialize()
{
...
decltype(auto) variational_formulation = GetNonCstVariationalFormulation();

variational_formulation.RunStaticCase();
variational_formulation.WriteSolution(time_manager, variational_formulation.GetNumberingSubset());
}

If you try to run the code, you will have an issue: the WriteSolution() method expects that the system solution is properly allocated, and this has to be done explicitly (because in complex models with several numbering subsets there is often no need to build all the possible configuration of matrices and vectors).

VariationalFormulation: allocating the system matrices and vectors

In VariationalFormulation.cpp, fill AllocateMatricesAndVectors():

void VariationalFormulation::AllocateMatricesAndVectors()
{
decltype(auto) numbering_subset = GetNumberingSubset();

parent::AllocateSystemMatrix(numbering_subset, numbering_subset);
parent::AllocateSystemVector(numbering_subset);
}

The code should now run properly.

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity/VariationalFormulation.cpp
git commit -m "#0 System linear algebra properly initialized; new method RunStaticCase() added in the VariationalFormulation (but still an empty shell at this stage)." 

VariationalFormulation: define both operators (\underline{\underline{\mathbb{K}}} and \underline{\mathbb{F}} ) required for the static case

Source operator

RHS for the system is just the surfacic source; we therefore need to define the related operator. This corresponds to the surfacic loading vector \displaystyle \underline{\mathbb{F}} = \int_{\Gamma_0^N} \underline{\underline{\mathbb{N}}}^T \cdot \underline{g}_0 \, \text{d}S_0:

In VariationalFormulation.hpp:

# include "OperatorInstances/VariationalOperator/LinearForm/TransientSource.hpp"

In the aliases at the top of the class:

private:

//! Alias to the type of the source operator.
using source_type_operator = GlobalVariationalOperatorNS::TransientSource<ParameterNS::Type::vector>;

And within the class:

private:

//! Surfacic source operator.
const source_type_operator& GetSurfacicSourceOperator() const noexcept;

private:

//! Surfacic source operator.
source_type_operator::const_unique_ptr surfacic_source_operator_ = nullptr;

In VariationalFormulation.hxx

inline const VariationalFormulation::source_type_operator& VariationalFormulation
::GetSurfacicSourceOperator() const noexcept
{
assert(!(!surfacic_source_operator_));
return *surfacic_source_operator_;
}

I also usually define a method InitializeOperators:

In VariationalFormulation.hpp (parameter will be used a bit later):

private:

/*!
* \brief Initialize the operators.
*
* \copydoc doxygen_hide_input_data_arg
*/
void InitializeOperators(const input_data_type& input_data);

In VariationalFormulation.cpp:

• Add a new line at the end of SupplInit:

void VariationalFormulation::SupplInit(const input_data_type& input_data)
{
... 

InitializeOperators(input_data);
}

• And provide the implementation for the new method:

void VariationalFormulation::InitializeOperators(const input_data_type& input_data)
{
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) unknown_manager = UnknownManager::GetInstance(__FILE__, __LINE__);

decltype(auto) displacement_ptr =
unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::displacement));

{
decltype(auto) felt_space_neumann =
god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::neumann));

decltype(auto) source_parameter = GetSourceParameter();

surfacic_source_operator_ =
std::make_unique<source_type_operator>(felt_space_neumann,
displacement_ptr,
source_parameter);
}
}

Stiffness operator

We also need to define the sitffness matrix. This corresponds to \displaystyle \underline{\underline{\mathbb{K}}} = \int_{\Omega_{0}}^{} \underline{\underline{\mathbb{B}}}^T \cdot \underline{\underline{\hat{\text{A}}}} \cdot \underline{\underline{\mathbb{B}}} \, \textrm{d}\Omega_0

In VariationalFormulation.hpp:

# include "OperatorInstances/VariationalOperator/BilinearForm/GradOnGradientBasedElasticityTensor.hpp"

private:

//! Get the stiffness operator.
const GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor& GetStiffnessOperator() const noexcept;

private:

//! Stiffness operator.
GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor::const_unique_ptr
stiffness_operator_ = nullptr;

In VariationalFormulation.hxx

inline const GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor&
VariationalFormulation::GetStiffnessOperator() const noexcept
{
assert(!(!stiffness_operator_));
return *stiffness_operator_;
}

In VariationalFormulation.cpp, complete InitializeOperators():

void VariationalFormulation::InitializeOperators(const input_data_type& input_data)
{
... 

decltype(auto) felt_space_highest_dimension =
god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::highest_dimension));

decltype(auto) mesh_dimension = god_of_dof.GetMesh().GetDimension();

const auto configuration =
ParameterNS::ReadGradientBasedElasticityTensorConfigurationFromFile(mesh_dimension,
input_data);

stiffness_operator_ = std::make_unique<GlobalVariationalOperatorNS
::GradOnGradientBasedElasticityTensor>(felt_space_highest_dimension,
displacement_ptr,
displacement_ptr,
GetYoungModulus(),
GetPoissonRatio(),
configuration);
}

In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 Add both operators required to run the static case."

Variational formulation: define a global matrix for stiffness

Additionally to the system matrices, it is convenient to add work matrices that might be reused and spare some recomputation. For instance it is useful to compute once and for all the stiffness matrix (given we consider in our model a constant volumic mass and time step); so let's add this matrix:

In VariationalFormulation.hpp:

private:

//! Accessor to the stiffness matrix.
const GlobalMatrix& GetStiffnessMatrix() const noexcept;

//! Non constant accessor to the stiffness matrix.
GlobalMatrix& GetNonCstStiffnessMatrix() noexcept;

private:

//! Stiffness matrix.
GlobalMatrix::unique_ptr stiffness_matrix_ = nullptr;

In VariationalFormulation.hxx

inline const GlobalMatrix& VariationalFormulation
::GetStiffnessMatrix() const noexcept
{
assert(!(!stiffness_matrix_));
return *stiffness_matrix_;
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstStiffnessMatrix() noexcept
{
return const_cast<GlobalMatrix&>(GetStiffnessMatrix());
}

In VariationalFormulation.cpp:

void VariationalFormulation
::AllocateMatricesAndVectors()
{
...

decltype(auto) system_matrix = parent::GetSystemMatrix(numbering_subset, numbering_subset);

stiffness_matrix_ = std::make_unique<GlobalMatrix>(system_matrix);
}

This uses up the copy constructor of the GlobalMatrix: system_matrix has been properly initialized beforehand and we reuse the result rather than calling again all the Petsc stuff to initialize properly this (Petsc) matrix.

Assembling the stiffness

Now we're able to assemble the stiffness operator into the stiffness matrix:

In VariationalFormulation.cpp:

void VariationalFormulation
::SupplInit(const InputData& input_data)
{
...

{
// Assemble the stiffness matrix.
GlobalMatrixWithCoefficient matrix(GetNonCstStiffnessMatrix(), 1.);
GetStiffnessOperator().Assemble(std::make_tuple(std::ref(matrix)));
}
}

The syntax is a tad heavy here, but it's the price to pay for flexibility and efficiency: we may assemble into several matrices in a single command (or matrices and vectors for non linear operators).

Static case (at last!)

We are now finally able to run the whole static case, sovling:

\begin{cases}
\displaystyle \underline{\underline{\mathbb{K}}} \cdot \underline{\mathbb{U}}_h  = \underline{\mathbb{F}} \\
\underline{y}(\underline{x}) = \underline{0} \quad \text{on} \quad \Gamma^D
\end{cases}

In VariationalFormulation.cpp:

void VariationalFormulation::RunStaticCase()
{
decltype(auto) numbering_subset = GetNumberingSubset();

// Assembling transient source into system RHS.
{
constexpr double irrelevant_time = 0.; // as this parameter has no time dependency. 
GlobalVectorWithCoefficient vector(GetNonCstSystemRhs(numbering_subset), 1.);
GetSurfacicSourceOperator().Assemble(std::make_tuple(std::ref(vector)), irrelevant_time);
}

parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset).Copy(GetStiffnessMatrix(), __FILE__, __LINE__);


ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_matrix_and_rhs>(numbering_subset,
numbering_subset);
SolveLinear<IsFactorized::no>(numbering_subset, numbering_subset, __FILE__, __LINE__);

parent::DebugPrintSolutionElasticWithOneUnknown(numbering_subset, 10);
}

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity/VariationalFormulation.cpp
git commit -m "#0 Static case implemented and working." 

Defining the mass operator

Here we are defining the mass matrix required for the dynamic part of the run \underline{\underline{\mathbb{M}}} = \int_{\Omega_0} \rho_0 \underline{\underline{\mathbb{N}}} ^T \cdot \underline{\underline{\mathbb{N}}} \, \text{d}\Omega_0

In VariationalFormulation.hpp:

# include "OperatorInstances/VariationalOperator/BilinearForm/Mass.hpp"

private:

//! Get the mass per square time step operator.
const GlobalVariationalOperatorNS::Mass& GetMassOperator() const noexcept;

private:

//! Mass operator.
GlobalVariationalOperatorNS::Mass::const_unique_ptr mass_operator_ = nullptr;

In VariationalFormulation.hxx:

inline const GlobalVariationalOperatorNS::Mass& VariationalFormulation
::GetMassOperator() const noexcept
{
assert(!(!mass_operator_));
return *mass_operator_;
}

In VariationalFormulation.cpp, complete InitializeOperators():

void VariationalFormulation::InitializeOperators(const InputData& input_data)
{
...
mass_operator_ = std::make_unique<GlobalVariationalOperatorNS::Mass>(felt_space_highest_dimension,
displacement_ptr,
displacement_ptr);
}

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity/VariationalFormulation.cpp
git commit -m "#0 Mass operator added." 

Additional global linear algebra

In our simple model, dynamic iterations do not need a call to a solver: basic matrix / vector operators are enough. Indeed, we just need to update the values of the velocity and displacement fields as the problem we are solving here has only one dynamic dependency (related to the inertia). This is done with the following time scheme (Newmark): \frac{\dot{\underline{y}}^{n+1} + \dot{\underline{y}}^n}{2} = \frac{\underline{y}^{n+1} - \underline{y}^n}{\Delta t}. We will define helpful matrices and vectors to limit at maximum recomputation (this will get a tad verbose with all accessors):

In VariationalFormulation.hpp:

private:

/// \name Accessors to vectors and matrices specific to the elastic problem.
///@{

//! Accessor to the \a GlobalVector which contains current displacement.
const GlobalVector& GetVectorCurrentDisplacement() const noexcept;

//! Non constant accessor to the \a GlobalVector which contains current displacement.
GlobalVector& GetNonCstVectorCurrentDisplacement() noexcept;

//! Accessor to the \a GlobalVector which contains current velocity.
const GlobalVector& GetVectorCurrentVelocity() const noexcept;

//! Non constant accessor to the \a GlobalVector which contains current velocity.
GlobalVector& GetNonCstVectorCurrentVelocity() noexcept;

//! Accessor to the \a GlobalMatrix used along displacement in the model.
const GlobalMatrix& GetMatrixCurrentDisplacement() const noexcept;

//! Non constant accessor to the \a GlobalMatrix used along displacement in the model.
GlobalMatrix& GetNonCstMatrixCurrentDisplacement() noexcept;

//! Accessor to the \a GlobalMatrix used along velocity in the model.
const GlobalMatrix& GetMatrixCurrentVelocity() const noexcept;

//! Non constant accessor to the \a GlobalMatrix used along velocity in the model.
GlobalMatrix& GetNonCstMatrixCurrentVelocity() noexcept;

//! Accessor to the mass matrix.
const GlobalMatrix& GetMassMatrix() const noexcept;

//! Non constant accessor to the mass matrix.
GlobalMatrix& GetNonCstMassMatrix() noexcept;

///@}

private:

//! Vector current displacement.
GlobalVector::unique_ptr vector_current_displacement_ = nullptr;

//! Vector current velocity.
GlobalVector::unique_ptr vector_current_velocity_ = nullptr;

//! Matrix current displacement.
GlobalMatrix::unique_ptr matrix_current_displacement_ = nullptr;

//! Matrix current velocity.
GlobalMatrix::unique_ptr matrix_current_velocity_ = nullptr;

//! Mass matrix.
GlobalMatrix::unique_ptr mass_matrix_ = nullptr;

In VariationalFormulation.hxx

inline const GlobalVector& VariationalFormulation
::GetVectorCurrentDisplacement() const noexcept
{
assert(!(!vector_current_displacement_));
return *vector_current_displacement_;
}


inline GlobalVector& VariationalFormulation
::GetNonCstVectorCurrentDisplacement() noexcept
{
return const_cast<GlobalVector&>(GetVectorCurrentDisplacement());
}


inline const GlobalVector& VariationalFormulation
::GetVectorCurrentVelocity() const noexcept
{
assert(!(!vector_current_velocity_));
return *vector_current_velocity_;
}


inline GlobalVector& VariationalFormulation
::GetNonCstVectorCurrentVelocity() noexcept
{
return const_cast<GlobalVector&>(GetVectorCurrentVelocity());
}


inline const GlobalMatrix& VariationalFormulation
::GetMatrixCurrentDisplacement() const noexcept
{
assert(!(!matrix_current_displacement_));
return *matrix_current_displacement_;
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstMatrixCurrentDisplacement() noexcept
{
return const_cast<GlobalMatrix&>(GetMatrixCurrentDisplacement());
}


inline const GlobalMatrix& VariationalFormulation
::GetMatrixCurrentVelocity() const noexcept
{
assert(!(!matrix_current_velocity_));
return *matrix_current_velocity_;
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstMatrixCurrentVelocity() noexcept
{
return const_cast<GlobalMatrix&>(GetMatrixCurrentVelocity());
}


inline const GlobalMatrix& VariationalFormulation
::GetMassMatrix() const noexcept
{
assert(!(!mass_matrix_));
return *mass_matrix_;
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstMassMatrix() noexcept
{
return const_cast<GlobalMatrix&>(GetMassMatrix());
}

In VariationalFormulation.cpp:

void VariationalFormulation
::AllocateMatricesAndVectors()
{
...

mass_matrix_ = std::make_unique<GlobalMatrix>(system_matrix);
matrix_current_displacement_ = std::make_unique<GlobalMatrix>(system_matrix);
matrix_current_velocity_ = std::make_unique<GlobalMatrix>(system_matrix);

decltype(auto) system_rhs = parent::GetSystemRhs(numbering_subset);    
vector_current_velocity_ = std::make_unique<GlobalVector>(system_rhs);
vector_current_displacement_ = std::make_unique<GlobalVector>(system_rhs);
}

In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 Add all linear algebra required to deal with dynamic iterations." 

PrepareDynamicRun()

In Model.cpp:

void Model::SupplInitialize()
{
...
variational_formulation.PrepareDynamicRuns();
}

In VariationalFormulation.hpp:

public:

/*!
* \brief Prepare dynamic runs.
*
* For instance for dynamic iterations the system matrix is always the same; compute it once and for all
* here.
*
* StaticOrDynamic rhs is what changes between two time iterations, but to compute it the same matrices are used at
* each time iteration; they are also computed there.
*/
void PrepareDynamicRuns();

private:

//! Compute all the matrices required for dynamic calculation.
void ComputeDynamicMatrices();

//! Update the displacement for the next time iteration.
void UpdateDisplacement();   

In VariationalFormulation.cpp:

void VariationalFormulation::PrepareDynamicRuns()
{
// Assemble once and for all the system matrix in dynamic case; intermediate matrices used
// to compute rhs at each time iteration are also computed there.
ComputeDynamicMatrices();
UpdateDisplacement();
}

void VariationalFormulation::ComputeDynamicMatrices()
{
decltype(auto) numbering_subset = GetNumberingSubset();
decltype(auto) system_matrix = this->GetNonCstSystemMatrix(numbering_subset, numbering_subset);
decltype(auto) stiffness = GetStiffnessMatrix();

{
GlobalMatrixWithCoefficient mass(GetNonCstMassMatrix(), 1.);
GetMassOperator().Assemble(std::make_tuple(std::ref(mass)));
}

decltype(auto) mass = GetMassMatrix();

{
// Compute the system matrix, which won't change afterwards!
system_matrix.Copy(stiffness, __FILE__, __LINE__);
system_matrix.Scale(0.5, __FILE__, __LINE__);

const auto coefficient =
2. * GetVolumicMass().GetConstantValue() / NumericNS::Square(parent::GetTimeManager().GetTimeStep());

#ifndef NDEBUG
AssertSameNumberingSubset(mass, system_matrix);
#endif // NDEBUG

Wrappers::Petsc::AXPY<NonZeroPattern::same>( coefficient,
mass,
system_matrix,
__FILE__, __LINE__);
}

{
// Displacement matrix.
decltype(auto) current_displacement_matrix = GetNonCstMatrixCurrentDisplacement();
current_displacement_matrix.Copy(mass, __FILE__, __LINE__);

const auto coefficient =
2. * GetVolumicMass().GetConstantValue() / NumericNS::Square(parent::GetTimeManager().GetTimeStep());

current_displacement_matrix.Scale(coefficient, __FILE__, __LINE__);

#ifndef NDEBUG
AssertSameNumberingSubset(stiffness, current_displacement_matrix);
#endif // NDEBUG

Wrappers::Petsc::AXPY<NonZeroPattern::same>( -.5,
stiffness,
current_displacement_matrix,
__FILE__, __LINE__);

}

{
// Velocity matrix.
decltype(auto) current_velocity_matrix = GetNonCstMatrixCurrentVelocity();
current_velocity_matrix.Copy(mass, __FILE__, __LINE__);
current_velocity_matrix.Scale(2. * GetVolumicMass().GetConstantValue() / parent::GetTimeManager().GetTimeStep(),
__FILE__, __LINE__);
}
}

void VariationalFormulation::UpdateDisplacement()
{
decltype(auto) numbering_subset = GetNumberingSubset();

GetNonCstVectorCurrentDisplacement().Copy(parent::GetSystemSolution(numbering_subset),
__FILE__, __LINE__);
}

In a terminal at the root of MoReFEM:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 DemoElasticity: quantities required for the dynamic runs have been computed in the initialization phase."

Model::Forward

This is where we are solving our dynamic linear system:

\begin{cases}
\underline{\underline{\mathbb{M}}} \cdot \underline{\dot{\mathbb{V}}}_h + \underline{\underline{\mathbb{K}}} \cdot \underline{\mathbb{U}}_h  = \underline{\mathbb{F}} \\
\underline{y}(\underline{x}) = \underline{0} \quad \text{on} \quad \Gamma^D
\end{cases}

In VariationalFormulation.hpp:

public:

//! At each time iteration, compute the system Rhs.
void ComputeDynamicSystemRhs();

In VariationalFormulation.cpp:

void VariationalFormulation::ComputeDynamicSystemRhs()
{
decltype(auto) numbering_subset = GetNumberingSubset();

decltype(auto) rhs = this->GetNonCstSystemRhs(numbering_subset);

// Compute the system RHS. The rhs is effectively zeroed through the first MatMult call.
decltype(auto) current_displacement_matrix = GetMatrixCurrentDisplacement();
decltype(auto) current_velocity_matrix = GetMatrixCurrentVelocity();

decltype(auto) current_displacement_vector = GetVectorCurrentDisplacement();
decltype(auto) current_velocity_vector = GetVectorCurrentVelocity();

Wrappers::Petsc::MatMult(current_displacement_matrix, current_displacement_vector, rhs, __FILE__, __LINE__);
Wrappers::Petsc::MatMultAdd(current_velocity_matrix, current_velocity_vector, rhs, rhs, __FILE__, __LINE__);
}

In Model.cpp, complete Forward():

void Model::Forward()
{
decltype(auto) formulation = this->GetNonCstVariationalFormulation();            

// Only Rhs is modified at each time iteration; compute it and solve the system.
formulation.ComputeDynamicSystemRhs();
decltype(auto) numbering_subset = formulation.GetNumberingSubset();

if (GetTimeManager().NtimeModified() == 1)
{
formulation.ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_matrix_and_rhs>(numbering_subset, numbering_subset);
formulation.SolveLinear<IsFactorized::no>(numbering_subset, numbering_subset, __FILE__, __LINE__);
}
else
{
formulation.ApplyEssentialBoundaryCondition<VariationalFormulationNS::On::system_rhs>(numbering_subset, numbering_subset);
formulation.SolveLinear<IsFactorized::yes>(numbering_subset, numbering_subset, __FILE__, __LINE__);
}
}

FinalizeStep

In Model.cpp, complete FinalizeStep():

void Model::SupplFinalizeStep()
{
// Update quantities for next iteration.
decltype(auto) formulation = this->GetNonCstVariationalFormulation();
decltype(auto) numbering_subset = formulation.GetNumberingSubset();

formulation.DebugPrintNorm(numbering_subset);
formulation.DebugPrintSolutionElasticWithOneUnknown(numbering_subset, 10);
formulation.WriteSolution(this->GetTimeManager(), numbering_subset);

formulation.UpdateDisplacementAndVelocity();
}

In VariationalFormulation.hpp:

public:
//! Update displacement and velocity for next time step.
void UpdateDisplacementAndVelocity();

private:
/*!
* \brief Update the velocity for the next time iteration.
*
* BEWARE: this method must be called BEFORE UpdateDisplacement(), as it relies upon the displacement
* that has been used in the last Ksp solve.
*/
void UpdateVelocity();    

In VariationalFormulation.cpp:

void VariationalFormulation::UpdateDisplacementAndVelocity()
{
// Ordering of both calls is capital here!
UpdateVelocity();
UpdateDisplacement();
}


void VariationalFormulation::UpdateVelocity()
{
decltype(auto) current_displacement_vector = GetVectorCurrentDisplacement();
decltype(auto) system_solution = parent::GetSystemSolution(GetNumberingSubset());
decltype(auto) current_velocity_vector = GetNonCstVectorCurrentVelocity();

assert(parent::GetTimeManager().GetStaticOrDynamic() == StaticOrDynamic::dynamic_);

{
// Update first the velocity.
Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_only> solution(system_solution, __FILE__, __LINE__);
Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_only> displacement_prev(current_displacement_vector, __FILE__, __LINE__);
Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_and_write> velocity(current_velocity_vector, __FILE__, __LINE__);

const unsigned int size = velocity.GetSize(__FILE__, __LINE__);
assert(size == solution.GetSize(__FILE__, __LINE__));
assert(size == displacement_prev.GetSize(__FILE__, __LINE__));

const double factor = 2. / parent::GetTimeManager().GetTimeStep();

for (unsigned int i = 0u; i < size; ++i)
{
velocity[i] *= -1.;
velocity[i] += factor * (solution.GetValue(i) - displacement_prev.GetValue(i));
}
}
}

In a terminal at the root of MoReFEM project:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 DemoElasticity executable fully implemented."

Add CMake file

Registering the new model in CMake

To enable compilation in command line, you should provide a CMakeLists.txt file which should be added itself in the ModelInstances CMakeLists.txt with the line:

include(${CMAKE_CURRENT_LIST_DIR}/DemoElasticity/CMakeLists.txt)

Building the libraries and executables

In your CMakeLists.txt:

Define the new library:

add_library(MoReFEM4DemoElasticity_lib ${LIBRARY_TYPE} "")

target_sources(MoReFEM4DemoElasticity_lib
    PRIVATE
        "${CMAKE_CURRENT_LIST_DIR}/Model.cpp" / 
        "${CMAKE_CURRENT_LIST_DIR}/VariationalFormulation.cpp"
	PUBLIC
		"${CMAKE_CURRENT_LIST_DIR}/Model.hpp" /
		"${CMAKE_CURRENT_LIST_DIR}/Model.hxx" /
		"${CMAKE_CURRENT_LIST_DIR}/InputData.hpp" /
 		"${CMAKE_CURRENT_LIST_DIR}/VariationalFormulation.hpp" /
		"${CMAKE_CURRENT_LIST_DIR}/VariationalFormulation.hxx"
)

Link to the MoReFEM library (or libraries)

target_link_libraries(MoReFEM4DemoElasticity_lib
                      ${ALL_LOAD_BEGIN_FLAG}
                      ${MOREFEM_MODEL}
                      ${ALL_LOAD_END_FLAG})

${MOREFEM_MODEL} mihght point to just the Model librry of MoReFEM one, depending whether the library was built as one monolithic one or several ones.

Create the executable

add_executable(MoReFEM4DemoElasticity ${CMAKE_CURRENT_LIST_DIR}/main.cpp)
target_link_libraries(MoReFEM4DemoElasticity
                      MoReFEM4DemoElasticity_lib)
                      
apply_lto_if_supported(MoReFEM4DemoElasticity)                            

This last command enables so-called LTO optimizations if it is supported but the operating system and the compiler.

Install the sources

morefem_install(MoReFEM4DemoElasticity MoReFEM4DemoElasticity_lib)

This line makes the make install (or ninja install) do its bidding: putting the executable and the library in the target directory outside of the sources.

EnsightOutput and tests

We usually provide with the models:

• An executable to produce an Ensight-compatible output in post-processing (which may be opened by Ensight or Paraview)
• Tests that checks the numerical results stay constant over time. To take into account numerical uncertainties and variability (especially with parallel runs!), what is compared is the output provided by Ensight, which are less precise than the original outputs.

Specific notes will be written to explain how to do so: this one is already long enough!

In a terminal at the root of MoReFEM project:

git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 CMake build has been introduced properly. You have reached the end iof this tutorial, congratulations\!"