model_tutorial.md 68 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11

<!-- toc -->

- [Introduction](#introduction)
- [Description of the mechanical problem](#description-of-the-mechanical-problem)
- [Resolution of the linear system](#resolution-of-the-linear-system)
  * [Spatial discretization](#spatial-discretization)
  * [Time discretization](#time-discretization)
- [Creating the new project](#creating-the-new-project)
  * [If you're not using XCode](#if-youre-not-using-xcode)
  * [With XCode](#with-xcode)
12
    + [Create the target and the scheme](#create-the-target-and-the-scheme)
13
14
15
16
17
18
19
20
21
22
- [Skeleton of the new model](#skeleton-of-the-new-model)
  * [Main file](#main-file)
  * [Model files](#model-files)
  * [Variational formulation files](#variational-formulation-files)
  * [InputData file](#inputdata-file)
  * [Fix header guards](#fix-header-guards)
  * [First compilation](#first-compilation)
- [Writing the model](#writing-the-model)
  * [InputData](#inputdata)
  * [The demo.lua file](#the-demolua-file)
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
    + [First run of the executable](#first-run-of-the-executable)
    + [Filling the Lua file](#filling-the-lua-file)
      - [Transient block](#transient-block)
      - [NumberingSubset1](#numberingsubset1)
      - [Unknown1](#unknown1)
      - [Mesh10](#mesh10)
      - [Domains](#domains)
        * [Domain1](#domain1)
        * [Domain2](#domain2)
        * [Domain3](#domain3)
        * [Domain10](#domain10)
      - [EssentialBoundaryCondition1.](#essentialboundarycondition1)
      - [FiniteElementSpaces](#finiteelementspaces)
        * [FiniteElementSpace1](#finiteelementspace1)
        * [FiniteElementSpace2](#finiteelementspace2)
      - [Solid](#solid)
        * [Volumic mass $` \rho_0 `$](#volumic-mass--rho_0-)
        * [Young modulus $`E`$](#young-modulus-e)
        * [Poisson ratio $`\nu`$](#poisson-ratio-nu)
        * [Transient source.](#transient-source)
        * [Result](#result)
44
45
46
47
48
49
50
  * [Running again... and commit!](#running-again-and-commit)
  * [Add variational_formulation object in the Model](#add-variational_formulation-object-in-the-model)
  * [Variational formulation: defining the Parameters](#variational-formulation-defining-the-parameters)
  * [VariationalFormulation: add the numbering subset as a data attribute](#variationalformulation-add-the-numbering-subset-as-a-data-attribute)
  * [Add method to run the static case](#add-method-to-run-the-static-case)
  * [VariationalFormulation: allocating the system matrices and vectors](#variationalformulation-allocating-the-system-matrices-and-vectors)
  * [VariationalFormulation: define both operators ($`\underline{\underline{\mathbb{K}}}`$ and $` \underline{\mathbb{F}}`$ ) required for the static case](#variationalformulation-define-both-operators-underlineunderlinemathbbk-and--underlinemathbbf---required-for-the-static-case)
51
52
    + [Source operator](#source-operator)
    + [Stiffness operator](#stiffness-operator)
53
54
55
56
57
58
59
60
61
62
63
  * [Variational formulation: define a global matrix for stiffness](#variational-formulation-define-a-global-matrix-for-stiffness)
  * [Assembling the stiffness](#assembling-the-stiffness)
  * [Static case (at last!)](#static-case-at-last)
  * [Defining the mass operator](#defining-the-mass-operator)
  * [Additional global linear algebra](#additional-global-linear-algebra)
  * [PrepareDynamicRun()](#preparedynamicrun)
  * [Model::Forward](#modelforward)
  * [FinalizeStep](#finalizestep)
- [Add CMake file](#add-cmake-file)
  * [Registering the new model in CMake](#registering-the-new-model-in-cmake)
  * [Building the libraries and executables](#building-the-libraries-and-executables)
64
65
66
67
68
    + [Define the new library:](#define-the-new-library)
    + [Link to the MoReFEM library (or libraries)](#link-to-the-morefem-library-or-libraries)
    + [Create the executable](#create-the-executable)
    + [Install the sources](#install-the-sources)
    + [EnsightOutput and tests](#ensightoutput-and-tests)
69
70
71

<!-- tocstop -->

72
73
74
75
# 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.  

76
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).  
77
78
79
80
81
82
83
84
85


The prerequisites are:

- Install [MoReFEM](https://gitlab.inria.fr/MoReFEM/CoreLibrary/MoReFEM).
- Initialize stuff such as git commit hook or installation of XCode templates:
````
python Scripts/init_morefem.py
````
86

87
- Create a branch named 0_demo_elasticity.
88

89
````
90
git checkout -b 0_demo_elasticity
91
92
93
94
95
96
97
````

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).


98
99
# Description of the mechanical problem

100
In this tutorial we are solving a simple <strong>2D problem</strong>: the deformation of a bar within the theory of <strong> linear elasticity </strong> 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.
101
102
103
104
105
106
107
108
109
110
111
112
113

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

\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}

```

114
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.
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

The weak formulation of the fundamental law of dynamics in the deformed configuration reads:
```math
\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:  

```math
\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 <strong>plane strain </strong> hypothesis, namely $` \varepsilon_{33} = 0`$ which implies: 

```math
\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:

```math
\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:

```math
\begin{cases}
DIAZ Jerome's avatar
DIAZ Jerome committed
162
163
\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 \\
164
165
166
167
168
169
170
\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
DIAZ Jerome's avatar
DIAZ Jerome committed
171
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: 
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237

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

```math
\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: 
```math
\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:
```math
\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:
```math
DIAZ Jerome's avatar
DIAZ Jerome committed
238
239
\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
240
241
242
243
244
245
```
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:

```math
DIAZ Jerome's avatar
DIAZ Jerome committed
246
247
\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]
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262

```
```math
\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:
```math
\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 <strong>mass matrix</strong>, $` \underline{\underline{\mathbb{K}}} `$  to the <strong>sitffness matrix</strong> and $` \underline{\mathbb{F}} `$ is the discretized <strong> right-hand side</strong> (corresponding to a surfacic load). Each of these 3 quantities correspond to different operators that will be assembled in the <strong>variational formulation</strong> 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}} `$:

```math
DIAZ Jerome's avatar
DIAZ Jerome committed
263
\int_{\Omega _{0}}^{} \rho_0 \, \underline{y}^{*} \cdot  \underline{\ddot{y}} \, \textrm{d}\Omega_0
264
265
266
267
268
269
```
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:
```math
\frac{\dot{\underline{y}}^{n+1} + \dot{\underline{y}}^n}{2} = \frac{\underline{y}^{n+1} - \underline{y}^n}{\Delta t}
```

270
271
272
273
274
275
276
277
278
279
280
# Creating the new project

We will create a new project in _Sources/ModelInstances_ called DemoElasticity.

## If you're not using XCode

Read 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:
281
282
283
284
285
* 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
286

287
* 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).
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303





## 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:
304
305
306
307
* 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:
308
````
309
--input_data ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/DemoElasticity/demo.lua    
310
311
312
313
314
315
316
317
318
319
320
321
322
````

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).

- 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_:
323
324
325
326
327
328
- 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).


329
330
## Model files

331
- 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: 
332

333
334
335
336
* 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.

337
338
## Variational formulation files

339
- Repeat mostly what was done for the main file, except you now choose _VariationalFormulation_ in MoReFEM section.
340
341
342
343


## InputData file

344
- Repeat mostly what was done for the main file, except now you choose _InputData_ in MoReFEM section. 
345
346
347
348
349

## 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).

350
<p><strong><font color="red">In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:</font></strong></p>
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370

````
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:

371
372
373
* 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.
374
375
376

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

377
<p><strong><font color="red">In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:</font></strong></p>
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401

````
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)
````

402
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.
403
404
405
406
407
408


````
//! \copydoc doxygen_hide_mesh_enum
enum class MeshIndex
{
409
mesh = 10 // only one mesh considered in current model; we arbitrarily index it by a 10.
410
411
412
413
414
415
};


//! \copydoc doxygen_hide_domain_enum
enum class DomainIndex
{
416
417
418
419
highest_dimension = 1,
neumann = 2,
dirichlet = 3,
full_mesh = 10
420
421
422
423
424
425
};


//! \copydoc doxygen_hide_felt_space_enum
enum class FEltSpaceIndex
{
426
427
highest_dimension = 1, // doesn't matter the name is the same or not as in DomainIndex
neumann = 2
428
429
430
431
432
433
};


//! \copydoc doxygen_hide_unknown_enum
enum class UnknownIndex
{
434
displacement = 1
435
436
437
438
439
440
};


//! \copydoc doxygen_hide_solver_enum
enum class SolverIndex
{
441
solver = 1
442
443
444
445
446
447
};


//! \copydoc doxygen_hide_numbering_subset_enum
enum class NumberingSubsetIndex
{
448
monolithic = 1
449
450
451
452
453
454
};


//! \copydoc doxygen_hide_source_enum
enum class SourceIndex
{
455
surfacic = 1
456
457
458
459
460
461
};


//! \copydoc doxygen_hide_boundary_condition_enum
enum class BoundaryConditionIndex
{
462
sole = 1
463
464
465
466
467
468
469
470
471
};

````

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
472
using InputDataTuple = std::tuple
473
<
474
InputDataNS::TimeManager,
475

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

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

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

482
483
484
485
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::highest_dimension)>,
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::neumann)>,
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::dirichlet)>,
InputDataNS::Domain<EnumUnderlyingType(DomainIndex::full_mesh)>,
486

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

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

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

494
495
496
497
InputDataNS::Solid::VolumicMass,
InputDataNS::Solid::YoungModulus,
InputDataNS::Solid::PoissonRatio,
InputDataNS::Solid::PlaneStressStrain,
498

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

501
InputDataNS::Result
502
503
504
>;
````

505
The code should compile as well, so:
506

507
<p><strong><font color="red">In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:</font></strong></p>
508
509
510
511
512
513
514
515
516
517
518
519
````
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:

520
521
522
* By pressing **Command** + **R** 
* Or through the menu: _Product_ > _Run_
* If you're using cmake, use the generated demo_elasticity executable:
523
524
525
526
527
528
529
530
531

````
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.
532
Exception caught from MoReFEMData<InputData>: Exception found at Sources/Core/MoReFEMData/MoReFEMData.hxx, line 150: Input parameter file to edit and complete!
533
534
535
536
537
538
539
540
541
542
543
````

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

544
__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).
545
546
547

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

548
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.
549
550
551
552


#### Transient block

553
This block sets the time interval over which the model is run. Choose `timeMax = 0.5`  for instance.
554
555
556

#### NumberingSubset1

557
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"`.
558
559
560

#### Unknown1

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

563
So choose "displacement" (or whatever name you want: the only constraint is that two different unknowns cannot share the same name) and "vectorial".
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582

````
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

##### Domain1

583
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 \} `$.
584
585
586
587
588
589
590
591
592
593

````
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

594
Second Domain is the area upon which the Neumann boundary condition is applied. This corresponds to $`\Gamma^N`$.
595
596
597
598
599
600
601
602
603
604

````
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

605
Third Domain is the area upon which Dirichlet boundary condition is applied. This corresponds to $` \Gamma^D`$.
606
607
608
609
610
611
612
613
614
615
616

````
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

617
Fourth domain is the whole mesh with nothing left aside. This corresponds to $`\Omega_0`$.
618
619
620
621
622
623
624
625

````
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
````

626
627
#### EssentialBoundaryCondition1. 
This corresponds to $` \underline{y}(\underline{x}) = \underline{0} \quad \text{on} \quad \Gamma^D`$
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665

````
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
````

#### 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.

````
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 } 

````

#### Solid 

666
##### Volumic mass $` \rho_0 `$ 
667
668
669
670
671
672

````
nature = "constant",
value = 1.3
````

673
##### Young modulus $`E`$
674
675
676
677

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

682
##### Poisson ratio $`\nu`$
683
684
685
686
687
688

````
nature = "constant",
value = 0.04
````

689
##### Transient source. 
DIAZ Jerome's avatar
DIAZ Jerome committed
690
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 `$ 
691
692
693
694
695
696
697
698
699
700
701


````
nature_x = "constant",
nature_y = "constant",
nature_z = "constant",
value_x = 0.,
value_y = 5.e-3,
value_z = 0.
````

702
##### Result
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729

````
output_directory = "${MOREFEM_RESULT_DIR}/DemoElasticity",
````

## Running again... and commit!

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

<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

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


## Add variational_formulation object in the Model

<p><strong><font color="green">In Model.hpp:</font></strong></p>

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

// At the bottom of the class definition
private:
730
731
732
733
734
735
736
737

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

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


738
private:
739
740
741

//! Underlying variational formulation.
VariationalFormulation::unique_ptr variational_formulation_ = nullptr;
742
743
744
745
746
747
748
749
750

````


<p><strong><font color="green">In Model.hxx:</font></strong></p>

````
inline const VariationalFormulation& Model::GetVariationalFormulation() const noexcept
{
751
752
assert(!(!variational_formulation_));
return *variational_formulation_;
753
754
755
756
757
}


inline VariationalFormulation& Model::GetNonCstVariationalFormulation() noexcept
{
758
return const_cast<VariationalFormulation&>(GetVariationalFormulation());
759
760
761
762
763
764
765
766
}
````

<p><strong><font color="green">In Model.cpp, fill SupplInitialize():</font></strong></p>

````
void Model::SupplInitialize()
{
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
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());
783
784
785
786
787
788
789
790
791
792
793
794
795
}
````


<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

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

## Variational formulation: defining the Parameters

796
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.
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817

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.

<p><strong><font color="green">In VariationalFormulation.hpp, in the class declaration:</font></strong></p>

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:

818
819
820
821
822
823
824
825
826
827
828
//! 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;
829
830
831
832


private:

833
834
835
836
837
//! Volumic mass.
ScalarParameter<>::unique_ptr volumic_mass_ = nullptr;

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

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

842
843
//! Source parameter.
typename source_parameter_type::unique_ptr source_parameter_ = nullptr;
844
845
846
847
848
849
850
851
````

<p><strong><font color="green">In VariationalFormulation.hxx:</font></strong></p>

````
inline const ScalarParameter<>& VariationalFormulation
::GetVolumicMass() const noexcept
{
852
853
assert(!(!volumic_mass_));
return *volumic_mass_;
854
855
856
857
858
859
}


inline const ScalarParameter<>& VariationalFormulation
::GetYoungModulus() const noexcept
{
860
861
assert(!(!young_modulus_));
return *young_modulus_;
862
863
864
865
866
867
}


inline const ScalarParameter<>& VariationalFormulation
::GetPoissonRatio() const noexcept
{
868
869
assert(!(!poisson_ratio_));
return *poisson_ratio_;
870
871
872
873
874
875
}


inline const VariationalFormulation::source_parameter_type& VariationalFormulation
::GetSourceParameter() const noexcept
{
876
877
assert(!(!source_parameter_));
return *source_parameter_;
878
879
880
881
882
883
884
885
886
887
888
889
}
````


<p><strong><font color="green">In VariationalFormulation.cpp, fill SupplInit():</font></strong></p>


````
#include "Parameters/InitParameterFromInputData/InitParameterFromInputData.hpp"

void VariationalFormulation::SupplInit(const input_data_type& input_data)
{
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
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);

921
922
923
924
925
926
}
````


<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

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

929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
````
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`:

<p><strong><font color="green">In VariationalFormulation.hpp:</font></strong></p>

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:

````
/*!
944
945
946
* \copydoc doxygen_hide_varf_constructor
* \param[in] numbering_subset The only \a NumberingSubset considered in this model.
*/
947
explicit VariationalFormulation(const morefem_data_type& morefem_data,
948
949
950
951
const GodOfDof& god_of_dof,
DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
const NumberingSubset& numbering_subset,
TimeManager& time_manager);
952
953
954
955
956
957
958
````

and add the data attribute and a public accessor:

````
public:

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

962
963
private:

964
965
//! Only numbering subset considered in this model.
const NumberingSubset& numbering_subset_;
966
967
968
969
970
971
972
973
974

````

<p><strong><font color="green">In VariationalFormulation.hxx:</font></strong></p>


````
inline const NumberingSubset& VariationalFormulation::GetNumberingSubset() const noexcept
{
975
return numbering_subset_;
976
977
978
979
980
981
982
}
````

<p><strong><font color="green">In VariationalFormulation.cpp, modify the constructor:</font></strong></p>

````
VariationalFormulation::VariationalFormulation(const morefem_data_type& morefem_data,
983
984
985
986
const GodOfDof& god_of_dof,
DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
const NumberingSubset& numbering_subset,
TimeManager& time_manager)
987
: parent(morefem_data,
988
989
990
time_manager,
god_of_dof,
std::move(boundary_condition_list)),
991
numbering_subset_(numbering_subset)
992
993
994
995
996
997
998
{ }
````

<p><strong><font color="green">In Model.cpp, modify the construction of the variational formulation in SupplInitialize():</font></strong></p>

````
decltype(auto) numbering_subset =
999
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::monolithic));
1000
1001

variational_formulation_ = std::make_unique<VariationalFormulation>(morefem_data,
1002
1003
1004
1005
1006
1007
god_of_dof,
std::move(bc_list),
numbering_subset,
time_manager);


1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
````

<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

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


## Add method to run the static case

1020
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}} `$ .
1021
1022
1023
1024
1025
1026

<p><strong><font color="green">In VariationalFormulation.hpp, in the class:</font></strong></p>

````
public:

1027
1028
//! Run the static case.
void RunStaticCase();
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
````
<p><strong><font color="green">In VariationalFormulation.cpp, in the class:</font></strong></p>

````
void VariationalFormulation::RunStaticCase()
{ }
````


<p><strong><font color="green">In Model.cpp:</font></strong></p>

````
void Model::SupplInitialize()
{
1043
...
1044
decltype(auto) variational_formulation = GetNonCstVariationalFormulation();
1045
1046
1047

variational_formulation.RunStaticCase();
variational_formulation.WriteSolution(time_manager, variational_formulation.GetNumberingSubset());
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
}
````

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

<p><strong><font color="green">In VariationalFormulation.cpp, fill AllocateMatricesAndVectors():</font></strong></p>

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

1062
1063
parent::AllocateSystemMatrix(numbering_subset, numbering_subset);
parent::AllocateSystemVector(numbering_subset);
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
}
````

The code should now run properly.


<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

````
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)." 
````

1077
## VariationalFormulation: define both operators ($`\underline{\underline{\mathbb{K}}}`$ and $` \underline{\mathbb{F}}`$ )  required for the static case
1078
1079
1080

### Source operator

DIAZ Jerome's avatar
DIAZ Jerome committed
1081
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 `$:
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092

<p><strong><font color="green">In VariationalFormulation.hpp:</font></strong></p>

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

In the aliases at the top of the class:

````
private:
1093
1094
1095

//! Alias to the type of the source operator.
using source_type_operator = GlobalVariationalOperatorNS::TransientSource<ParameterNS::Type::vector>;
1096
1097
1098
1099
1100
1101
1102
1103
````

And within the class:


````
private:

1104
1105
//! Surfacic source operator.
const source_type_operator& GetSurfacicSourceOperator() const noexcept;
1106
1107
1108

private:

1109
1110
//! Surfacic source operator.
source_type_operator::const_unique_ptr surfacic_source_operator_ = nullptr;
1111
1112
1113
1114
1115
1116
1117
1118
````

<p><strong><font color="green">In VariationalFormulation.hxx</font></strong></p>

````
inline const VariationalFormulation::source_type_operator& VariationalFormulation
::GetSurfacicSourceOperator() const noexcept
{
1119
1120
assert(!(!surfacic_source_operator_));
return *surfacic_source_operator_;
1121
1122
1123
1124
1125
}
````

I also usually define a method `InitializeOperators`:

1126
<p><strong><font color="green">In VariationalFormulation.hpp (parameter will be used a bit later):</font></strong></p> 
1127
1128
1129
1130
1131


````
private:

1132
1133
1134
1135
1136
1137
/*!
* \brief Initialize the operators.
*
* \copydoc doxygen_hide_input_data_arg
*/
void InitializeOperators(const input_data_type& input_data);
1138
1139
1140
1141
````

<p><strong><font color="green">In VariationalFormulation.cpp:</font></strong></p>

1142
- Add a new line at the end of SupplInit:
1143
1144
1145
1146

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

1149
InitializeOperators(input_data);
1150
1151
1152
1153
1154
1155
1156
1157
}
````

- And provide the implementation for the new method:

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

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

1164
1165
1166
{
decltype(auto) felt_space_neumann =
god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::neumann));
1167

1168
decltype(auto) source_parameter = GetSourceParameter();
1169

1170
1171
1172
1173
1174
surfacic_source_operator_ =
std::make_unique<source_type_operator>(felt_space_neumann,
displacement_ptr,
source_parameter);
}
1175
1176
1177
1178
1179
}
````

### Stiffness operator

DIAZ Jerome's avatar
DIAZ Jerome committed
1180
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 `$
1181
1182
1183
1184
1185
1186
1187
1188
1189
<p><strong><font color="green">In VariationalFormulation.hpp:</font></strong></p>

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

````
private:

1190
1191
//! Get the stiffness operator.
const GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor& GetStiffnessOperator() const noexcept;
1192
1193
1194

private:

1195
1196
1197
//! Stiffness operator.
GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor::const_unique_ptr
stiffness_operator_ = nullptr;
1198
1199
1200
1201
1202
1203
````


<p><strong><font color="green">In VariationalFormulation.hxx</font></strong></p>

````
1204
1205
1206
1207
1208
1209
inline const GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor&
VariationalFormulation::GetStiffnessOperator() const noexcept
{
assert(!(!stiffness_operator_));
return *stiffness_operator_;
}
1210
1211
1212
1213
1214
1215
1216
````

<p><strong><font color="green">In VariationalFormulation.cpp, complete InitializeOperators():</font></strong></p>

````
void VariationalFormulation::InitializeOperators(const input_data_type& input_data)
{
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
... 

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);
1235
1236
1237
1238
}
````


1239
<p><strong><font color="red">In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:</font></strong></p>
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254

````
git add Sources/ModelInstances/DemoElasticity
git commit -m "#0 Add both operators required to write 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:

<p><strong><font color="green">In VariationalFormulation.hpp:</font></strong></p>

````
private:

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

1258
1259
//! Non constant accessor to the stiffness matrix.
GlobalMatrix& GetNonCstStiffnessMatrix() noexcept;
1260
1261
1262
1263
1264
````      

```` 
private:

1265
1266
1267
1268
1269
1270
/*!
* \brief Stiffness matrix.
*
* \internal <b><tt>[internal]</tt></b> Actually we could dodge this one entirely but the code is more readable by making it explicit.
*/
GlobalMatrix::unique_ptr stiffness_matrix_ = nullptr;
1271
1272
1273
1274
1275
1276
1277
````


<p><strong><font color="green">In VariationalFormulation.hxx</font></strong></p>

````
inline const GlobalMatrix& VariationalFormulation
1278
::GetStiffnessMatrix() const noexcept
1279
{
1280
1281
assert(!(!stiffness_matrix_));
return *stiffness_matrix_;
1282
1283
1284
1285
1286
1287
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstStiffnessMatrix() noexcept
{
1288
return const_cast<GlobalMatrix&>(GetStiffnessMatrix());
1289
1290
1291
1292
1293
1294
1295
1296
1297
}
````

<p><strong><font color="green">In VariationalFormulation.cpp:</font></strong></p>

````
void VariationalFormulation
::AllocateMatricesAndVectors()
{
1298
1299
1300
1301
1302
...

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

stiffness_matrix_ = std::make_unique<GlobalMatrix>(system_matrix);
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
}
````

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:

<p><strong><font color="green">In VariationalFormulation.cpp:</font></strong></p>


````
void VariationalFormulation
1318
::SupplInit(const InputData& input_data)
1319
{
1320
1321
1322
1323
1324
1325
1326
...

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

1330
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).
1331
1332
1333

## Static case (at last!)

1334
1335
1336
1337
1338
1339
1340
We are now finally able to run the whole static case, sovling:
```math
\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}
```
1341
1342
1343
1344
1345
1346

<p><strong><font color="green">In VariationalFormulation.cpp:</font></strong></p>

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

1349
1350
1351
1352
1353
1354
// 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);
}
1355

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


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

1363
parent::DebugPrintSolutionElasticWithOneUnknown(numbering_subset, 10);
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
}
````


<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

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

## Defining the mass operator

DIAZ Jerome's avatar
DIAZ Jerome committed
1377
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  `$
1378

1379
1380
1381
1382
1383
1384
1385
1386
1387
<p><strong><font color="green">In VariationalFormulation.hpp:</font></strong></p>

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

````
private:

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

private:

1393
1394
//! Mass operator.
GlobalVariationalOperatorNS::Mass::const_unique_ptr mass_operator_ = nullptr;
1395
1396
1397
1398
1399
1400
1401
1402
````    

<p><strong><font color="green">In VariationalFormulation.hxx:</font></strong></p>

````
inline const GlobalVariationalOperatorNS::Mass& VariationalFormulation
::GetMassOperator() const noexcept
{
1403
1404
assert(!(!mass_operator_));
return *mass_operator_;
1405
1406
1407
1408
1409
1410
}
````

<p><strong><font color="green">In VariationalFormulation.cpp, complete InitializeOperators():</font></strong></p>

````
1411
void VariationalFormulation::InitializeOperators(const InputData& input_data)
1412
{
1413
1414
1415
1416
...
mass_operator_ = std::make_unique<GlobalVariationalOperatorNS::Mass>(felt_space_highest_dimension,
displacement_ptr,
displacement_ptr);
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
}
````


<p><strong><font color="red">In a terminal at the root of MoReFEM:</font></strong></p>

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


## Additional global linear algebra

1431
1432
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):
1433
1434
1435
1436
1437
1438

<p><strong><font color="green">In VariationalFormulation.hpp:</font></strong></p>

````
private:

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

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

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

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

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

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

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

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

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

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

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

1472
///@}
1473
1474
1475

private:

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

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

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

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

1488
1489
//! Mass matrix.
GlobalMatrix::unique_ptr mass_matrix_ = nullptr;
1490
1491
1492
1493
1494
1495
1496
1497
1498
`````

<p><strong><font color="green">In VariationalFormulation.hxx</font></strong></p>


````
inline const GlobalVector& VariationalFormulation
::GetVectorCurrentDisplacement() const noexcept
{
1499
1500
assert(!(!vector_current_displacement_));
return *vector_current_displacement_;
1501
1502
1503
1504
1505
1506
}


inline GlobalVector& VariationalFormulation
::GetNonCstVectorCurrentDisplacement() noexcept
{
1507
return const_cast<GlobalVector&>(GetVectorCurrentDisplacement());
1508
1509
1510
1511
1512
1513
}


inline const GlobalVector& VariationalFormulation
::GetVectorCurrentVelocity() const noexcept
{
1514
1515
assert(!(!vector_current_velocity_));
return *vector_current_velocity_;
1516
1517
1518
1519
1520
1521
}


inline GlobalVector& VariationalFormulation
::GetNonCstVectorCurrentVelocity() noexcept
{
1522
return const_cast<GlobalVector&>(GetVectorCurrentVelocity());
1523
1524
1525
1526
1527
1528
}


inline const GlobalMatrix& VariationalFormulation
::GetMatrixCurrentDisplacement() const noexcept
{
1529
1530
assert(!(!matrix_current_displacement_));
return *matrix_current_displacement_;
1531
1532
1533
1534
1535
1536
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstMatrixCurrentDisplacement() noexcept
{
1537
return const_cast<GlobalMatrix&>(GetMatrixCurrentDisplacement());
1538
1539
1540
1541
1542
1543
}


inline const GlobalMatrix& VariationalFormulation
::GetMatrixCurrentVelocity() const noexcept
{
1544
1545
assert(!(!matrix_current_velocity_));
return *matrix_current_velocity_;
1546
1547
1548
1549
1550
1551
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstMatrixCurrentVelocity() noexcept
{
1552
return const_cast<GlobalMatrix&>(GetMatrixCurrentVelocity());
1553
1554
1555
1556
1557
1558
}


inline const GlobalMatrix& VariationalFormulation
::GetMassMatrix() const noexcept
{
1559
1560
assert(!(!mass_matrix_));
return *mass_matrix_;
1561
1562
1563
1564
1565
1566
}


inline GlobalMatrix& VariationalFormulation
::GetNonCstMassMatrix() noexcept
{
1567
return const_cast<GlobalMatrix&>(GetMassMatrix());
1568
1569
1570
1571
1572
1573
1574
1575
1576
}
````

<p><strong><font color="green">In VariationalFormulation.cpp:</font></strong></p>

````
void VariationalFormulation
::AllocateMatricesAndVectors()
{
1577
1578
1579
1580
1581
1582
1583
1584
1585
...

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);
1586
1587
1588
1589
}
`````


1590
<p><strong><font color="red">In a terminal at the root of MoReFEM project (${HOME}/Codes/MoReFEM/CoreLibrary), after a check compilation works:</font></strong></p>
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603

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

## PrepareDynamicRun()

<p><strong><font color="green">In Model.cpp:</font></strong></p>

````
void Model::SupplInitialize()
{
1604
1605
...
variational_formulation.PrepareDynamicRuns();
1606
1607
1608
1609
1610