diff --git a/doc/dev-doc/manual/diffusion.rst b/doc/dev-doc/manual/diffusion.rst
index c0f86910b..705c78d18 100644
--- a/doc/dev-doc/manual/diffusion.rst
+++ b/doc/dev-doc/manual/diffusion.rst
@@ -1,177 +1,181 @@
+.. _sect-dm:
+
Diffusion Model
---------------
The heat transfer model is a specific implementation of the :cpp:class:`Model
<akantu::Model>` interface dedicated to handle the dynamic heat equation.
Theory
``````
The strong form of the dynamic heat equation
can be expressed as:
.. math::
\rho c_v \dot{T} + \nabla \cdot \vec{\kappa} \nabla T = b
with :math:`T` the scalar temperature field, :math:`c_v` the specific heat capacity, :math:`\rho`
the mass density, :math:`\mat{\kappa}` the conductivity tensor, and :math:`b` the heat
generation per unit of volume. The discretized weak form with a finite number of
elements is
.. math::
\forall i \quad
\sum_j \left( \int_\Omega \rho c_v N_j N_i d\Omega \right) \dot{T}_j
- \sum_j \left( \int_\Omega \vec{\kappa} \nabla N_j \nabla N_i d\Omega \right) T_j =
- \int_{\Gamma} N_i \vec{q} \cdot \vec{n} d\Gamma + \int_\Omega b N_i d\Omega
with :math:`i` and :math:`j` the node indices, :math:`\vec{n}` the normal field to the surface
:math:`\Gamma = \partial \Omega`.
To simplify, we can define the capacity and the conductivity matrices as
.. math::
C_{ij} = \int_\Omega \rho c_v N_j N_i d\Omega \qquad \textrm{and} \qquad
K_{ij} = - \int_\Omega \vec{\kappa} \nabla N_j \nabla N_i d\Omega
and the system to solve can be written:
.. math::
\mat{C} \cdot \vec{\dot{T}} = \vec{Q}^{\text{ext}} -\mat{K} \cdot \vec{T}~,
with :math:`\vec{Q}^{\text{ext}}` the consistent heat generated.
The diffusion model is meant as a base to implement different diffusive processes, heat diffusion, chemical diffusion, simple flow problems, etc. Currently only one of this physics is implemented, the heat diffusion, in a model named HeatTransferModel
+.. _sect-dm-using:
+
Using the Heat Transfer Model
`````````````````````````````
A material file name has to be provided during initialization.
Currently, the :cpp:class:`HeatTransferModel <akantu::HeatTransferModel>` object uses dynamic analysis
with an explicit time integration scheme. It can simply be created
like this
.. code-block:: c++
HeatTransferModel model(mesh, spatial_dimension);
while an existing mesh has been used (see \ref{sect:common:mesh}).
Then the model object can be initialized with:
.. code-block:: c++
model.initFull()
This function will load the material properties, and allocate / initialize the nodes and element :cpp:class:`Arrays <akantu::Array>`
More precisely, the heat transfer model contains 4 :cpp:class:`Arrays <akantu::Array>`:
- **temperature** contains the nodal temperature :math:`T` (zero by default after the initialization).
- **temperature_rate** contains the variations of temperature :math:`\dot{T}` (zero by default after the initialization).
- **blocked_dofs** contains a Boolean value for each degree of freedom specifying whether the degree is blocked or not. A Dirichlet boundary condition (:math:`T_d`) can be prescribed by setting the **blocked_dofs** value of a degree of freedom to ``true``. The **temperature** and the **temperature_rate** are computed for all degrees of freedom where the **blocked_dofs** value is set to ``false``. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept.
- **external_heat_rate** contains the external heat generations. :math:`\vec{Q^{ext}}` on the nodes.
- **internal_heat_rate** contains the internal heat generations. :math:`\vec{Q^{int}} = -\mat{K} \cdot \vec{T}` on the nodes.
Only a single material can be specified on the domain. A material text file (*e.g.* material.dat) provides the material properties as follows:
.. code-block:: python
model heat_transfer_model [
constitutive_law heat_diffusion [
name = %\emph{XXX}%
capacity = %\emph{XXX}%
density = %\emph{XXX}%
conductivity = [%\emph{XXX}% ... %\emph{XXX}%]
]
]
where the ``capacity`` and ``density`` are scalars, and the ``conductivity`` is
specified as a :math:`3\times 3` tensor.
Explicit Dynamic
````````````````
The explicit time integration scheme in ``Akantu`` uses a lumped capacity matrix
:math:`\mat{C}` (reducing the computational cost, see Chapter :ref:`sect-smm`).
This matrix is assembled by distributing the capacity of each element onto its
nodes. Therefore, the resulting :math:`\mat{C}` is a diagonal matrix stored in
the ``capacity`` :cpp:class:`Array <akantu::Array>` of the model.
.. code-block:: c++
model.assembleCapacityLumped();
.. note::
Currently, only the explicit time integration with lumped capacity
matrix is implemented within ``Akantu``.
The explicit integration scheme is *Forward Euler* :cite:`curnier92a`.
- Predictor: :math:`\vec{T}_{n+1} = \vec{T}_{n} + \Delta t \dot{\vec{T}}_{n}`
- Update residual: :math:`\vec{R}_{n+1} = \left( \vec{Q^{ext}_{n+1}} - \vec{K}\vec{T}_{n+1} \right)`
- Corrector : :math:`\dot{\vec{T}}_{n+1} = \mat{C}^{-1} \vec{R}_{n+1}`
The explicit integration scheme is conditionally stable. The time step has to be
smaller than the stable time step, and it can be obtained in ``Akantu`` as
follows:
.. code-block:: c++
time_step = model.getStableTimeStep();
The stable time step is defined as:
.. math::
\Delta t_{\st{crit}} = 2 \Delta x^2 \frac{\rho c_v}{\mid\mid \mat{\kappa} \mid\mid^\infty}
:label: eqn:htm:explicit:stabletime
where :math:`\Delta x` is the characteristic length (*e.g* the in-radius in the
case of linear triangle element), :math:`\rho` is the density,
:math:`\mat{\kappa}` is the conductivity tensor, and :math:`c_v` is the specific
heat capacity. It is necessary to impose a time step which is smaller than the
stable time step, for instance, by multiplying the stable time step by a safety
factor smaller than one.
.. code-block:: c++
const Real safety_time_factor = 0.1;
Real applied_time_step = time_step * safety_time_factor;
model.setTimeStep(applied_time_step);
The following loop allows, for each time step, to update the ``temperature``,
``residual`` and ``temperature_rate`` fields following the previously described
integration scheme.
.. code-block:: c++
for (Int s = 1; (s-1)*applied_time_step < total_time; ++s) {
model.solveStep();
}
An example of explicit dynamic heat propagation is presented in
``examples/heat_transfer/explicit_heat_transfer.cc``. This example consists
of a square 2D plate of :math:`1 \text{m}^2` having an initial temperature of
:math:`100 \text{K}` everywhere but a none centered hot point maintained at
:math:`300 \text{K}`. :numref:`fig:htm:explicit:dynamic-1` presents the geometry
of this case. The material used is a linear fictitious elastic material with a
density of :math:`8940 \text{kg}/\text{m}^3`, a conductivity of
:math:`401 \text{W}/\text{m}/\text{K}` and a specific heat capacity of
:math:`385 \text{J}/\text{K}/\text{kg}`. The time step used is
:math:`0.12 \text{s}`.
.. _fig:htm:explicit:dynamic-1:
.. figure:: figures/hot-point-1.png
:align: center
Initial temperature field
.. _fig:htm:explicit:dynamic-2:
.. figure:: figures/hot-point-2.png
:align: center
Temperature field after 15000 time steps = 30 minutes. The lines represent iso-surfaces.
diff --git a/examples/c++/diffusion_model/README.rst b/examples/c++/diffusion_model/README.rst
index c15e13dc6..c6d69df27 100644
--- a/examples/c++/diffusion_model/README.rst
+++ b/examples/c++/diffusion_model/README.rst
@@ -1,46 +1,52 @@
diffusion_model
'''''''''''''''
In ``diffusion_model``, examples of the ``HeatTransferModel`` are presented.
An example of a static heat propagation is presented in
``heat_diffusion_static_2d.cc``. This example consists of a square 2D plate of
:math:`1 \text{m}^2` having an initial temperature of :math:`100 \text{K}`
everywhere but a non centered hot point maintained at
:math:`300 \text{K}`. :numref:`fig-ex-diffusion_static` presents the geometry
of this case (left) and the results (right). The material used is a linear
fictitious elastic material with a density of :math:`8940 \text{kg}/\text{m}^3`,
a conductivity of :math:`401 \text{W}/\text{m}/\text{K}` and a specific heat
capacity of :math:`385 \text{J}/\text{K}/\text{kg}`.
+The simulation is set following the procedure described in :ref:`_sect-dm-using`
+
.. _fig-ex-diffusion_static:
.. figure:: examples/c++/diffusion_model/images/diffusion_static.png
:align: center
:width: 70%
Initial (left) and final (right) temperature field
In ``heat_diffusion_dynamics_2d.cc``, the same example is solved dynamically
-using an explicit time scheme. The time step used is :math:`0.12 \text{s}`.
+using an explicit time scheme. The time step used is :math:`0.12 \text{s}`. The only main difference with the previous example lies in the model initiation::
+
+ model.initFull(_analysis_method = _explicit_lumped_mass);
.. _fig-ex-diffusion_explicit:
.. figure:: examples/c++/diffusion_model/images/hot-point-2.png
:align: center
:width: 60%
Temperature field after 15000 time steps = 30 minutes. The lines represent
iso-surfaces.
In ``heat_diffusion_dynamics_3d.cc``, a 3D explicit dynamic heat propagation
problem is solved. It consists of a cube having an initial temperature of
:math:`100 \text{K}` everywhere but a centered sphere maintained at
-:math:`300 \text{K}`. :numref:`fig-ex-diffusion_3d` presents the resulting
-temperature field evolution.
+:math:`300 \text{K}`.
+The simulation is set exactly as ``heat_diffusion_dynamics_2d.cc`` except that the mesh is now a 3D mesh and that the heat source has a third coordinate and is placed at the cube center.
+
+:numref:`fig-ex-diffusion_3d` presents the resulting temperature field evolution.
.. _fig-ex-diffusion_3d:
.. figure:: examples/c++/diffusion_model/images/diffusion_3d.gif
:align: center
:width: 70%
Temperature field evolution.