diff --git a/doc/dev-doc/conf.py.in b/doc/dev-doc/conf.py.in
index fb87b1f96..8f53ac5b9 100644
--- a/doc/dev-doc/conf.py.in
+++ b/doc/dev-doc/conf.py.in
@@ -1,212 +1,225 @@
# -*- coding: utf-8 -*-
#
# Configuration file for the Sphinx documentation builder.
#
# This file does only contain a selection of the most common options. For a
# full list see the documentation:
# http://www.sphinx-doc.org/en/master/config
# -- Path setup --------------------------------------------------------------
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
import os
# import subprocess
# import sys
# sys.path.insert(0, os.path.abspath('.'))
# -- Project information -----------------------------------------------------
project = '@PROJECT_NAME@'
copyright = '@AKANTU_COPYRING@'
author = '@AKANTU_MAINTAINER@'
version = '@AKANTU_VERSION_MAJOR@.@AKANTU_VERSION_MINOR@'
# The full version, including alpha/beta/rc tags
release = '@AKANTU_VERSION@'
# -- General configuration ---------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#
# needs_sphinx = '1.0'
+# Number figures
+numfig = True
+
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = [
'sphinx.ext.autodoc',
'sphinx.ext.intersphinx',
'sphinx.ext.coverage',
'sphinx.ext.mathjax',
'sphinx.ext.ifconfig',
'sphinx.ext.viewcode',
'breathe',
]
read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True'
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# The suffix(es) of source filenames.
# You can specify multiple suffix as a list of string:
#
# source_suffix = ['.rst', '.md']
source_suffix = '.rst'
# The master toctree document.
master_doc = 'index'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#
# This is also used if you do content translation via gettext catalogs.
# Usually you set "language" from the command line for these cases.
language = None
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This pattern also affects html_static_path and html_extra_path.
exclude_patterns = ["CMakeLists.txt"]
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
primary_domain = 'cpp'
highlight_language = 'cpp'
# -- Options for HTML output -------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
if read_the_docs_build:
html_theme = 'default'
else:
html_theme = 'sphinx_rtd_theme'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#
# html_theme_options = {}
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# Custom sidebar templates, must be a dictionary that maps document names
# to template names.
#
# The default sidebars (for documents that don't match any pattern) are
# defined by theme itself. Builtin themes are using these templates by
# default: ``['localtoc.html', 'relations.html', 'sourcelink.html',
# 'searchbox.html']``.
#
# html_sidebars = {}
html_sidebars = {
'**': [
'relations.html', # needs 'show_related': True theme option to display
'searchbox.html',
]}
+# MathJax configuration
+mathjax_config = {
+ 'TeX': {
+ 'Macros': {
+ 'st': [r'\mathrm{#1}', 1],
+ 'mat': [r'\mathbf{#1}', 1],
+ },
+ },
+}
+
# -- Options for HTMLHelp output ---------------------------------------------
# Output file base name for HTML help builder.
htmlhelp_basename = 'Akantudoc'
# -- Options for LaTeX output ------------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#
# 'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#
# 'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#
'preamble': r'''\usepackage{amsmath}''',
# Latex figure (float) alignment
#
# 'figure_align': 'htbp',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title,
# author, documentclass [howto, manual, or own class]).
latex_documents = [
(master_doc, 'Akantu.tex', 'Akantu Documentation',
'Nicolas Richart', 'manual'),
]
# -- Options for manual page output ------------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
(master_doc, 'akantu', 'Akantu Documentation',
[author], 1)
]
# -- Options for Texinfo output ----------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
(master_doc, 'Akantu', 'Akantu Documentation',
author, 'Akantu', 'One line description of project.',
'Miscellaneous'),
]
# -- Options for Epub output -------------------------------------------------
# Bibliographic Dublin Core info.
epub_title = project
epub_author = author
epub_publisher = author
epub_copyright = copyright
# The unique identifier of the text. This can be a ISBN number
# or the project homepage.
#
epub_identifier = ''
# A unique identification for the text.
#
# epub_uid = ''
# A list of files that should not be packed into the epub file.
epub_exclude_files = ['search.html']
# -- Extension configuration -------------------------------------------------
if read_the_docs_build:
akantu_path = "."
else:
akantu_path = "@CMAKE_CURRENT_BINARY_DIR@"
os.makedirs("@CMAKE_CURRENT_BINARY_DIR@/_static", exist_ok=True)
# print("akantu_path = '{}'".format(akantu_path))
breathe_projects = {"Akantu": os.path.join(akantu_path, "xml")}
breathe_default_project = "Akantu"
# -- Options for intersphinx extension ---------------------------------------
intersphinx_mapping = {
'numpy': ('https://docs.scipy.org/doc/numpy/', None),
'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
}
diff --git a/doc/dev-doc/index.rst b/doc/dev-doc/index.rst
index 9856890aa..2d0b0164d 100644
--- a/doc/dev-doc/index.rst
+++ b/doc/dev-doc/index.rst
@@ -1,28 +1,29 @@
.. Akantu documentation master file, created by
sphinx-quickstart on Fri Apr 17 16:35:46 2020.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
Akantu: a FEM library
=====================
.. toctree::
:maxdepth: 2
:caption: User Manual
./manual/getting_started.rst
./manual/fe_engine.rst
+ ./manual/solidmechanicsmodel.rst
.. toctree::
:maxdepth: 2
:caption: API Reference
./reference.rst
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`
diff --git a/doc/dev-doc/manual/figures/bc_and_ic_example.svg b/doc/dev-doc/manual/figures/bc_and_ic_example.svg
new file mode 100644
index 000000000..a9ae246da
--- /dev/null
+++ b/doc/dev-doc/manual/figures/bc_and_ic_example.svg
@@ -0,0 +1,1014 @@
+
+
+
+
diff --git a/doc/dev-doc/manual/figures/dirichlet.svg b/doc/dev-doc/manual/figures/dirichlet.svg
new file mode 100644
index 000000000..3d7782d4c
--- /dev/null
+++ b/doc/dev-doc/manual/figures/dirichlet.svg
@@ -0,0 +1,335 @@
+
+
+
+
diff --git a/doc/dev-doc/manual/figures/dynamic_analysis.png b/doc/dev-doc/manual/figures/dynamic_analysis.png
new file mode 100644
index 000000000..0e058f70d
Binary files /dev/null and b/doc/dev-doc/manual/figures/dynamic_analysis.png differ
diff --git a/doc/dev-doc/manual/figures/explicit.svg b/doc/dev-doc/manual/figures/explicit.svg
new file mode 100644
index 000000000..bac613bf4
--- /dev/null
+++ b/doc/dev-doc/manual/figures/explicit.svg
@@ -0,0 +1,2330 @@
+
+
diff --git a/doc/dev-doc/manual/figures/implicit_dynamic.svg b/doc/dev-doc/manual/figures/implicit_dynamic.svg
new file mode 100644
index 000000000..2404b812e
--- /dev/null
+++ b/doc/dev-doc/manual/figures/implicit_dynamic.svg
@@ -0,0 +1,381 @@
+
+
+
+
diff --git a/doc/dev-doc/manual/figures/problem_domain.svg b/doc/dev-doc/manual/figures/problem_domain.svg
new file mode 100644
index 000000000..c87491d98
--- /dev/null
+++ b/doc/dev-doc/manual/figures/problem_domain.svg
@@ -0,0 +1,423 @@
+
+
+
+
\ No newline at end of file
diff --git a/doc/dev-doc/manual/figures/static.svg b/doc/dev-doc/manual/figures/static.svg
new file mode 100644
index 000000000..39c08a593
--- /dev/null
+++ b/doc/dev-doc/manual/figures/static.svg
@@ -0,0 +1,658 @@
+
+
+
+
diff --git a/doc/dev-doc/manual/figures/static_analysis.png b/doc/dev-doc/manual/figures/static_analysis.png
new file mode 100644
index 000000000..b4cc12be3
Binary files /dev/null and b/doc/dev-doc/manual/figures/static_analysis.png differ
diff --git a/doc/dev-doc/manual/figures/tangent.svg b/doc/dev-doc/manual/figures/tangent.svg
new file mode 100644
index 000000000..6e6710e1f
--- /dev/null
+++ b/doc/dev-doc/manual/figures/tangent.svg
@@ -0,0 +1,185 @@
+
+
+
+
diff --git a/doc/dev-doc/manual/solidmechanicsmodel.rst b/doc/dev-doc/manual/solidmechanicsmodel.rst
new file mode 100644
index 000000000..82f44b2ed
--- /dev/null
+++ b/doc/dev-doc/manual/solidmechanicsmodel.rst
@@ -0,0 +1,1218 @@
+Solid Mechanics Model
+=====================
+
+The solid mechanics model is a specific implementation of the ``Model``
+interface dedicated to handle the equations of motion or equations of
+equilibrium. The model is created for a given mesh. It will create its own
+``FEEngine`` object to compute the interpolation, gradient, integration and
+assembly operations. A :cpp:class:`SolidMechanicsModel
+` object can simply be created like this::
+
+ SolidMechanicsModel model(mesh);
+
+where ``mesh`` is the mesh for which the equations are to be
+solved. A second parameter called ``spatial_dimension`` can be
+added after ``mesh`` if the spatial dimension of the problem is
+different than that of the mesh.
+
+This model contains at least the following six ``Arrays``:
+
+blocked_dofs
+ contains a Boolean value for each degree of freedom specifying whether that
+ degree is blocked or not. A Dirichlet boundary condition can be prescribed
+ by setting the **blocked_dofs** value of a degree of freedom to
+ ``true``. A Neumann boundary condition can be applied by setting the
+ **blocked_dofs** value of a degree of freedom to ``false``. The
+ **displacement**, **velocity** and **acceleration** are
+ computed for all degrees of freedom for which the **blocked_dofs**
+ value is set to ``false``. For the remaining degrees of freedom, the imposed
+ values (zero by default after initialization) are kept.
+
+displacement
+ contains the displacements of all degrees of freedom. It can be either a
+ computed displacement for free degrees of freedom or an imposed displacement
+ in case of blocked ones (:math:`\vec{u}` in the following).
+
+velocity
+ contains the velocities of all degrees of freedom. As **displacement**,
+ it contains computed or imposed velocities depending on the nature of the
+ degrees of freedom (:math:`\dot{\vec{u}}` in the following).
+
+acceleration
+ contains the accelerations of all degrees of freedom. As **displacement**,
+ it contains computed or imposed accelerations depending on the nature of the
+ degrees of freedom (:math:`\ddot{\vec{u}}` in the following).
+
+external_force
+ contains the external forces applied on the nodes
+ (:math:`\vec{f}_{\st{ext}}` in the following).
+
+internal_force
+ contains the internal forces on the nodes (:math:`\vec{f}_{\mathrm{int}}` in
+ the following).
+
+
+Some examples to help to understand how to use this model will be
+presented in the next sections.
+
+
+Model Setup
+-----------
+
+
+Setting Initial Conditions
+``````````````````````````
+
+For a unique solution of the equations of motion, initial
+displacements and velocities for all degrees of freedom must be
+specified:
+
+.. math::
+ \vec{u}(t=0) & = \vec{u}_0\\
+ \dot{\vec u}(t=0) & = \vec{v}_0
+
+The solid mechanics model can be initialized as
+follows::
+
+ model.initFull()
+
+This function initializes the internal arrays and sets them to
+zero. Initial displacements and velocities that are not equal to zero
+can be prescribed by running a loop over the total number of
+nodes. Here, the initial displacement in :math:`x`-direction and the
+initial velocity in :math:`y`-direction for all nodes is set to :math:`0.1` and :math:`1`,
+respectively::
+
+ auto & disp = model.getDisplacement();
+ auto & velo = model.getVelocity();
+
+ for (UInt node = 0; node < mesh.getNbNodes(); ++node) {
+ disp(node, 0) = 0.1;
+ velo(node, 1) = 1.;
+ }
+
+Setting Boundary Conditions
+```````````````````````````
+
+This section explains how to impose Dirichlet or Neumann boundary
+conditions. A Dirichlet boundary condition specifies the values that
+the displacement needs to take for every point :math:`x` at the boundary
+(:math:`\Gamma_u`) of the problem domain (:numref:`fig:smm:boundaries`):
+
+.. math::
+ \vec{u} = \bar{\vec u} \quad \forall \vec{x}\in \Gamma_{u}
+
+
+A Neumann boundary condition imposes the value of the gradient of the
+solution at the boundary :math:`\Gamma_t` of the problem domain
+(:numref:`fig:smm:boundaries`):
+
+.. math::
+ \vec{t} = \mat{\sigma} \vec{n} = \bar{\vec t} \quad
+ \forall \vec{x}\in \Gamma_{t}
+
+.. _fig:smm:boundaries:
+.. figure:: figures/problem_domain.svg
+ :align: center
+ :width: 75%
+
+ Problem domain :math:`\Omega` with boundary in three dimensions. The
+ Dirchelet and the Neumann regions of the boundary are denoted with
+ :math:`\Gamma_u` and :math:`\Gamma_t`, respecitvely.
+
+Different ways of imposing these boundary conditions exist. A basic
+way is to loop over nodes or elements at the boundary and apply local
+values. A more advanced method consists of using the notion of the
+boundary of the mesh. In the following both ways are presented.
+
+Starting with the basic approach, as mentioned, the Dirichlet boundary
+conditions can be applied by looping over the nodes and assigning the
+required values. :numref:`fig:smm:dirichlet_bc` shows a beam with a
+fixed support on the left side. On the right end of the beam, a load
+is applied. At the fixed support, the displacement has a given
+value. For this example, the displacements in both the :math:`x` and the
+:math:`y`-direction are set to zero. Implementing this displacement boundary
+condition is similar to the implementation of initial displacement
+conditions described above. However, in order to impose a displacement
+boundary condition for all time steps, the corresponding nodes need to
+be marked as boundary nodes using the function ``blocked``. While,
+in order to impose a load on the right side, the nodes are not marked.
+The detail codes are shown as follows::
+
+ auto & blocked = model.getBlockedDOFs();
+ const auto & pos = mesh.getNodes();
+
+ UInt nb_nodes = mesh.getNbNodes();
+
+ for (UInt node = 0; node < nb_nodes; ++node) {
+ if(Math::are_float_equal(pos(node, _x), 0)) {
+ blocked(node, _x) = true; // block dof in x-direction
+ blocked(node, _y) = true; // block dof in y-direction
+ disp(node, _x) = 0.; // fixed displacement in x-direction
+ disp(node, _y) = 0.; // fixed displacement in y-direction
+ } else if (Math::are_float_equal(pos(node, _y), 0)) {
+ blocked(node, _x) = false; // unblock dof in x-direction
+ forces(node, _x) = 10.; // force in x-direction
+ }
+ }
+
+
+.. _fig:smm:dirichlet_bc:
+.. figure:: figures/dirichlet.svg
+ :align: center
+
+ Beam with fixed support and load.
+
+
+For the more advanced approach, one needs the notion of a boundary in
+the mesh. Therefore, the boundary should be created before boundary
+condition functors can be applied. Generally the boundary can be
+specified from the mesh file or the geometry. For the first case, the
+function ``createGroupsFromMeshData`` is called. This function
+can read any types of mesh data which are provided in the mesh
+file. If the mesh file is created with Gmsh, the function takes one
+input strings which is either ``tag_0``, ``tag_1`` or
+``physical_names``. The first two tags are assigned by Gmsh to
+each element which shows the physical group that they belong to. In
+Gmsh, it is also possible to consider strings for different groups of
+elements. These elements can be separated by giving a string
+``physical_names`` to the function
+``createGroupsFromMeshData``::
+
+ mesh.createGroupsFromMeshData("physical_names").
+
+Boundary conditions support can also be created from the geometry by calling
+``createBoundaryGroupFromGeometry``. This function gathers all the elements on
+the boundary of the geometry.
+
+To apply the required boundary conditions, the function ``applyBC`` needs to be
+called on a :cpp:class:`SolidMechanicsModel `. This
+function gets a Dirichlet or Neumann functor and a string which specifies the
+desired boundary on which the boundary conditions is to be applied. The functors
+specify the type of conditions to apply. Three built-in functors for Dirichlet
+exist: ``FlagOnly, FixedValue,`` and ``IncrementValue``. The functor
+``FlagOnly`` is used if a point is fixed in a given direction. Therefore, the
+input parameter to this functor is only the fixed direction. The ``FixedValue``
+functor is used when a displacement value is applied in a fixed direction. The
+``IncrementValue`` applies an increment to the displacement in a given
+direction. The following code shows the utilization of three functors for the
+top, bottom and side surface of the mesh which were already defined in the Gmsh
+file::
+
+ model.applyBC(BC::Dirichlet::FixedValue(13.0, _y), "Top");
+ model.applyBC(BC::Dirichlet::FlagOnly(_x), "Bottom");
+ model.applyBC(BC::Dirichlet::IncrementValue(13.0, _x), "Side");
+
+To apply a Neumann boundary condition, the applied traction or stress should be
+specified before. In case of specifying the traction on the surface, the functor
+``FromTraction`` of Neumann boundary conditions is called. Otherwise, the
+functor ``FromStress`` should be called which gets the stress tensor as an input
+parameter::
+
+ Vector surface_traction = {0., 0., 1.};
+ auto surface_stress(3, 3) = Matrix::eye(3);
+
+ model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom");
+ model.applyBC(BC::Neumann::FromStress(surface_stress), "Top");
+
+If the boundary conditions need to be removed during the simulation, a
+functor is called from the Neumann boundary condition to free those
+boundary conditions from the desired boundary::
+
+ model.applyBC(BC::Neumann::FreeBoundary(), "Side");
+
+User specified functors can also be implemented. A full example for
+setting both initial and boundary conditions can be found in
+``examples/boundary_conditions.cc``. The problem solved
+in this example is shown in Fig.~\ref{fig:smm:bc_and_ic}. It consists
+of a plate that is fixed with movable supports on the left and bottom
+side. On the right side, a traction, which increases linearly with the
+number of time steps, is applied. The initial displacement and
+velocity in :math:`x`-direction at all free nodes is zero and two
+respectively.
+
+.. _fig:smm:bc_and_ic:
+.. figure:: figures/bc_and_ic_example.svg
+ :align: center
+ :width: 75%
+
+ Plate on movable supports.
+
+..
+ \begin{figure}[!htb]
+ \centering
+ \includegraphics[scale=0.8]{figures/bc_and_ic_example}
+ \caption{Plate on movable supports.\label{fig:smm:bc_and_ic}}
+ \end{figure}
+
+As it is mentioned in Section \ref{sect:common:groups}, node and
+element groups can be used to assign the boundary conditions. A
+generic example is given below with a Dirichlet boundary condition::
+
+ // create a node group
+ NodeGroup & node_group = mesh.createNodeGroup("nodes_fix");
+
+ /* fill the node group with the nodes you want */
+
+ // create an element group using the existing node group
+ mesh.createElementGroupFromNodeGroup("el_fix",
+ "nodes_fix",
+ spatial_dimension-1);
+
+ // boundary condition can be applied using the element group name
+ model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "el_fix");
+
+Material Selector
+`````````````````
+
+If the user wants to assign different materials to different
+finite elements groups in \akantu, a material selector has to be
+used. By default, \akantu assigns the first valid material in the
+material file to all elements present in the model (regular continuum
+materials are assigned to the regular elements and cohesive materials
+are assigned to cohesive elements or element facets).
+
+To assign different materials to specific elements, mesh data
+information such as tag information or specified physical names can be
+used. ``MeshDataMaterialSelector`` class uses this information to
+assign different materials. With the proper physical name or tag name
+and index, different materials can be assigned as demonstrated in the
+examples below::
+
+ auto mat_selector =
+ std::make_shared>("physical_names",
+ model);
+ model.setMaterialSelector(mat_selector);
+
+In this example the physical names specified in a GMSH geometry file will by
+used to match the material names in the input file.
+
+Another example would be to use the first (``tag_0``) or the second
+(``tag_1``) tag associated to each elements in the mesh::
+
+ auto mat_selector = std::make_shared>(
+ "tag_1", model, first_index);
+ model.setMaterialSelector(*mat_selector);
+
+where ``first_index`` (default is 1) is the value of ``tag_1`` that will
+be associated to the first material in the material input file. The following
+values of the tag will be associated with the following materials.
+
+There are four different material selectors pre-defined in
+\akantu. ``MaterialSelector`` and ``DefaultMaterialSelector`` is
+used to assign a material to regular elements by default. For the
+regular elements, as in the example above,
+``MeshDataMaterialSelector`` can be used to assign different
+materials to different elements.
+
+Apart from the \akantu's default material selectors, users can always
+develop their own classes in the main code to tackle various
+multi-material assignment situations.
+
+Static Analysis
+---------------
+
+The :cpp:class:`SolidMechanicsModel ` class can
+handle different analysis methods, the first one being presented is the static
+case. In this case, the equation to solve is
+
+.. math::
+ \mat{K} \vec{u} = \vec{f}_{\mathrm{ext}}
+
+where :math:`\mat{K}` is the global stiffness matrix, :math:`\vec{u}` the
+displacement vector and :math:`\vec{f}_{\st{ext}}` the vector of external
+forces applied to the system.
+
+To solve such a problem, the static solver of the
+:cpp:class:`SolidMechanicsModel ` object is used.
+First, a model has to be created and initialized. To create the model, a mesh
+(which can be read from a file) is needed, as explained in
+Section~\ref{sect:common:mesh}. Once an instance of a
+:cpp:class:`SolidMechanicsModel ` is obtained, the
+easiest way to initialize it is to use the ``initFull`` method by giving the
+``SolidMechanicsModelOptions``. These options specify the type of analysis to be
+performed and whether the materials should be initialized with ``initMaterials``
+or not::
+
+ SolidMechanicsModel model(mesh);
+ model.initFull(_analysis_method = _static);
+
+Here, a static analysis is chosen by passing the argument
+``_static`` to the method. By default, the Boolean for no
+initialization of the materials is set to false, so that they are
+initialized during the ``initFull``. The method ``initFull``
+also initializes all appropriate vectors to zero. Once the model is
+created and initialized, the boundary conditions can be set as
+explained in Section~\ref{sect:smm:boundary}. Boundary conditions
+will prescribe the external forces for some free degrees of freedom
+:math:`\vec{f}_{\st{ext}}` and displacements for some others. At this point
+of the analysis, the function
+``solveStep``\index{SolidMechanicsModel!solveStep} can be called::
+
+ model.solveStep<_scm_newton_raphson_tangent_modified,
+ SolveConvergenceCriteria::_residual>(1e-4, 1);
+
+This function is templated by the solving method and the convergence criterion
+and takes two arguments: the tolerance and the maximum number of iterations (100
+by default), which are :math:`10^{-4}` and :math:`1` for this example. The
+modified Newton-Raphson method is chosen to solve the system. In this method,
+the equilibrium equation (\ref{eqn:smm:static}) is modified in order to apply a
+Newton-Raphson convergence algorithm:
+
+.. math::
+ \mat{K}^{i+1}\delta\vec{u}^{i+1} &= \vec{r} \\
+ &= \vec{f}_{\st{ext}} -\vec{f}_{\st{int}}\\
+ &= \vec{f}_{\st{ext}} - \mat{K}^{i} \vec{u}^{i}\\
+ \vec{u}^{i+1} &= \vec{u}^{i} + \delta\vec{u}^{i+1}~,
+
+where :math:`\delta\vec{u}` is the increment of displacement to be added from
+one iteration to the other, and :math:`i` is the Newton-Raphson iteration
+counter. By invoking the ``solveStep`` method in the first step, the global
+stiffness matrix :math:`\mat{K}` from Equation~(\ref{eqn:smm:static}) is
+automatically assembled. A Newton-Raphson iteration is subsequently started,
+:math:`\mat{K}` is updated according to the displacement computed at the
+previous iteration and one loops until the forces are balanced
+(``SolveConvergenceCriteria::_residual``), i.e. :math:`||\vec{r}|| <
+\texttt{SolveConvergenceCriteria::_residual}`. One can also iterate until the
+increment of displacement is zero (``SolveConvergenceCriteria::_increment``)
+which also means that the equilibrium is found. For a linear elastic problem,
+the solution is obtained in one iteration and therefore the maximum number of
+iterations can be set to one. But for a non-linear case, one needs to iterate as
+long as the norm of the residual exceeds the tolerance threshold and therefore
+the maximum number of iterations has to be higher, e.g. :math:`100`::
+
+ model.solveStep<_scm_newton_raphson_tangent_modified,
+ SolveConvergenceCriteria::_residual>(1e-4, 100)
+
+At the end of the analysis, the final solution is stored in the
+**displacement** vector. A full example of how to solve a static
+problem is presented in the code ``examples/static/static.cc``.
+This example is composed of a 2D plate of steel, blocked with rollers
+on the left and bottom sides as shown in :numref:`fig:smm:static`.
+The nodes from the right side of the sample are displaced by :math:`0.01\%`
+of the length of the plate.
+
+.. _fig:smm:static:
+.. figure:: figures/static.svg
+ :align: center
+ :width: 75%
+
+ Numerical setup.
+
+The results of this analysis is depicted in
+:numref:`fig:smm:implicit:static_solution`.
+
+..
+ \begin{figure}[!htb]
+ \centering
+ \includegraphics[width=.7\linewidth]{figures/static_analysis}
+ \caption{Solution of the static analysis. Left: the initial
+ condition, right: the solution (deformation magnified 50 times)}
+ \label{fig:smm:implicit:static_solution}
+ \end{figure}
+
+.. _fig:smm:implicit:static_solution:
+.. figure:: figures/static_analysis.png
+ :align: center
+ :width: 75%
+
+ Solution of the static analysis. Left: the initial condition, right:
+ the solution (deformation magnified 50 times).
+
+Static implicit analysis with dynamic insertion of cohesive elements
+````````````````````````````````````````````````````````````````````
+
+
+In order to solve problems with the extrinsic cohesive method in the
+static implicit solution scheme, the function ``solveStepCohesive``
+has to be used::
+
+ model.solveStepCohesive<_scm_newton_raphson_tangent,
+ SolveConvergenceCriteria::_increment>(1e-13, error,
+ 25, false, 1e5,
+ true);
+
+
+in which the arguments are: tolerance, error, max_iteration,
+load_reduction, tol_increase_factor, do_not_factorize. This
+function, first applies the Newton-Raphson procedure to solve the
+problem. Then, it calls the method ``checkCohesiveStress`` to
+check if cohesive elements have to be inserted. Since the approach is
+implicit, only one element is added, the most stressed one (see
+Section \ref{extrinsic_insertion}). After insertion, the
+Newton-Raphson procedure is applied again to solve the same
+incremental loading step, with the new inserted cohesive element. The
+procedure loops in this way since no new cohesive elements have to be
+inserted. At that point, the solution is saved, and the simulation
+can advance to the next incremental loading step. In case the
+convergence is not reached, the obtained solution is not saved and the
+simulation return to the main file with the error given by the
+solution saved in the argument of the function *error*. In this
+way, the user can intervene in the simulation in order to find anyhow
+convergence. A possibility is, for instance, to reduce the last
+incremental loading step. The variable *load_reduction* can be
+used to identify if the load has been already reduced or not. At the
+same time, with the variable *tol_increase_factor* it is
+possible to increase the tolerance by a factor defined by the user in
+the main file, in order to accept a solution even with an error bigger
+than the tolerance set at the beginning. It is possible to increase
+the tolerance only in the phase of loading reduction, i.e., when
+load_reduction = true. A not converged solution is never saved. In
+case the convergence is not reached even after the loading reduction
+procedure, the displacement field is not updated and remains the one
+of the last converged incremental steps. Also, cohesive elements are
+inserted only if convergence is reached. An example of the extrinsic
+cohesive method in the static implicit solution scheme is presented in
+``examples/cohesive_element/cohesive_extrinsic_implicit``.
+
+Dynamic Methods
+---------------
+
+Different ways to solve the equations of motion are implemented in the
+solid mechanics model. The complete equations that should be solved
+are:
+
+.. _eqn:equation-motion:
+.. math::
+ \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} =
+ \vec{f}_{\mathrm{ext}}
+
+where :math:`\mat{M}`, :math:`\mat{C}` and :math:`\mat{K}` are the mass,
+damping and stiffness matrices, respectively.
+
+In the previous section, it has already been discussed how to solve this
+equation in the static case, where :math:`\ddot{\vec{u}} = \dot{\vec{u}} = 0`. Here
+the method to solve this equation in the general case will be presented. For
+this purpose, a time discretization has to be specified. The most common
+discretization method in solid mechanics is the Newmark-:math:`\beta` method, which is
+also the default in Akantu.
+
+For the Newmark-:math:`\beta` method, (:numref:`eqn:equation-motion`) becomes a
+system of three equations (see \cite{curnier92a} \cite{hughes-83a} for
+more details):
+
+.. _eqn:finite-difference:
+.. math::
+ \mat{M} \ddot{\vec{u}}_{n+1} + \mat{C}\dot{\vec{u}}_{n+1} + \mat{K} \vec{u}_{n+1} &={\vec{f}_{\st{ext}}}_{\, n+1}
+ \label{eqn:equation-motion-discret} \\
+ \vec{u}_{n+1} &=\vec{u}_{n} + \left(1 - \alpha\right) \Delta t \dot{\vec{u}}_{n} +
+ \alpha \Delta t \dot{\vec{u}}_{n+1} + \left(\frac{1}{2} -
+ \alpha\right) \Delta t^2
+ \ddot{\vec{u}}_{n} \label{eqn:finite-difference-1}\\
+ \dot{\vec{u}}_{n+1} &= \dot{\vec{u}}_{n} + \left(1 - \beta\right)
+ \Delta t \ddot{\vec{u}}_{n} + \beta \Delta t
+ \ddot{\vec{u}}_{n+1}
+
+
+In these new equations, :math:`\ddot{\vec{u}}_{n}`, :math:`\dot{\vec{u}}_{n}`
+and :math:`\vec{u}_{n}` are the approximations of :math:`\ddot{\vec{u}}(t_n)`,
+:math:`\dot{\vec{u}}(t_n)` and :math:`\vec{u}(t_n)`.
+Equation~(\ref{eqn:equation-motion-discret}) is the equation of motion
+discretized in space (finite-element discretization), and the equations above
+are discretized in both space and time (Newmark discretization). The
+:math:`\alpha` and :math:`\beta` parameters determine the stability and the
+accuracy of the algorithm. Classical values for :math:`\alpha` and :math:`\beta`
+are usually :math:`\beta = 1/2` for no numerical damping and :math:`0 < \alpha <
+1/2`.
+
++-------------------+-------------+----------+
+|:math:`\alpha` |Method |Type |
+| |(:math:`\beta| |
+| |= 1/2`) | |
++-------------------+-------------+----------+
+|:math:`0` |central |explicit |
+| |difference | |
++-------------------+-------------+----------+
+|:math:`\frac{1}{6}`|Fox-Goodwin |implicit |
+| |(royal road) | |
++-------------------+-------------+----------+
+|:math:`\frac{1}{3}`|Linear |implicit |
+| |acceleration | |
++-------------------+-------------+----------+
+|:math:`\frac{1}{2}`|Average |implicit |
+| |acceleration | |
+| |(trapeziodal | |
+| |rule) | |
++-------------------+-------------+----------+
+
+
+The solution of this system of equations,
+(\ref{eqn:equation-motion-discret})-(\ref{eqn:finite-difference-2}) is
+split into a predictor and a corrector system of equations. Moreover,
+in the case of a non-linear equations, an iterative algorithm such as
+the Newton-Raphson method is applied. The system of equations can be
+written as:
+
+- *Predictor:*
+.. math::
+ \vec{u}_{n+1}^{0} &= \vec{u}_{n} + \Delta t
+ \dot{\vec{u}}_{n} + \frac{\Delta t^2}{2} \ddot{\vec{u}}_{n} \\
+ \dot{\vec{u}}_{n+1}^{0} &= \dot{\vec{u}}_{n} + \Delta t
+ \ddot{\vec{u}}_{n} \\
+ \ddot{\vec{u}}_{n+1}^{0} &= \ddot{\vec{u}}_{n}
+- *Solve:*
+.. math::
+
+ \left(c \mat{M} + d \mat{C} + e \mat{K}_{n+1}^i\right)
+ \vec{w} = {\vec{f}_{\st{ext}}}_{\,n+1} - {\vec{f}_{\st{int}}}_{\,n+1}^i -
+ \mat{C} \dot{\vec{u}}_{n+1}^i - \mat{M} \ddot{\vec{u}}_{n+1}^i = \vec{r}_{n+1}^i
+- *Corrector:*
+.. math::
+ \ddot{\vec{u}}_{n+1}^{i+1} &= \ddot{\vec{u}}_{n+1}^{i} +c \vec{w} \\
+ \dot{\vec{u}}_{n+1}^{i+1} &= \dot{\vec{u}}_{n+1}^{i} + d\vec{w} \\
+ \vec{u}_{n+1}^{i+1} &= \vec{u}_{n+1}^{i} + e \vec{w}
+
+where :math:`i` is the Newton-Raphson iteration counter and :math:`c`, :math:`d` and :math:`e`
+are parameters depending on the method used to solve the equations
+
++--------------------+----------------------------+------------------------+----------------------+----------------------+
+| |:math:`\vec{w}` |:math:`e` |:math:`d` |:math:`c` |
+| | | | | |
+| | | | | |
++--------------------+----------------------------+------------------------+----------------------+----------------------+
+|in acceleration |:math:`\delta\ddot{\vec{u}}`|:math:`\alpha\beta\Delta|:math:`\beta\Delta |:math:`1` |
+| | |t^2` |t` | |
+| | | | | |
++--------------------+----------------------------+------------------------+----------------------+----------------------+
+|in velocity |:math:`\delta\dot{\vec{u}}` |:math:`\alpha\Delta t` |:math:`1` |:math:`\frac{1}{\beta |
+| | | | |\Delta t}` |
+| | | | | |
++--------------------+----------------------------+------------------------+----------------------+----------------------+
+|in displacement |:math:`\delta\vec{u}` |:math:`1` |:math:`\frac{1}{\alpha|:math:`\frac{1}{\alpha|
+| | | |\Delta t}` |\beta \Delta t^2}` |
+| | | | | |
++--------------------+----------------------------+------------------------+----------------------+----------------------+
+
+
+..
+ \begin{center}
+ \begin{tabular}{lcccc}
+ \toprule
+ & :math:`\vec{w}` & :math:`e` & :math:`d` & :math:`c`\\
+ \midrule
+ in acceleration &:math:` \delta\ddot{\vec{u}}` & :math:`\alpha \beta\Delta t^2` &:math:`\beta \Delta t` &:math:`1`\\
+ in velocity & :math:` \delta\dot{\vec{u}}`& :math:`\alpha\Delta t` & :math:`1` & :math:`\frac{1}{\beta \Delta t}`\\
+ in displacement &:math:`\delta\vec{u}` & :math:` 1` & :math:`\frac{1}{\alpha \Delta t}` & :math:`\frac{1}{\alpha \beta \Delta t^2}`\\
+ \bottomrule
+ \end{tabular}
+ \end{center}
+
+..
+ % \note{If you want to use the implicit solver \akantu should be compiled at
+ % least with one sparse matrix solver such as Mumps\cite{mumps}.}
+
+
+Implicit Time Integration
+`````````````````````````
+
+To solve a problem with an implicit time integration scheme, first a
+:cpp:class:`SolidMechanicsModel ` object has to be
+created and initialized. Then the initial and boundary conditions have to be
+set. Everything is similar to the example in the static case
+(Section~\ref{sect:smm:static}), however, in this case the implicit dynamic
+scheme is selected at the initialization of the model::
+
+ SolidMechanicsModel model(mesh);
+ model.initFull(_analysis_method = _implicit_dynamic);
+
+Because a dynamic simulation is conducted, an integration time step
+:math:`\Delta t` has to be specified. In the case of implicit simulations,
+\akantu implements a trapezoidal rule by default. That is to say
+:math:`\alpha = 1/2` and :math:`\beta = 1/2` which is unconditionally
+stable. Therefore the value of the time step can be chosen arbitrarily
+within reason::
+
+ model.setTimeStep(time_step);
+
+Since the system has to be solved for a given amount of time steps, the
+method ``solveStep()``, (which has already been used in the static
+example in Section~\ref{sect:smm:static}), is called inside a time
+loop::
+
+ /// time loop
+ Real time = 0.;
+
+ auto & solver = model.getNonLinearSolver();
+ solver.set("max_iterations", 100);
+ solver.set("threshold", 1e-12);
+ solver.set("convergence_type", SolveConvergenceCriteria::_solution);
+
+ for (UInt s = 1; time ` class. In the initialization, the explicit
+scheme is selected using the ``_explicit_lumped_mass`` constant::
+
+ SolidMechanicsModel model(mesh);
+ model.initFull(_analysis_method = _explicit_lumped_mass);
+
+
+.. note::
+ Writing ``model.initFull()`` or ``model.initFull();`` is
+ equivalent to use the ``_explicit_lumped_mass`` keyword, as this
+ is the default case.
+
+The explicit time integration scheme implemented in \akantu uses a
+lumped mass matrix :math:`\mat{M}` (reducing the computational cost). This
+matrix is assembled by distributing the mass of each element onto its
+nodes. The resulting :math:`\mat{M}` is therefore a diagonal matrix stored
+in the **mass** vector of the model.
+
+The explicit integration scheme is conditionally stable. The time step
+has to be smaller than the stable time step which is obtained in
+Akantu as follows::
+
+ critical_time_step = model.getStableTimeStep();
+
+The stable time step corresponds to the time the fastest wave (the compressive
+wave) needs to travel the characteristic length of the mesh:
+
+.. math::
+ \Delta t_{\st{crit}} = \frac{\Delta x}{c}
+
+where :math:`\Delta x` is a characteristic length (\eg the inradius in the case
+of linear triangle element) and :math:`c` is the celerity of the fastest wave in
+the material. It is generally the compressive wave of celerity :math:`c =
+\sqrt{\frac{2 \mu + \lambda}{\rho}}`, :math:`\mu` and :math:`\lambda` are the
+first and second Lame's coefficients and :math:`\rho` is the density. However,
+it is recommended to impose a time step that is smaller than the stable time
+step, for instance, by multiplying the stable time step by a safety factor
+smaller than one::
+
+ const Real safety_time_factor = 0.8;
+ Real applied_time_step = critical_time_step * safety_time_factor;
+ model.setTimeStep(applied_time_step);
+
+The initial displacement and velocity fields are, by default, equal to zero if
+not given specifically by the user (see \ref{sect:smm:initial_condition}).
+
+Like in implicit dynamics, a time loop is used in which the
+displacement, velocity and acceleration fields are updated at each
+time step. The values of these fields are obtained from the
+Newmark:math:`-\beta` equations with :math:`\beta=1/2` and :math:`\alpha=0`. In \akantu
+these computations at each time step are invoked by calling the
+function ``solveStep``::
+
+ for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) {
+ model.solveStep();
+ }
+
+The method ``solveStep`` wraps the four following functions:
+
+- ``model.explicitPred()`` allows to compute the displacement
+ field at :math:`t+1` and a part of the velocity field at :math:`t+1`, denoted by
+ :math:`\vec{\dot{u}^{\st{p}}}_{n+1}`, which will be used later in the method
+ ``model.explicitCorr()``. The equations are:
+
+ .. math::
+ \vec{u}_{n+1} &= \vec{u}_{n} + \Delta t
+ \vec{\dot{u}}_{n} + \frac{\Delta t^2}{2} \vec{\ddot{u}}_{n}\\
+ \vec{\dot{u}^{\st{p}}}_{n+1} &= \vec{\dot{u}}_{n} + \Delta t
+ \vec{\ddot{u}}_{n}
+
+- ``model.updateResidual()`` and ``model.updateAcceleration()`` compute the acceleration increment
+ :math:`\delta \vec{\ddot{u}}`:
+
+ .. math::
+ \left(\mat{M} + \frac{1}{2} \Delta t \mat{C}\right)
+ \delta \vec{\ddot{u}} = \vec{f_{\st{ext}}} - \vec{f}_{\st{int}\, n+1}
+ - \mat{C} \vec{\dot{u}^{\st{p}}}_{n+1} - \mat{M} \vec{\ddot{u}}_{n}
+
+ The internal force :math:`\vec{f}_{\st{int}\, n+1}` is computed from the
+ displacement :math:`\vec{u}_{n+1}` based on the constitutive law.
+
+- ``model.explicitCorr()`` computes the velocity and
+ acceleration fields at :math:`t+1`:
+
+ .. math::
+ \vec{\dot{u}}_{n+1} &= \vec{\dot{u}^{\st{p}}}_{n+1} + \frac{\Delta t}{2}
+ \delta \vec{\ddot{u}} \\ \vec{\ddot{u}}_{n+1} &=
+ \vec{\ddot{u}}_{n} + \delta \vec{\ddot{u}}
+
+The use of an explicit time integration scheme is illustrated by the example:
+``examples/explicit/explicit_dynamic.cc``. This example models the propagation
+of a wave in a steel beam. The beam and the applied displacement in the
+:math:`x` direction are shown in :numref:`fig:smm:explicit`.
+
+
+.. _fig:smm:explicit:
+.. figure:: figures/explicit.svg
+ :align: center
+ :width: 90%
+
+ Numerical setup.
+
+..
+ \begin{figure}[!htb] \centering
+ \begin{tikzpicture}
+ \coordinate (c) at (0,2);
+ \draw[shift={(c)},thick, color=blue] plot [id=x, domain=-5:5, samples=50] ({\x, {(40 * sin(0.1*pi*3*\x) * exp(- (0.1*pi*3*\x)*(0.1*pi*3*\x) / 4))}});
+ \draw[shift={(c)},-latex] (-6,0) -- (6,0) node[right, below] {:math:`x`};
+ \draw[shift={(c)},-latex] (0,-0.7) -- (0,1) node[right] {:math:`u`};
+ \draw[shift={(c)}] (-0.1,0.6) node[left] {:math:`A`}-- (1.5,0.6);
+
+ \coordinate (l) at (0,0.6);
+ \draw[shift={(0,-0.7)}] (-5, 0) -- (5,0) -- (5, 1) -- (-5, 1) -- cycle;
+ \draw[shift={(l)}, latex-latex] (-5,0)-- (5,0) node [midway, above] {:math:`L`};
+ \draw[shift={(l)}] (5,0.2)-- (5,-0.2);
+ \draw[shift={(l)}] (-5,0.2)-- (-5,-0.2);
+
+ \coordinate (h) at (5.3,-0.7);
+ \draw[shift={(h)}, latex-latex] (0,0)-- (0,1) node [midway, right] {:math:`h`};
+ \draw[shift={(h)}] (-0.2,1)-- (0.2,1);
+ \draw[shift={(h)}] (-0.2,0)-- (0.2,0);
+ \end{tikzpicture}
+
+ \caption{Numerical setup \label{fig:smm:explicit}}
+ \end{figure}
+
+The length and height of the beam are :math:`L={10}\textrm{m}` and :math:`h =
+{1}\textrm{m}`, respectively. The material is linear elastic, homogeneous and
+isotropic (density: \SI{7800}{\kilo\gram\per\cubic\metre}, Young's modulus:
+\SI{210}{\giga\pascal} and Poisson's ratio: :math:`0.3`). The imposed
+displacement follow a Gaussian function with a maximum amplitude of :math:`A =
+{0.01}\textrm{m}`. The potential, kinetic and total energies are computed. The
+safety factor is equal to :math:`0.8`.
+
+ \input{manual-constitutive-laws}
+
+Adding a New Constitutive Law
+-----------------------------
+
+There are several constitutive laws in \akantu as described in the
+previous Section~\ref{sect:smm:CL}. It is also possible to use a
+user-defined material for the simulation. These materials are referred
+to as local materials since they are local to the example of the user
+and not part of the \akantu library. To define a new local material,
+two files (``material_XXX.hh`` and ``material_XXX.cc``) have
+to be provided where ``XXX`` is the name of the new material. The
+header file ``material_XXX.hh`` defines the interface of your
+custom material. Its implementation is provided in the
+``material_XXX.cc``. The new law must inherit from the
+``Material`` class or any other existing material class. It is
+therefore necessary to include the interface of the parent material
+in the header file of your local material and indicate the inheritance
+in the declaration of the class::
+
+ /* ---------------------------------------------------------------------- */
+ #include "material.hh"
+ /* ---------------------------------------------------------------------- */
+
+ #ifndef __AKANTU_MATERIAL_XXX_HH__
+ #define __AKANTU_MATERIAL_XXX_HH__
+
+ namespace akantu {
+
+ class MaterialXXX : public Material {
+
+ /// declare here the interface of your material
+
+ };
+
+In the header file the user also needs to declare all the members of the new
+material. These include the parameters that a read from the
+material input file, as well as any other material parameters that will be
+computed during the simulation and internal variables.
+
+
+In the following the example of adding a new damage material will be
+presented. In this case the parameters in the material will consist of the
+Young's modulus, the Poisson coefficient, the resistance to damage and the
+damage threshold. The material will then from these values compute its Lamé
+coefficients and its bulk modulus. Furthermore, the user has to add a new
+internal variable ``damage`` in order to store the amount of damage at each
+quadrature point in each step of the simulation. For this specific material the
+member declaration inside the class will look as follows::
+
+ class LocalMaterialDamage : public Material {
+
+ /// declare constructors/destructors here
+
+ /// declare methods and accessors here
+
+ /* -------------------------------------------------------------------- */
+ /* Class Members */
+ /* -------------------------------------------------------------------- */
+
+ AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real);
+ private:
+
+ /// the young modulus
+ Real E;
+
+ /// Poisson coefficient
+ Real nu;
+
+ /// First Lame coefficient
+ Real lambda;
+
+ /// Second Lame coefficient (shear modulus)
+ Real mu;
+
+ /// resistance to damage
+ Real Yd;
+
+ /// damage threshold
+ Real Sd;
+
+ /// Bulk modulus
+ Real kpa;
+
+ /// damage internal variable
+ InternalField damage;
+
+ };
+
+In order to enable to print the material parameters at any point in
+the user's example file using the standard output stream by typing::
+
+ for (UInt m = 0; m < model.getNbMaterials(); ++m)
+ std::cout << model.getMaterial(m) << std::endl;
+
+the standard output stream operator has to be redefined. This should be done at the end of the header file::
+
+ class LocalMaterialDamage : public Material {
+
+ /// declare here the interace of your material
+
+ }:
+ /* ---------------------------------------------------------------------- */
+ /* inline functions */
+ /* ---------------------------------------------------------------------- */
+ /// standard output stream operator
+ inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this)
+ {
+ _this.printself(stream);
+ return stream;
+ }
+
+However, the user still needs to register the material parameters that
+should be printed out. The registration is done during the call of the
+constructor. Like all definitions the implementation of the
+constructor has to be written in the ``material_XXX.cc``
+file. However, the declaration has to be provided in the
+``material_XXX.hh`` file::
+
+ class LocalMaterialDamage : public Material {
+ /* -------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* -------------------------------------------------------------------- */
+ public:
+
+ LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = "");
+ };
+
+The user can now define the implementation of the constructor in the
+``material_XXX.cc`` file::
+
+ /* ---------------------------------------------------------------------- */
+ #include "local_material_damage.hh"
+ #include "solid_mechanics_model.hh"
+
+ namespace akantu {
+
+ /* ---------------------------------------------------------------------- */
+ LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model,
+ const ID & id) :
+ Material(model, id),
+ damage("damage", *this) {
+ AKANTU_DEBUG_IN();
+
+ this->registerParam("E", E, 0., _pat_parsable, "Young's modulus");
+ this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio");
+ this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient");
+ this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient");
+ this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient");
+ this->registerParam("Yd", Yd, 50., _pat_parsmod);
+ this->registerParam("Sd", Sd, 5000., _pat_parsmod);
+
+ damage.initialize(1);
+
+ AKANTU_DEBUG_OUT();
+ }
+
+During the intializer list the reference to the model and the material id are
+assigned and the constructor of the internal field is called. Inside the scope
+of the constructor the internal values have to be initialized and the
+parameters, that should be printed out, are registered with the function:
+``registerParam``::
+
+ void registerParam(name of the parameter (key in the material file),
+ member variable,
+ default value (optional parameter),
+ access permissions,
+ description);
+
+The available access permissions are as follows:
+- ``_pat_internal``: Parameter can only be output when the material is printed.
+- ``_pat_writable``: User can write into the parameter. The parameter is output when the material is printed.
+- ``_pat_readable``: User can read the parameter. The parameter is output when the material is printed.
+- ``_pat_modifiable``: Parameter is writable and readable.
+- ``_pat_parsable``: Parameter can be parsed, *i.e.* read from the input file.
+- ``_pat_parsmod``: Parameter is modifiable and parsable.
+
+In order to implement the new constitutive law the user needs to
+specify how the additional material parameters, that are not
+defined in the input material file, should be calculated. Furthermore,
+it has to be defined how stresses and the stable time step should be
+computed for the new local material. In the case of implicit
+simulations, in addition, the computation of the tangent stiffness needs
+to be defined. Therefore, the user needs to redefine the following
+functions of the parent material::
+
+ void initMaterial();
+
+ // for explicit and implicit simulations void
+ computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
+
+ // for implicit simulations
+ void computeTangentStiffness(const ElementType & el_type,
+ Array & tangent_matrix,
+ GhostType ghost_type = _not_ghost);
+
+ // for explicit and implicit simulations
+ Real getStableTimeStep(Real h, const Element & element);
+
+In the following a detailed description of these functions is provided:
+- ``initMaterial``: This method is called after the material file is fully read
+ and the elements corresponding to each material are assigned. Some of the
+ frequently used constant parameters are calculated in this method. For
+ example, the Lam\'{e} constants of elastic materials can be considered as
+ such parameters.
+
+- ``computeStress``: In this method, the stresses are computed based on the
+ constitutive law as a function of the strains of the quadrature points. For
+ example, the stresses for the elastic material are calculated based on the
+ following formula:
+ .. math::
+ \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon}
+
+ Therefore, this method contains a loop on all quadrature points assigned to
+ the material using the two macros:
+ ``MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN`` and
+ ``MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END``
+
+ .. code::
+
+ MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type);
+
+ // sigma <- f(grad_u)
+
+ MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
+
+ The strain vector in Akantu contains the values of :math:`\nabla \vec{u}`,
+ i.e. it is really the *displacement gradient*,
+
+- ``computeTangentStiffness``: This method is called when the tangent to the
+ stress-strain curve is desired (see Fig \ref {fig:smm:AL:K}). For example,
+ it is called in the implicit solver when the stiffness matrix for the
+ regular elements is assembled based on the following formula:
+
+ .. math::
+ \label{eqn:smm:constitutive_elasc} \mat{K }
+ =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}}
+
+ Therefore, in this method, the ``tangent`` matrix (\mat{D}) is
+ computed for a given strain.
+
+ The ``tangent`` matrix is a :math:`4^{th}` order tensor which is stored as
+ a matrix in Voigt notation.
+
+ .. _fig:smm:AL:K:
+ .. figure:: figures/tangent.svg
+ :align: center
+ :width: 60%
+
+ Tangent to the stress-strain curve.
+
+..
+ \begin{figure}[!htb]
+ \begin{center}
+ \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf}
+ \caption{Tangent to the stress-strain curve.}
+ \label{fig:smm:AL:K}
+ \end{center}
+ \end{figure}
+
+- ``getCelerity``: The stability criterion of the explicit integration scheme
+ depend on the fastest wave celerity~\eqref{eqn:smm:explicit:stabletime}. This
+ celerity depend on the material, and therefore the value of this velocity
+ should be defined in this method for each new material. By default, the
+ fastest wave speed is the compressive wave whose celerity can be defined in ``getPushWaveSpeed``.
+
+Once the declaration and implementation of the new material has been
+completed, this material can be used in the user's example by including the header file::
+
+ #include "material_XXX.hh"
+
+For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by
+default, the materials are initialized inside the method
+``initFull``. If a local material should be used instead, the
+initialization of the material has to be postponed until the local
+material is registered in the model. Therefore, the model is
+initialized with the boolean for skipping the material initialization
+equal to true::
+
+ /// model initialization
+ model.initFull(_analysis_method = _explicit_lumped_mass);
+
+Once the model has been initialized, the local material needs
+to be registered in the model::
+
+ model.registerNewCustomMaterials("name_of_local_material");
+
+Only at this point the material can be initialized::
+
+ model.initMaterials();
+
+A full example for adding a new damage law can be found in
+``examples/new_material``.
+
+Adding a New Non-Local Constitutive Law
+```````````````````````````````````````
+
+In order to add a new non-local material we first have to add the local
+constitutive law in Akantu (see above). We can then add the non-local version
+of the constitutive law by adding the two files (``material_XXX_non_local.hh``
+and ``material_XXX_non_local.cc``) where ``XXX`` is the name of the
+corresponding local material. The new law must inherit from the two classes,
+non-local parent class, such as the ``MaterialNonLocal`` class, and from the
+local version of the constitutive law, *i.e.* ``MaterialXXX``. It is therefore
+necessary to include the interface of those classes in the header file of your
+custom material and indicate the inheritance in the declaration of the class::
+
+ /* ---------------------------------------------------------------------- */
+ #include "material_non_local.hh" // the non-local parent
+ #include "material_XXX.hh"
+ /* ---------------------------------------------------------------------- */
+
+ #ifndef __AKANTU_MATERIAL_XXX_HH__
+ #define __AKANTU_MATERIAL_XXX_HH__
+
+ namespace akantu {
+
+ class MaterialXXXNonLocal : public MaterialXXX,
+ public MaterialNonLocal {
+
+ /// declare here the interface of your material
+
+ };
+
+As members of the class we only need to add the internal fields to store the
+non-local quantities, which are obtained from the averaging process::
+
+ /* -------------------------------------------------------------------------- */
+ /* Class members */
+ /* -------------------------------------------------------------------------- */
+ protected:
+ InternalField grad_u_nl;
+
+The following four functions need to be implemented in the non-local material::
+
+ /// initialization of the material
+ void initMaterial();
+ /// loop over all element and invoke stress computation
+ virtual void computeNonLocalStresses(GhostType ghost_type);
+ /// compute stresses after local quantities have been averaged
+ virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type)
+ /// compute all local quantities
+ void computeStress(ElementType el_type, GhostType ghost_type);
+
+In the intialization of the non-local material we need to register the local
+quantity for the averaging process. In our example the internal field
+*grad_u_nl* is the non-local counterpart of the gradient of the displacement
+field (*grad_u_nl*)::
+
+ void MaterialXXXNonLocal::initMaterial() {
+ MaterialXXX::initMaterial();
+ MaterialNonLocal::initMaterial();
+ /// register the non-local variable in the manager
+ this->model->getNonLocalManager().registerNonLocalVariable(
+ this->grad_u.getName(),
+ this->grad_u_nl.getName(),
+ spatial_dimension * spatial_dimension);
+ }
+
+The function to register the non-local variable takes as parameters the name of
+the local internal field, the name of the non-local counterpart and the number
+of components of the field we want to average. In the *computeStress* we now
+need to compute all the quantities we want to average. We can then write a loop
+for the stress computation in the function *computeNonLocalStresses* and then
+provide the constitutive law on each integration point in the function
+*computeNonLocalStress*.