diff --git a/.gitignore b/.gitignore
index 050156e04..8074ed792 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,35 +1,36 @@
_skbuild
build*
.dir-locals.el
TAGS
third-party/*/
!third-party/cmake/*
!third-party/akantu_iterators
!third-party/iohelper
+doc/dev-doc/auto_examples/*
*~
release
.*.swp
*.tar.gz
*.tgz
*.tbz
*.tar.bz2
.idea
__pycache__
.mailmap
paraview/*
*.vtu
*.pvd
*.pvtu
*.vtk
compile_commands.json
.clangd
.iwyu.imp
.cache
setup.cfg
.vscode
.auctex*
.clangd
.ccls-cache
.ccls
VERSION
diff --git a/doc/dev-doc/CMakeLists.txt b/doc/dev-doc/CMakeLists.txt
index aaafcce20..c3dfbb244 100644
--- a/doc/dev-doc/CMakeLists.txt
+++ b/doc/dev-doc/CMakeLists.txt
@@ -1,78 +1,78 @@
#===============================================================================
# Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This file is part of Akantu
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with Akantu. If not, see .
#
#===============================================================================
# configured documentation tools and intermediate build results
set(BINARY_BUILD_DIR "${CMAKE_CURRENT_BINARY_DIR}/_build")
# Sphinx cache with pickled ReST documents
set(SPHINX_CACHE_DIR "${CMAKE_CURRENT_BINARY_DIR}/_doctrees")
# HTML output directory
set(SPHINX_HTML_DIR "${CMAKE_CURRENT_BINARY_DIR}/html")
set(SPHINX_OUTPUT "${SPHINX_HTML_DIR}/index.html")
set(SPHINX_INPUT "${CMAKE_CURRENT_BINARY_DIR}/conf.py")
# ---------------------------------------------------------------------------- #
# Sphinx #
# ---------------------------------------------------------------------------- #
find_package(Sphinx REQUIRED)
set(SPHINX_VERBOSE_FLAG "-q")
if(CMAKE_VERBOSE_MAKEFILE)
set(SPHINX_VERBOSE_FLAG)
endif(CMAKE_VERBOSE_MAKEFILE)
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/conf.py"
"${SPHINX_INPUT}"
@ONLY)
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/manual/manual-bibliography.bib"
"${CMAKE_CURRENT_BINARY_DIR}/manual/manual-bibliography.bib"
COPYONLY)
set(SPHINX_PARALLEL_FLAG)
if (SPHINX_VERSION VERSION_GREATER_EQUAL 1.7.0)
set(SPHINX_PARALLEL_FLAG -j auto)
endif()
set(_sphinx_command ${SPHINX_BUILD_EXECUTABLE}
${SPHINX_PARALLEL_FLAG}
- ${SPHINX_VERBOSE_FLAG} -b html
+ ${SPHINX_VERBOSE_FLAG} -D plot_gallery=0 -b html
-c "${CMAKE_CURRENT_BINARY_DIR}"
-d "${SPHINX_CACHE_DIR}"
"${CMAKE_CURRENT_SOURCE_DIR}"
"${SPHINX_HTML_DIR}"
)
file(GLOB_RECURSE _SPHINX_SRCS
"*.rst")
set(ENV(RUNNING_IN_CMAKE) "True")
add_custom_command(
OUTPUT ${SPHINX_OUTPUT}
COMMAND ${_sphinx_command}
DEPENDS ${SPHINX_INPUT} ${_SPHINX_SRCS} "${CMAKE_CURRENT_SOURCE_DIR}/akantu.dox.j2"
COMMENT "Building HTML documentation with Sphinx in ${SPHINX_HTML_DIR}"
)
add_custom_target(sphinx-doc ALL
DEPENDS ${SPHINX_OUTPUT})
diff --git a/doc/dev-doc/akantu.dox.j2 b/doc/dev-doc/akantu.dox.j2
index 321e974ca..5b0b0ddb7 100644
--- a/doc/dev-doc/akantu.dox.j2
+++ b/doc/dev-doc/akantu.dox.j2
@@ -1,54 +1,55 @@
PROJECT_NAME = Akantu
#PROJECT_NUMBER = {{ akantu_version }}
STRIP_FROM_PATH = {{ akantu_source_path }}
STRIP_FROM_INC_PATH = {{ akantu_source_path }}
TAB_SIZE = 4
ALIASES = "rst=\verbatim embed:rst" \
"endrst=\endverbatim"
QUIET = NO
WARN_IF_UNDOCUMENTED = NO
WARN_IF_DOC_ERROR = NO
WARN_AS_ERROR = NO
INPUT = {{ akantu_source_path }}/src \
{{ akantu_source_path }}/test/ci/includes_for_ci
FILE_PATTERNS = *.c *.cc *.hh *.py
EXCLUDE = {{ akantu_source_path }}/src/common/aka_fwd.hh \
- {{ akantu_source_path }}/src/io/dumper/dumpable_dummy.hh
- {{ akantu_source_path }}/src/synchronizer/communicator_dummy.hh
+ {{ akantu_source_path }}/src/io/dumper/dumpable_dummy.hh \
+ {{ akantu_source_path }}/src/synchronizer/communicator_dummy.hh \
{{ akantu_source_path }}/src/common/aka_config.hh.in
RECURSIVE = YES
EXAMPLE_PATH = {{ akantu_source_path }}/examples
EXAMPLE_RECURSIVE = YES
SOURCE_BROWSER = NO
GENERATE_HTML = NO
GENERATE_HTMLHELP = NO
USE_MATHJAX = YES
GENERATE_LATEX = NO
GENERATE_XML = YES
XML_OUTPUT = xml
-DOXYGEN_INCLUDE_PATH = {{ akantu_source_path }}/src/common \
+SEARCH_INCLUDES = YES
+INCLUDE_PATH = {{ akantu_source_path }}/src/common \
{{ akantu_source_path }}/src/fe_engine \
{{ akantu_source_path }}/src/mesh \
{{ akantu_source_path }}/src/model \
{{ akantu_source_path }}/test/ci/includes_for_ci
ENABLE_PREPROCESSING = YES
MACRO_EXPANSION = YES
PREDEFINED = DOXYGEN \
__cplusplus=201402L \
"AKANTU_GET_MACRO(name, value, type) = type get##name() const;" \
"AKANTU_GET_MACRO_NOT_CONST(name, value, type) = type get##name();"
EXPAND_AS_DEFINED = __BEGIN_AKANTU__ \
__END_AKANTU__ \
__BEGIN_AKANTU_DUMPER__ \
__END_AKANTU_DUMPER__ \
AKANTU_SET_MACRO \
AKANTU_GET_MACRO_DEREF_PTR \
AKANTU_GET_MACRO_BY_ELEMENT_TYPE \
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST \
AKANTU_GET_MACRO_AUTO
COLLABORATION_GRAPH = NO
UML_LOOK = YES
TEMPLATE_RELATIONS = YES
CALL_GRAPH = YES
CALLER_GRAPH = YES
LOOKUP_CACHE_SIZE = 0
diff --git a/doc/dev-doc/conf.py b/doc/dev-doc/conf.py
index fb06f77eb..b1ad2b5ce 100644
--- a/doc/dev-doc/conf.py
+++ b/doc/dev-doc/conf.py
@@ -1,331 +1,300 @@
# -*- 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
__copyright__ = (
"Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)"
"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__ = "LGPLv3"
# -- 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 shutil
import jinja2
# import git
# import re
import subprocess
# -- 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",
"sphinxcontrib.bibtex",
"breathe",
"myst_parser",
+ "sphinx_gallery.gen_gallery",
+ "sphinx_copybutton"
]
read_the_docs_build = os.environ.get("READTHEDOCS", None) == "True"
cmake_configure = os.environ.get("RUNNING_IN_CMAKE", None) == "True"
if read_the_docs_build:
akantu_path = "."
akantu_source_path = "../../"
elif cmake_configure:
akantu_path = "@CMAKE_CURRENT_BINARY_DIR@"
akantu_source_path = "@CMAKE_SOURCE_DIR@"
else: # most probably running by hand
akantu_path = "build-doc"
akantu_source_path = "../../"
try:
os.mkdir(akantu_path)
except FileExistsError:
pass
# 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 = 'en'
# 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",
"manual/appendix/elements.rst",
"manual/appendix/material-parameters.rst",
"manual/constitutive-laws.rst",
"manual/new-constitutive-laws.rst",
]
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = "sphinx"
primary_domain = "cpp"
highlight_language = "cpp"
bibtex_bibfiles = ["manual/manual-bibliography.bib"]
# -- Project information -----------------------------------------------------
project = "Akantu"
copyright = (
"2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)"
+ " Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
author = "Nicolas Richart"
-# try:
-# tag_prefix = 'v'
-# git_repo = git.Repo(akantu_source_path)
-
-# git_describe = git_repo.git.describe('--tags', '--dirty', '--always',
-# '--long',
-# '--match', '{}*'.format(tag_prefix))
-
-# print("GIT Describe: {}".format(git_describe))
-
-# # git describe to PEP404 version
-# describe_matches = re.search(
-# (r'^{}(?P.+?)' +
-# r'(?:-(?P\d+)-g(?P[0-9a-f]+)' +
-# r'(?:-(?Pdirty))?)?$').format(tag_prefix),
-# git_describe)
-
-# if describe_matches:
-# describe_matches = describe_matches.groupdict()
-
-# release = describe_matches['version']
-# if describe_matches['distance']:
-# release += '.' if '+' in release else '+'
-# release += '{distance}.{sha}'.format(**describe_matches)
-# if describe_matches['dirty']:
-# release += '.dirty'
-# else:
-# count = git_repo.git.rev_list('HEAD', '--count')
-# describe_matches = re.search(
-# (r'^(?P[0-9a-f]+)' +
-# r'(?:-(?Pdirty))?$').format(tag_prefix),
-# git_describe).groupdict()
-# release = '{}.{}+{}'.format(file_release, count,
-# describe_matches['sha'])
-
-# except git.InvalidGitRepositoryError:
-# with open(os.path.join(akantu_source_path, 'VERSION'), 'r') as fh:
-# version_file = fh.readlines()
-# file_release = version_file[0].strip()
-# release = file_release
-
-
-# print("Release: {} - Version: {}".format(release, version))
-
# -- 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"]
html_logo = "_static/logo_only_akantu.svg"
# 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",
]
}
math_eqref_format = "Eq. {number}"
# MathJax configuration
mathjax3_config = {
"tex": {
"macros": {
"st": ["\\mathrm{#1}", 1],
"mat": ["\\mathbf{#1}", 1],
"half": "\\frac{1}{2}",
},
"packages": {"[+]": ["ams"]},
},
"loader": {"load": ["[tex]/ams"]},
}
# -- 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 -------------------------------------------------
j2_args = {}
if read_the_docs_build or not cmake_configure:
j2_template_path = "."
elif cmake_configure:
j2_template_path = "@CMAKE_CURRENT_SOURCE_DIR@"
os.makedirs(os.path.join(akantu_path, "_static"), exist_ok=True)
shutil.copyfile(
os.path.join("@CMAKE_CURRENT_SOURCE_DIR@", html_logo),
os.path.join(akantu_path, html_logo),
)
j2_args = {
"akantu_source_path": akantu_source_path,
}
print(akantu_path)
j2_env = jinja2.Environment(
loader=jinja2.FileSystemLoader(j2_template_path),
undefined=jinja2.DebugUndefined
)
j2_template = j2_env.get_template("akantu.dox.j2")
with open(os.path.join(akantu_path, "akantu.dox"), "w") as fh:
fh.write(j2_template.render(j2_args))
subprocess.run(["doxygen", "akantu.dox"], cwd=akantu_path)
# print("akantu_path = '{}'".format(akantu_path))
breathe_projects = {"Akantu": os.path.join(akantu_path, "xml")}
breathe_default_project = "Akantu"
breathe_default_members = ("members", "undoc-members")
breathe_implementation_filename_extensions = [".c", ".cc", ".cpp", ".hh"]
breathe_show_enumvalue_initializer = True
-breathe_debug_trace_directives = True
+breathe_debug_trace_directives = False
+breathe_short_warning = True
+
+# -- Gallery ------------------------------------------------------------------
+sphinx_gallery_conf = {
+ 'examples_dirs': os.path.join(akantu_source_path, 'examples'),
+ 'gallery_dirs': os.path.join(akantu_source_path, 'doc', 'dev-doc', 'auto_examples'),
+ 'download_all_examples': False,
+ 'plot_gallery': 'False',
+ 'only_warn_on_example_error': True,
+ 'log_level': {'backreference_missing': 'debug'},
+}
-# -- Options for intersphinx extension ---------------------------------------
+# -- 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/examples.rst b/doc/dev-doc/examples.rst
new file mode 100644
index 000000000..671c7b10c
--- /dev/null
+++ b/doc/dev-doc/examples.rst
@@ -0,0 +1,5 @@
+.. include:: auto_examples/index.rst
+
+
+Tutorials
+=========
diff --git a/doc/dev-doc/index.rst b/doc/dev-doc/index.rst
index 4f7009f4b..653e2c178 100644
--- a/doc/dev-doc/index.rst
+++ b/doc/dev-doc/index.rst
@@ -1,70 +1,71 @@
.. 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
=====================
|license| |joss|
.. |joss| image:: https://joss.theoj.org/papers/3abf3c7945cc9a016a946ce9e02e357f/status.svg
:target: https://joss.theoj.org/papers/3abf3c7945cc9a016a946ce9e02e357f
.. |license| image:: https://img.shields.io/badge/license-LGPLv3-green
:target: https://www.gnu.org/licenses/lgpl-3.0.en.html
**Akantu** means a little element in *Kinyarwanda*, a Bantu language. From now
on it is also an opensource object-oriented *Finite-Element* library which has
the ambition to be generic and efficient. **Akantu** is developed within the
`LSMS `_ (Computational Solid Mechanics Laboratory)*, where
research is conducted at the interface of mechanics, material science, and
scientific computing. The open-source philosophy is important for any scientific
software project evolution. The collaboration permitted by shared codes enforces
sanity when users (and not only developers) can criticize the implementation
details. Akantu was born with the vision to associate genericity, robustness and
efficiency while benefiting the open-source visibility.
.. toctree::
:maxdepth: 2
:caption: User Manual
./manual/getting_started.rst
./manual/basic_types.rst
- ./manual/solidmechanicsmodel.rst
- ./manual/heattransfermodel.rst
- ./manual/contactmechanicsmodel.rst
- ./manual/phasefieldmodel.rst
- ./manual/structuralmechanicsmodel.rst
+ ./manual/models.rst
./manual/io.rst
.. toctree::
- :maxdepth: 2
+ :maxdepth: 1
+ :caption: Examples & Tutorials
+
+ ./auto_examples/index.rst
+
+.. toctree::
+ :maxdepth: 1
:caption: Changelog
./changelog.md
.. toctree::
:maxdepth: 2
:caption: API Reference
./reference.rst
.. toctree::
- :maxdepth: 2
+ :maxdepth: 1
:caption: Appendix
./manual/appendix.rst
.. toctree::
:maxdepth: 2
:caption: Bibliography
./manual/bibliography.rst
-Indices and tables
-==================
+.. toctree::
+ :maxdepth: 1
+ :caption: Indices and tables
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
+ ./manual/extra.rst
diff --git a/doc/dev-doc/manual/constitutive-laws.rst b/doc/dev-doc/manual/constitutive-laws.rst
index 6b76d5957..1a28ed49d 100644
--- a/doc/dev-doc/manual/constitutive-laws.rst
+++ b/doc/dev-doc/manual/constitutive-laws.rst
@@ -1,1407 +1,1407 @@
.. _sect-smm-cl:
Constitutive Laws
------------------
+`````````````````
In order to compute an element’s response to deformation, one needs to
use an appropriate constitutive relationship. The constitutive law is
used to compute the element’s stresses from the element’s strains.
In the finite-element discretization, the constitutive formulation is
applied to every quadrature point of each element. When the implicit
formulation is used, the tangent matrix has to be computed.
| The chosen materials for the simulation have to be specified in the
mesh file or, as an alternative, they can be assigned using the at
``element_material`` vector. For
every material assigned to the problem one has to specify the material
characteristics (constitutive behavior and material properties) using
the text input file (see :ref:`sect-io-material`).
| In order to conveniently store values at each quadrature in a material point
``Akantu`` provides a special data structure, the at :cpp:class:`InternalField
`. The internal fields are inheriting from the at
:cpp:class:`ElementTypeMapArray `. Furthermore,
it provides several functions for initialization, auto-resizing and auto
removal of quadrature points.
The constitutive law is precised within an input file. The dedicated material section
is then read by :cpp:func:`initFull `
method of :cpp:class:`SolidMechanicsModel ` which
initializes the different materials specified with the following convention
.. code-block:: python
material constitutive_law [
name = value
rho = value
...
]
where *constitutive_law* is the adopted constitutive law, followed by
the material properties listed one by line in the bracket (*e.g.*,
``name`` and density :math:``rho``. Some constitutive laws can also
have an *optional flavor*.
For certain materials, it is possible to activate the large deformation
strain and stress evaluations. Internally the strain measure becomes
the right Cauchy–Green deformation tensor and the evaluated stress
measure becomes the Piola-Kirchhoff stress tensor. This is activated by using:
.. code-block:: python
material constitutive_law [
finite_deformation = true # Activates the large deformation routines (bool)
...
]
Sometimes it is also desired to generate random distributions of
internal parameters. An example might be the critical stress at which
the material fails. To generate such a field, in the text input file, a
random quantity needs be added to the base value:
All parameters are real numbers. For the uniform distribution, minimum
and maximum values have to be specified. Random parameters are defined
as a :math:`base` value to which we add a random number that follows the
chosen distribution.
The
`Uniform `__
distribution is gives a random values between in :math:`[min, max)`. The
`Weibull `__
distribution is characterized by the following cumulative distribution
function:
.. math:: F(x) = 1- e^{-\left({x/\lambda}\right)^m}
which depends on :math:`m` and :math:`\lambda`, which are the shape
parameter and the scale parameter. These random distributions are
different each time the code is executed. In order to obtain always the
same one, it possible to manually set the *seed* that is the number from
which these pseudo-random distributions are created. This can be done by
adding the following line to the input file *outside* the material
parameters environments:
.. code-block::
seed = 1.0
where the value 1.0 can be substituted with any number. Currently
``Akantu`` can reproduce always the same distribution when the seed is
specified *only* in serial. The value of the *seed* can be also
specified directly in the code (for instance in the main file) with the
command:
.. code-block::
RandGenerator::seed(1.0)
The same command, with empty brackets, can be used to check the value of
the *seed* used in the simulation.
The following sections describe the constitutive models implemented in
``Akantu``.
-----
Elastic
-```````
+'''''''
The elastic law is a commonly used constitutive relationship that can be
used for a wide range of engineering materials (*e.g.*, metals,
concrete, rock, wood, glass, rubber, etc.) provided that the strains
remain small (*i.e.*, small deformation and stress lower than yield
strength).
The elastic laws are often expressed as
:math:`\boldsymbol{\sigma} =
\boldsymbol{C}:\boldsymbol{\varepsilon}` with
where :math:`\boldsymbol{\sigma}` is the Cauchy stress
tensor, :math:`\boldsymbol{\varepsilon}` represents the
infinitesimal strain tensor and :math:`\boldsymbol{C}` is
the elastic modulus tensor.
.. _sect-smm-linear-elastic-isotropic:
Linear isotropic
-''''''''''''''''
+""""""""""""""""
Keyword: **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material elastic [
name = steel
rho = 7800 # density (Real)
E = 2.1e11 # young's modulus (Real)
nu = 0.3 # poisson's ratio (Real)
Plane_stress = false # Plane stress simplification (only 2D problems) (bool)
finite_deformation = false # activates the evaluation of strains with green's tensor (bool)
]
Available energies Energies:
- ``potential``: elastic potential energy
-----
The linear isotropic elastic behavior is described by Hooke’s law, which
states that the stress is linearly proportional to the applied strain
(material behaves like an ideal spring), as illustrated in
:numref:`fig:smm:cl:el`.
.. figure:: figures/cl/stress_strain_el.svg
:alt: Elastic strain-stress curve
:name: fig:smm:cl:el
:width: 60.0%
Stress-strain curve of elastic material and schematic representation of
Hooke's law, denoted as a spring.
The equation that relates the strains to the displacements is: point)
from the displacements as follows:
.. math::
\label{eqn:smm:strain_inf}
\boldsymbol{\varepsilon} =
\frac{1}{2} \left[ \nabla_0 \boldsymbol{u}+\nabla_0 \boldsymbol{u}^T \right]
where :math:`\boldsymbol{\varepsilon}` represents the
infinitesimal strain tensor,
:math:`\nabla_{0}\boldsymbol{u}` the displacement gradient
tensor according to the initial configuration. The constitutive equation
for isotropic homogeneous media can be expressed as:
.. math::
\label{eqn:smm:material:constitutive_elastic}
\boldsymbol{\sigma } =\lambda\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{I}+2 \mu\boldsymbol{\varepsilon}
where :math:`\boldsymbol{\sigma}` is the Cauchy stress
tensor (:math:`\lambda` and :math:`\mu` are the the first and second
Lame’s coefficients).
In Voigt notation this correspond to
.. math::
\begin{aligned}
\left[\begin{array}{c}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{33}\\
\sigma_{23}\\
\sigma_{13}\\
\sigma_{12}\\
\end{array}\right]
&= \frac{E}{(1+\nu)(1-2\nu)}\left[
\begin{array}{cccccc}
1-\nu & \nu & \nu & 0 & 0 & 0\\
\nu & 1-\nu & \nu & 0 & 0 & 0\\
\nu & \nu & 1-\nu & 0 & 0 & 0\\
0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\
0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \\
\end{array}\right]
\left[\begin{array}{c}
\varepsilon_{11}\\
\varepsilon_{22}\\
\varepsilon_{33}\\
2\varepsilon_{23}\\
2\varepsilon_{13}\\
2\varepsilon_{12}\\
\end{array}\right]\end{aligned}
This formulation is not sufficient to represent all elastic material
behavior. Some materials have characteristic orientation that have to be
taken into account. To represent this anisotropy a more general
stress-strain law has to be used, as shown below.
-----
.. _sect-smm-linear-elastic-anisotropic:
Linear anisotropic
-''''''''''''''''''
+""""""""""""""""""
Keyword: **elastic_anisotropic**
Material description with input file:
.. code-block:: python
#input.dat
material elastic_anisotropic [
name = aluminum
rho = 1.6465129043044597 # density (Real)
C11 = 105.092023 # Coefficient ij of material tensor C (Real)
C12 = 59.4637759 # all the 36 values
C13 = 59.4637759 # in Voigt notation can be entered
C14 = 0 # zero coefficients can be omited
C15 = 0
C16 = 0
C22 = 105.092023
C23 = 59.4637759
C24 = 0
C25 = 0
C26 = 0
C33 = 105.092023
C34 = 0
C35 = 0
C36 = 0
C44 = 30.6596356
C45 = 0
C46 = 0
C55 = 30.6596356
C56 = 0
C66 = 30.6596356
n1 = [-1, 1, 0] # Direction of first material axis (Vector)
n2 = [ 1, 1, 1] # Direction of second material axis (Vector)
n3 = [ 1, 1, -2] # Direction of thrid material axis (Vector)
]
----
We define the elastic modulus tensor as follows:
.. math::
\begin{aligned}
\left[\begin{array}{c}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{33}\\
\sigma_{23}\\
\sigma_{13}\\
\sigma_{12}\\
\end{array}\right]
&= \left[
\begin{array}{cccccc}
c_{11} & c_{12} & c_{13} & c_{14} & c_{15} & c_{16}\\
c_{21} & c_{22} & c_{23} & c_{24} & c_{25} & c_{26}\\
c_{31} & c_{32} & c_{33} & c_{34} & c_{35} & c_{36}\\
c_{41} & c_{42} & c_{43} & c_{44} & c_{45} & c_{46}\\
c_{51} & c_{52} & c_{53} & c_{54} & c_{55} & c_{56}\\
c_{61} & c_{62} & c_{63} & c_{64} & c_{65} & c_{66}\\
\end{array}\right]
\left[\begin{array}{c}
\varepsilon_{11}\\
\varepsilon_{22}\\
\varepsilon_{33}\\
2\varepsilon_{23}\\
2\varepsilon_{13}\\
2\varepsilon_{12}\\
\end{array}\right]\end{aligned}
To simplify the writing of input files the :math:`\boldsymbol{C}` tensor
is expressed in the material basis. And this basis as to be given too.
This basis :math:`\Omega_{{\mathrm{mat}}}
= \{\boldsymbol{n_1}, \boldsymbol{n_2}, \boldsymbol{n_3}\}`
is used to define the rotation :math:`R_{ij} =
\boldsymbol{n_j} . \boldsymbol{e_i}`. And
:math:`\boldsymbol{C}` can be rotated in the global basis
:math:`\Omega
= \{\boldsymbol{e_1}, \boldsymbol{e_2}, \boldsymbol{e_3}\}`
as follow:
.. math::
\begin{aligned}
\boldsymbol{C}_{\Omega} &= \boldsymbol{R}_1 \boldsymbol{C}_{\Omega_{{\mathrm{mat}}}} \boldsymbol{R}_2\\
\boldsymbol{R}_1 &= \left[
\begin{array}{cccccc}
R_{11} R_{11} & R_{12} R_{12} & R_{13} R_{13} & R_{12} R_{13} & R_{11} R_{13} & R_{11} R_{12}\\
R_{21} R_{21} & R_{22} R_{22} & R_{23} R_{23} & R_{22} R_{23} & R_{21} R_{23} & R_{21} R_{22}\\
R_{31} R_{31} & R_{32} R_{32} & R_{33} R_{33} & R_{32} R_{33} & R_{31} R_{33} & R_{31} R_{32}\\
R_{21} R_{31} & R_{22} R_{32} & R_{23} R_{33} & R_{22} R_{33} & R_{21} R_{33} & R_{21} R_{32}\\
R_{11} R_{31} & R_{12} R_{32} & R_{13} R_{33} & R_{12} R_{33} & R_{11} R_{33} & R_{11} R_{32}\\
R_{11} R_{21} & R_{12} R_{22} & R_{13} R_{23} & R_{12} R_{23} & R_{11} R_{23} & R_{11} R_{22}\\
\end{array}\right]\\
\boldsymbol{R}_2 &= \left[
\begin{array}{cccccc}
R_{11} R_{11} & R_{21} R_{21} & R_{31} R_{31} & R_{21} R_{31} & R_{11} R_{31} & R_{11} R_{21}\\
R_{12} R_{12} & R_{22} R_{22} & R_{32} R_{32} & R_{22} R_{32} & R_{12} R_{32} & R_{12} R_{22}\\
R_{13} R_{13} & R_{23} R_{23} & R_{33} R_{33} & R_{23} R_{33} & R_{13} R_{33} & R_{13} R_{23}\\
R_{12} R_{13} & R_{22} R_{23} & R_{32} R_{33} & R_{22} R_{33} & R_{12} R_{33} & R_{12} R_{23}\\
R_{11} R_{13} & R_{21} R_{23} & R_{31} R_{33} & R_{21} R_{33} & R_{11} R_{33} & R_{11} R_{23}\\
R_{11} R_{12} & R_{21} R_{22} & R_{31} R_{32} & R_{21} R_{32} & R_{11} R_{32} & R_{11} R_{22}\\
\end{array}\right]\\\end{aligned}
-----
.. _sect-smm-linear-elastic-orthotropic:
Linear orthotropic
-''''''''''''''''''
+""""""""""""""""""
Keyword: **elastic_orthotropic**
Inherits from **elastic_anisotropic**
Material description with input file:
.. code-block:: python
#input.dat
material elastic_orthotropic [
name = test_mat_1
rho = 1 # density
n1 = [-1, 1, 0] # Direction of first material axis (Vector)
n2 = [ 1, 1, 1] # Direction of second material axis (Vector)
n3 = [ 1, 1, -2] # Direction of thrid material axis (Vector)
E1 = 1 # Young's modulus in direction n1 (Real)
E2 = 2 # Young's modulus in direction n2 (Real)
E3 = 3 # Young's modulus in direction n3 (Real)
nu12 = 0.1 # Poisson's ratio 12 (Real)
nu13 = 0.2 # Poisson's ratio 13 (Real)
nu23 = 0.3 # Poisson's ratio 23 (Real)
G12 = 0.5 # Shear modulus 12 (Real)
G13 = 1 # Shear modulus 13 (Real)
G23 = 2 # Shear modulus 23 (Real)
]
-----
A particular case of anisotropy is when the material basis is orthogonal
in which case the elastic modulus tensor can be simplified and rewritten
in terms of 9 independents material parameters.
.. math::
\begin{aligned}
\left[\begin{array}{c}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{33}\\
\sigma_{23}\\
\sigma_{13}\\
\sigma_{12}\\
\end{array}\right]
&= \left[
\begin{array}{cccccc}
c_{11} & c_{12} & c_{13} & 0 & 0 & 0 \\
& c_{22} & c_{23} & 0 & 0 & 0 \\
& & c_{33} & 0 & 0 & 0 \\
& & & c_{44} & 0 & 0 \\
& \text{sym.} & & & c_{55} & 0 \\
& & & & & c_{66}\\
\end{array}\right]
\left[\begin{array}{c}
\varepsilon_{11}\\
\varepsilon_{22}\\
\varepsilon_{33}\\
2\varepsilon_{23}\\
2\varepsilon_{13}\\
2\varepsilon_{12}\\
\end{array}\right]\end{aligned}
.. math::
\begin{aligned}
c_{11} &= E_1 (1 - \nu_{23}\nu_{32})\Gamma \qquad c_{22} = E_2 (1 - \nu_{13}\nu_{31})\Gamma \qquad c_{33} = E_3 (1 - \nu_{12}\nu_{21})\Gamma\\
c_{12} &= E_1 (\nu_{21} - \nu_{31}\nu_{23})\Gamma = E_2 (\nu_{12} - \nu_{32}\nu_{13})\Gamma\\
c_{13} &= E_1 (\nu_{31} - \nu_{21}\nu_{32})\Gamma = E_2 (\nu_{13} - \nu_{21}\nu_{23})\Gamma\\
c_{23} &= E_2 (\nu_{32} - \nu_{12}\nu_{31})\Gamma = E_3 (\nu_{23} - \nu_{21}\nu_{13})\Gamma\\
c_{44} &= \mu_{23} \qquad c_{55} = \mu_{13} \qquad c_{66} = \mu_{12} \\
\Gamma &= \frac{1}{1 - \nu_{12} \nu_{21} - \nu_{13} \nu_{31} - \nu_{32} \nu_{23} - 2 \nu_{21} \nu_{32} \nu_{13}}\end{aligned}
The Poisson ratios follow the rule
:math:`\nu_{ij} = \nu_{ji} E_i / E_j`.
-----
.. _sect-smm-cl-neohookean:
Neo-Hookean
-'''''''''''
+"""""""""""
Keyword: **neohookean**
Inherits from **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material neohookean [
name = material_name
rho = 7800 # density (Real)
E = 2.1e11 # young's modulus (Real)
nu = 0.3 # poisson's ratio (Real)
]
-----
The hyperelastic Neo-Hookean constitutive law results from an extension
of the linear elastic relationship (Hooke’s Law) for large deformation.
Thus, the model predicts nonlinear stress-strain behavior for bodies
undergoing large deformations.
.. figure:: figures/cl/stress_strain_neo.svg
:alt: Neo-hookean Stress-strain curve.
:name: fig:smm:cl:neo_hookean
:width: 40.0%
Neo-hookean Stress-strain curve.
As illustrated in :numref:`fig:smm:cl:neo_hookean`, the behavior
is initially linear and the mechanical behavior is very close to the
corresponding linear elastic material. This constitutive relationship,
which accounts for compressibility, is a modified version of the one
proposed by Ronald Rivlin :cite:`Belytschko:2000`.
The strain energy stored in the material is given by:
.. math::
\label{eqn:smm:constitutive:neohookean_potential}
\Psi(\boldsymbol{C}) = \frac{1}{2}\lambda_0\left(\ln J\right)^2-\mu_0\ln J+\frac{1}{2}
\mu_0\left(\mathrm{tr}(\boldsymbol{C})-3\right)
where :math:`\lambda_0` and :math:`\mu_0` are, respectively, Lamé’s
first parameter and the shear modulus at the initial configuration.
:math:`J` is the jacobian of the deformation gradient
(:math:`\boldsymbol{F}=\nabla_{\!\!\boldsymbol{X}}\boldsymbol{x}`):
:math:`J=\text{det}(\boldsymbol{F})`. Finally
:math:`\boldsymbol{C}` is the right Cauchy-Green
deformation tensor.
Since this kind of material is used for large deformation problems, a
finite deformation framework should be used. Therefore, the Cauchy
stress (:math:`\boldsymbol{\sigma}`) should be computed
through the second Piola-Kirchhoff stress tensor
:math:`\boldsymbol{S}`:
.. math:: \boldsymbol{\sigma } = \frac{1}{J}\boldsymbol{F}\boldsymbol{S}\boldsymbol{F}^T
Finally the second Piola-Kirchhoff stress tensor is given by:
.. math::
\boldsymbol{S} = 2\frac{\partial\Psi}{\partial\boldsymbol{C}} = \lambda_0\ln J
\boldsymbol{C}^{-1}+\mu_0\left(\boldsymbol{I}-\boldsymbol{C}^{-1}\right)
The parameters to indicate in the material file are the same as those
for the elastic case: ``E`` (Young’s modulus), ``nu`` (Poisson’s ratio).
-----
.. _sect-smm-cl-ve:
Visco-Elastic
-`````````````
+'''''''''''''
.. _sect-smm-cl-sls:
Standard-Linear Solid
-'''''''''''''''''''''
+"""""""""""""""""""""
Keyword: **sls_deviatoric**
Inherits from **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material sls_deviatoric [
name = material_name
rho = 1000 # density (Real)
E = 2.1e9 # young's modulus (Real)
nu = 0.4 # poisson's ratio (Real)
Eta = 1. # Viscosity (Real)
Ev = 0.5 # Stiffness of viscous element (Real)
Plane_stress = false # Plane stress simplification (bool, only 2D problems)
]
Since this material inherits from :cpp:class:`MaterialElastic
` the parameter :math:`E` preceeds :math:`E_{\mathrm{inf}}`. So only :math:`E`
and :math:`E_v` can be set, :math:`E_{\mathrm{inf}}` is deduced by :math:`E_{\mathrm{inf}} = E - E_{v}`.
-----
Visco-elasticity is characterized by strain rate dependent behavior.
Moreover, when such a material undergoes a deformation it dissipates
energy. This dissipation results in a hysteresis loop in the
stress-strain curve at every loading cycle (see
:numref:`fig:smm:cl:visco-elastic:hyst`).
In principle, it can be applied to many materials, since all materials
exhibit a visco-elastic behavior if subjected to particular conditions
(such as high temperatures).
.. figure:: figures/cl/stress_strain_visco.svg
:name: fig:smm:cl:visco-elastic:hyst
:align: center
:width: 40.0%
Characteristic stress-strain behavior of a visco-elastic material with hysteresis loop
.. figure:: figures/cl/visco_elastic_law.svg
:name: fig:smm:cl:visco-elastic:model
:align: center
:width: 40.0%
Schematic representation of the standard rheological linear solid visco-elastic model
The standard rheological linear solid model (see Sections 10.2 and 10.3
of :cite:`simo92`) has been implemented in ``Akantu``. This
model results from the combination of a spring mounted in parallel with
a spring and a dashpot connected in series, as illustrated in
:numref:`fig:smm:cl:visco-elastic:model`.
The advantage of this model is that it allows to account for creep or
stress relaxation. The equation that relates the stress to the strain is
(in 1D):
.. math::
\frac{d\varepsilon(t)}{dt} = \left ( E_{\mathrm{inf}} + E_{v} \right ) ^ {-1} \cdot \left [ \frac{d\sigma(t)}{dt} + \frac{E_{v}}{\eta}\sigma(t) - \frac{E_{\mathrm{inf}}E_V}{\eta}\varepsilon(t) \right ]
where :math:`\eta` is the viscosity. The equilibrium condition is unique and is
attained in the limit, as :math:`t \to \infty`. At this stage, the response is
elastic and depends on the Young’s modulus :math:`E`. The mandatory parameters
for the material file are the following: ``rho`` (density), ``E`` (Young’s
modulus), ``nu`` (Poisson’s ratio), ``Plane_Stress`` (if set to zero plane
strain, otherwise plane stress), ``eta`` (dashpot viscosity) and ``Ev``
(stiffness of the viscous element).
Note that the current standard linear solid model is applied only on the
deviatoric part of the strain tensor. The spheric part of the strain
tensor affects the stress tensor like an linear elastic material.
-----
.. _sect-smm-cl-maxwell:
Maxwell Chain Visco-Elasticity
-''''''''''''''''''''''''''''''
+""""""""""""""""""""""""""""""
Keyword: **viscoelastic_maxwell**
Inherits from **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material viscoelastic_maxwell [
name = material_name
rho = 1000 # density (Real)
Einf = 5.e9 # Infinite time Young's modulus (Real)
nu = 0.4 # poisson's ratio (Real)
Ev = [1.e9, 2.e9, 3.e9] # Maxwell elements' stiffness values (Vector)
Eta = [1.e14, 2.e16, 3.e16] # Dashpot elements' viscosity values (Vector)
Plane_stress = false # Plane stress simplification (bool, only 2D problems)
]
----
.. figure:: figures/cl/maxwell_chain.png
:name: fig-smm-cl-visco-elastic-maxwell
:align: center
:width: 40.0%
Schematic representation of the Maxwell chain
A different visco-elastic rheological model available to users is the
generalized Maxwell chain (see :cite:`de_borst_finiteelement_1994` and Section 46.7.4 of :cite:`diana_manual`).
It consists of a series of sequential spring-dashpots
(Maxwell elements) placed in parallel with one single spring (see
:numref:`fig-smm-cl-visco-elastic-maxwell`). The relation between stresses and
strain comes from
.. math:: \sigma \left ( t \right ) = \int_{-\infty}^{t} E \left ( t, \tau \right ) \mathbf{D} \dot{\varepsilon} d\tau
where :math:`E(t,\tau)` is the time-dependent relaxation function, :math:`\tau`
is the loading age, and :math:`\mathbf{D}` is the dimensionless matrix relating
a 3D deformation state to a 1D relaxation function. The relaxation function is
expanded in the exponential series
.. math::
:label: eqn-relaxation-function
E \left ( t, \tau \right ) = E_{0} + \sum_{\alpha=1}^{n} E_{\alpha} e^{- \frac{t- \tau}{\lambda_{\alpha}}}
where the relaxation time of each Maxwell element is defined as
:math:`\lambda_{\alpha}=\eta_{\alpha} / E_{\alpha}` with :math:`\eta_{\alpha}`
being the viscosity of a dash-pot. Assuming a constant strain rate within each
time step, the analytical integration of the right-hand side of
:eq:`eqn-relaxation-function` leads to the following form
.. math::
\sigma \left ( t + \Delta t \right ) = E_{0} \mathbf{D} \varepsilon +
\sum_{\alpha=1}^n \left ( \left ( 1 - e^{\frac{- \Delta t}{\lambda_{\alpha}}} \right )
\frac{E_{\alpha} \lambda_{\alpha}}{\Delta t} \mathbf{D} \delta \varepsilon +
e^{\frac{-\Delta t}{\lambda_{\alpha}}} \sigma_{\alpha} \left ( t \right ) \right )
with :math:`\sigma_{\alpha}(t)` being the internal stress within each Maxwell
element, defined as
.. math:: \sigma_{\alpha} \left ( t \right ) = \mathbf{D} \int_0^t E_{\alpha} e^{\frac{-t- \tau}{\lambda_{\alpha}}} \dot{\varepsilon} d \tau
The first term under the sum sign in above equation could be seen as the
effective stiffness of a single Maxwell element multiplied by the matrix
:math:`\mathbf{D}` and the strain increment :math:`\Delta \varepsilon`:
.. math::
E_{\alpha}^{ef} = \left ( 1- e^{\frac{-Δt}{λ_α}} \right ) \frac{E_α λ_α}{Δt}
Time increment :math:`Δt` controls the rate dependency of the effective
stiffness. By limit analysis, we find the limiting values of the effective
stiffness which are equal to :math:`E_0` for infinitely slow loading (:math:`Δt`
tending to infinity) and :math:`E_0+ΣE_α` for infinitely fast (:math:`Δt` tending to
0). At the end of each converged time step, the internal stress
:math:`σ_α(t)` is updated according to
.. math::
σ_α \left ( t \right ) = σ_α \left ( t - Δt \right ) e^{\frac{-Δt}{λ_α}} + E_α^{ef} \mathbf{D} Δε
The mandatory parameters for the material file are the following: ``rho``
(density), ``nu`` (Poisson’s ratio), ``Plane_Stress`` (if set to zero plane
strain, otherwise plane stress), ``Einf`` (infinite time Young’s modulus),
``Ev`` (Maxwell elements' stiffness values stored in a vector), ``Eta``
(dashpots' viscosity values stored in a vector).
The Maxwell model is applied on the entire strain tensor and does not
distinguish between its deviatoric and hydrostatic components. Note that the
time step has to be specified for the model using current material both for
static and dynamic simulations:
.. code-block:: c++
model.setTimeStep(time_step_value);
.. _sect-smm-cl-plastic:
Plastic
-```````
+'''''''
Small-Deformation Plasticity
-''''''''''''''''''''''''''''
+""""""""""""""""""""""""""""
Keyword: **plastic_linear_isotropic_hardening**
Inherits from **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material plastic_linear_isotropic_hardening [
name = material_name
rho = 1000 # density (Real)
E = 2.1e9 # young's modulus (Real)
nu = 0.4 # poisson's ratio (Real)
h = 0.1 # Hardening modulus (Real)
sigma_y = 1e6 # Yield stress (Real)
]
Energies:
- ``potential``: elastic part of the potential energy
- ``plastic``: dissipated plastic energy (integrated over time)
-----
The small-deformation plasticity is a simple plasticity material
formulation which accounts for the additive decomposition of strain into
elastic and plastic strain components. This formulation is applicable to
infinitesimal deformation where the additive decomposition of the strain
is a valid approximation. In this formulation, plastic strain is a
shearing process where hydrostatic stress has no contribution to
plasticity and consequently plasticity does not lead to volume change.
:numref:`fig:smm:cl:Lin-strain-hard` shows the linear strain
hardening elasto-plastic behavior according to the additive
decomposition of strain into the elastic and plastic parts in
infinitesimal deformation as
.. math::
\boldsymbol{\varepsilon} &= \boldsymbol{\varepsilon}^e +\boldsymbol{\varepsilon}^p\\
\boldsymbol{\sigma} &= 2G(\boldsymbol{\varepsilon}^e) + \lambda \mathrm{tr}(\boldsymbol{\varepsilon}^e)\boldsymbol{I}
.. figure:: figures/cl/isotropic_hardening_plasticity.svg
:name: fig:smm:cl:Lin-strain-hard
:align: center
Stress-strain curve for the small-deformation plasticity with linear isotropic hardening.
In this class, the von Mises yield criterion is used. In the von Mises
yield criterion, the yield is independent of the hydrostatic stress.
Other yielding criteria such as Tresca and Gurson can be easily
implemented in this class as well.
In the von Mises yield criterion, the hydrostatic stresses have no
effect on the plasticity and consequently the yielding occurs when a
critical elastic shear energy is achieved.
.. math::
\label{eqn:smm:constitutive:von Mises}
f = \sigma_{{\mathrm{eff}}} - \sigma_y = \left(\frac{3}{2} {\boldsymbol{\sigma}}^{{\mathrm{tr}}} : {\boldsymbol{\sigma}}^{{\mathrm{tr}}}\right)^\frac{1}{2}-\sigma_y (\boldsymbol{\varepsilon}^p)
.. math::
\label{eqn:smm:constitutive:yielding}
f < 0 \quad \textrm{Elastic deformation,} \qquad f = 0 \quad \textrm{Plastic deformation}
where :math:`\sigma_y` is the yield strength of the material which can
be function of plastic strain in case of hardening type of materials and
:math:`{\boldsymbol{\sigma}}^{{\mathrm{tr}}}` is the
deviatoric part of stress given by
.. math::
\label{eqn:smm:constitutive:deviatoric stress}
{\boldsymbol{\sigma}}^{{\mathrm{tr}}}=\boldsymbol{\sigma} - \frac{1}{3} \mathrm{tr}(\boldsymbol{\sigma}) \boldsymbol{I}
After yielding :math:`(f = 0)`, the normality hypothesis of plasticity
determines the direction of plastic flow which is normal to the tangent
to the yielding surface at the load point. Then, the tensorial form of
the plastic constitutive equation using the von Mises yielding criterion
(see equation 4.34) may be written as
.. math::
\label{eqn:smm:constitutive:plastic contitutive equation}
\Delta {\boldsymbol{\varepsilon}}^p = \Delta p \frac {\partial{f}}{\partial{\boldsymbol{\sigma}}}=\frac{3}{2} \Delta p \frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}}
In these expressions, the direction of the plastic strain increment (or
equivalently, plastic strain rate) is given by
:math:`\frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}}`
while the magnitude is defined by the plastic multiplier
:math:`\Delta p`. This can be obtained using the *consistency condition*
which impose the requirement for the load point to remain on the
yielding surface in the plastic regime.
Here, we summarize the implementation procedures for the
small-deformation plasticity with linear isotropic hardening:
#. Compute the trial stress:
.. math:: {\boldsymbol{\sigma}}^{{\mathrm{tr}}} = {\boldsymbol{\sigma}}_t + 2G\Delta \boldsymbol{\varepsilon} + \lambda \mathrm{tr}(\Delta \boldsymbol{\varepsilon})\boldsymbol{I}
#. Check the Yielding criteria:
.. math:: f = (\frac{3}{2} {\boldsymbol{\sigma}}^{{\mathrm{tr}}} : {\boldsymbol{\sigma}}^{{\mathrm{tr}}})^{1/2}-\sigma_y (\boldsymbol{\varepsilon}^p)
#. Compute the Plastic multiplier:
.. math::
\begin{aligned}
d \Delta p &= \frac{\sigma^{tr}_{eff} - 3G \Delta P^{(k)}- \sigma_y^{(k)}}{3G + h}\\
\Delta p^{(k+1)} &= \Delta p^{(k)}+ d\Delta p\\
\sigma_y^{(k+1)} &= (\sigma_y)_t+ h\Delta p
\end{aligned}
#. Compute the plastic strain increment:
.. math:: \Delta {\boldsymbol{\varepsilon}}^p = \frac{3}{2} \Delta p \frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}}
#. Compute the stress increment:
.. math:: {\Delta \boldsymbol{\sigma}} = 2G(\Delta \boldsymbol{\varepsilon}-\Delta \boldsymbol{\varepsilon}^p) + \lambda \mathrm{tr}(\Delta \boldsymbol{\varepsilon}-\Delta \boldsymbol{\varepsilon}^p)\boldsymbol{I}
#. Update the variables:
.. math::
\begin{aligned}
{\boldsymbol{\varepsilon^p}} &= {\boldsymbol{\varepsilon}}^p_t+{\Delta {\boldsymbol{\varepsilon}}^p}\\
{\boldsymbol{\sigma}} &= {\boldsymbol{\sigma}}_t+{\Delta \boldsymbol{\sigma}}
\end{aligned}
We use an implicit integration technique called *the radial return method* to
obtain the plastic multiplier. This method has the advantage of being
unconditionally stable, however, the accuracy remains dependent on the step
size. The plastic parameters to indicate in the material file are:
:math:`\sigma_y` (Yield stress) and ``h`` (Hardening modulus). In addition, the
elastic parameters need to be defined as previously mentioned: ``E`` (Young’s
modulus), ``nu`` (Poisson’s ratio).
-----
Damage
-``````
+""""""
In the simplified case of a linear elastic and brittle material,
isotropic damage can be represented by a scalar variable :math:`d`,
which varies from :math:`0` to :math:`1` for no damage to fully broken
material respectively. The stress-strain relationship then becomes:
.. math:: \boldsymbol{\sigma} = (1-d)\, \boldsymbol{C}:\boldsymbol{\varepsilon}
where :math:`\boldsymbol{\sigma}`,
:math:`\boldsymbol{\varepsilon}` are the Cauchy stress and
strain tensors, and :math:`\boldsymbol{C}` is the elastic
stiffness tensor. This formulation relies on the definition of an
evolution law for the damage variable. In ``Akantu``, many possibilities
exist and they are listed below.
----
.. _sect-smm-cl-damage-marigo:
Marigo
-''''''
+""""""
Keyword: **marigo**
Inherits from **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material marigo [
name = material_name
rho = 1000 # density (Real)
E = 2.1e9 # young's modulus (Real)
nu = 0.4 # poisson's ratio (Real)
Plane_stress = false # Plane stress simplification (bool, only 2D problems)
Yd = 0.1 # Hardening modulus (Random)
Sd = 1. # Damage energy (Real)
]
Energies:
- ``dissipated``: energy dissipated in damage
-----
This damage evolution law is energy based as defined by Marigo
:cite:`marigo81a`, :cite:`lemaitre96a`. It is an isotropic damage law.
.. math::
\begin{aligned}
Y &= \frac{1}{2}\boldsymbol{\varepsilon}:\boldsymbol{C}:\boldsymbol{\varepsilon}\\
F &= Y - Y_d - S d\\
d &= \left\{
\begin{array}{l l}
\mathrm{min}\left(\frac{Y-Y_d}{S},\;1\right) & \mathrm{if}\; F > 0\\
\mathrm{unchanged} & \mathrm{otherwise}
\end{array}
\right.\end{aligned}
In this formulation, :math:`Y` is the strain energy release rate,
:math:`Y_d` the rupture criterion and :math:`S` the damage energy. The
non-local version of this damage evolution law is constructed by
averaging the energy :math:`Y`.
.. _sect-smm-cl-damage-mazars:
Mazars
''''''
Keyword: **mazars**
Inherits from **elastic**
Material description with input file:
.. code-block:: python
#input.dat
material mazars [
name = concrete
rho = 3000 # density (Real)
E = 32e9 # young's modulus (Real)
nu = 0.2 # poisson's ratio (Real)
K0 = 9.375e-5 # Damage threshold (Real)
At = 1.15 # Parameter damage traction 1 (Real)
Bt = 10000 # Parameter damage traction 2 (Real)
Ac = 0.8 # Parameter damage compression 1 (Real)
Bc = 1391.3 # Parameter damage compression 2 (Real)
beta = 1.00 # Parameter for shear (Real)
]
Energies:
- ``dissipated``: energy dissipated in damage
-----
This law introduced by Mazars :cite:`mazars84a` is a
behavioral model to represent damage evolution in concrete. This model
does not rely on the computation of the tangent stiffness, the damage is
directly evaluated from the strain.
The governing variable in this damage law is the equivalent strain
:math:`\varepsilon_{{\mathrm{eq}}} =
\sqrt{<\boldsymbol{\varepsilon}>_+:<\boldsymbol{\varepsilon}>_+}`,
with :math:`<.>_+` the positive part of the tensor. This part is defined
in the principal coordinates (I, II, III) as
:math:`\varepsilon_{{\mathrm{eq}}} =
\sqrt{<\boldsymbol{\varepsilon_I}>_+^2 + <\boldsymbol{\varepsilon_{II}}>_+^2 + <\boldsymbol{\varepsilon_{III}}>_+^2}`.
The damage is defined as:
.. math::
\begin{aligned}
D &= \alpha_t^\beta D_t + (1-\alpha_t)^\beta D_c\\
D_t &= 1 - \frac{\kappa_0 (1- A_t)}{\varepsilon_{{\mathrm{eq}}}} - A_t \exp^{-B_t(\varepsilon_{{\mathrm{eq}}}-\kappa_0)}\\
D_c &= 1 - \frac{\kappa_0 (1- A_c)}{\varepsilon_{{\mathrm{eq}}}} - A_c
\exp^{-B_c(\varepsilon_{{\mathrm{eq}}}-\kappa_0)}\\
\alpha_t &= \frac{\sum_{i=1}^3<\varepsilon_i>_+\varepsilon_{{\mathrm{nd}}\;i}}{\varepsilon_{{\mathrm{eq}}}^2}\end{aligned}
With :math:`\kappa_0` the damage threshold, :math:`A_t` and :math:`B_t`
the damage parameter in traction, :math:`A_c` and :math:`B_c` the damage
parameter in compression, :math:`\beta` is the shear parameter.
:math:`\alpha_t` is the coupling parameter between traction and
compression, the :math:`\varepsilon_i` are the eigenstrain and the
:math:`\varepsilon_{{\mathrm{nd}}\;i}` are the eigenvalues of the strain
if the material were undamaged.
The coefficients :math:`A` and :math:`B` are the post-peak asymptotic
value and the decay shape parameters.
.. _sect:smm:CLNL:
Non-Local Constitutive Laws
-```````````````````````````
+'''''''''''''''''''''''''''
Continuum damage modeling of quasi-brittle materials undergo significant
softening after the onset of damage. This fast growth of damage causes a loss of
ellipticity of partial differential equations of equilibrium. Therefore, the
numerical simulation results won't be objective anymore, because the dissipated
energy will depend on mesh size used in the simulation. One way to avoid this
effect is the use of non-local damage formulations. In this approach a local
quantity such as the strain is replaced by its non-local average, where the size
of the domain, over which the quantitiy is averaged, depends on the underlying
material microstructure. ``Akantu`` provides non-local versions of many
constitutive laws for damage. Examples are for instance the material
:ref:`sect-smm-cl-damage-mazars` and the material
:ref:`sect-smm-cl-damage-marigo`, that can be used in a non-local context. In
order to use the corresponding non-local formulation the user has to define the
non-local material he wishes to use in the text input file:
.. code-block:: none
material constitutive_law_non_local [
name = material_name
rho = $value$
...
]
where ``constitutive_law_non_local`` is the name of the non-local constitutive law, *e.g.* `marigo_non_local`.
In addition to the material the non-local neighborhood, that should be used for the averaging process needs to be defined in the material file as well:
.. code-block:: none
non_local neighborhood_name weight_function_type [
radius = $value$
...
weight_function weight_parameter [
damage_limit = $value$
...
]
]
for the non-local averaging, *e.g.* ``base_wf``, followed by the properties of the non-local neighborhood, such as the radius, and the weight function parameters. It is important to notice that the non-local neighborhood must have the same name as the material to which the neighborhood belongs!
The following two sections list the non-local constitutive laws and different type of weight functions available in ``Akantu``.
\subsection{Non-local constitutive laws}
Let us consider a body having a volume :math:`V` and a boundary :math:`\Gamma`. The stress-strain relation for a non-local damage model can be described as follows:
.. _eq:non-local-const:
.. math:: \vec{\sigma} = (1-\bar{d}) \vec{D}:\epsilon
with :math:`\vec{D}` the elastic moduli tensor, :math:`\sigma` the stress tensor, :math:`\epsilon` the strain tensor and :math:`\bar{d}` the non-local damage variable. Note that this stres-strain relationship is similar to the relationship defined in Damage model except :math:`\bar{d}`. The non-local damage model can be extended to the damage constitutive laws: :ref:`sect-smm-cl-damage-marigo` and :ref:`sect-smm-cl-damage-mazars`.
The non-local damage variable :math:`\bar{d}` is defined as follows:
.. _eq:non-local-damage:
.. math:: \bar{d}(\vec{x}) = \int_{V}W(\vec{x}, \vec{y}) d(\vec{y}) dV(\vec{y})
with :math:`W(\vec{x},\vec{y})` the weight function which averages local damage variables to describe the non-local interactions. A list of available weight functions and its functionalities in \akantu are explained in the next section.
Non-local weight functions
-''''''''''''''''''''''''''
+""""""""""""""""""""""""""
The available weight functions in ``Akantu`` are follows:
- ``base_weight_function``: This weight function averages local damage variables by using a bell-shape function on spatial dimensions.
- ``damaged_weight_function``: A linear-shape weight function is applied to average local damage variables. Its slope is determined by damage variables. For example, the damage variables for an element which is highly damaged are averaged over large spatial dimension (linear function including a small slope).
- ``remove_damaged_weight_function``: This weight function averages damage values by using a bell-shape function as ``base_weight_function``, but excludes elements which are fully damaged.
- ``remove_damaged_with_damage_rate_weight_function``: A bell-shape function is applied to average local damage variables for elements having small damage rates.
- ``stress_based_weight_function``: Non local integral takes stress states, and use the states to construct weight function: an ellipsoid shape. Detailed explanations of this weight function are given in Giry et al. :cite:`giry13a`.
-----
.. _sec-cohesive-laws:
Cohesive Constitutive laws
-``````````````````````````
+''''''''''''''''''''''''''
.. _ssect-smm-cl-coh-snozzi:
Linear Irreversible Law
-'''''''''''''''''''''''
+"""""""""""""""""""""""
Keyword: **cohesive_linear**
Material description with input file:
.. code-block:: python
#input.dat
material cohesive_linear [
name = cohesive
sigma_c = 0.1 # critical stress sigma_c (default: 0)
G_c = 1e-2 # Mode I fracture energy
beta = 0. # weighting parameter for sliding and normal opening (default: 0)
penalty = 0. # stiffness in compression to prevent penetration (α in the text)
kappa = 1. # ration between mode-I and mode-II fracture energy (Gc_II/Gc_I)
contact_after_breaking = true # Activation of contact when the elements are fully damaged
max_quad_stress_insertion = false # Insertion of cohesive element when stress is high
# enough just on one quadrature point
# if false the average stress on facet's quadrature points is used
]
-----
.. figure:: figures/cl/linear_cohesive_law.svg
:alt: Irreversible cohesive laws for explicit simulations.
:name: fig:smm:coh:linear_cohesive_law
:align: center
:width: 60.0%
Irreversible cohesive laws for explicit simulations.
`Akantu` includes the Snozzi-Molinari :cite:`snozzi_cohesive_2013`
linear irreversible cohesive law (see
:numref:`fig:smm:coh:linear_cohesive_law`). It is an extension to
the Camacho-Ortiz :cite:`camacho_computational_1996` cohesive law in
order to make dissipated fracture energy path-dependent. The concept
of free potential energy is dropped and a new independent parameter
:math:`\kappa` is introduced:
.. math::
\kappa = \frac{G_\mathrm{c, II}}{G_\mathrm{c, I}}
where :math:`G_\mathrm{c, I}` and :math:`G_\mathrm{c, II}` are the
necessary works of separation per unit area to open completely a
cohesive zone under mode I and mode II, respectively. Their model yields to the
following equation for cohesive tractions :math:`\vec{T}` in case of crack
opening :math:`{\delta}`:
.. math::
\vec{T} = \left( \frac{\beta^2}{\kappa} \Delta_\mathrm{t} \vec{t} +
\Delta_\mathrm{n} \vec{n} \right)
\frac{\sigma_\mathrm{c}}{\delta}
\left( 1- \frac{\delta}{\delta_\mathrm{c}} \right)
= \hat{\vec T}\,
\frac{\sigma_\mathrm{c}}{\delta}
\left( 1- \frac{\delta}{\delta_\mathrm{c}} \right)
:label: eq-smm-coh-tractions
where :math:`\sigma_\mathrm{c}` is the material strength along the fracture,
:math:`\delta_\mathrm{c}` the critical effective displacement after which
cohesive tractions are zero (complete decohesion), :math:`\Delta_\mathrm{t}`
and :math:`\Delta_\mathrm{n}` are the tangential and normal components of
the opening displacement vector :math:`\vec{\Delta}`, respectively. The
parameter :math:`\beta` is a weight that indicates how big the tangential
opening contribution is. The effective opening displacement is:
.. math::
\delta = \sqrt{\frac{\beta^2}{\kappa^2} \Delta_\mathrm{t}^2 + \Delta_\mathrm{n}^2}
In case of unloading or reloading :math:`\delta < \delta_\mathrm{max}`,
tractions are calculated as:
.. math::
\begin{eqnarray}
T_\mathrm{n} &= \Delta_\mathrm{n}\, \frac{\sigma_\mathrm{c}}{\delta_\mathrm{max}} \left( 1- \frac{\delta_\mathrm{max}}{\delta_\mathrm{c}} \right) \\
T_\mathrm{t} &= \frac{\beta^2}{\kappa}\, \Delta_\mathrm{t}\, \frac{\sigma_\mathrm{c}}{\delta_\mathrm{max}} \left( 1- \frac{\delta_\mathrm{max}}{\delta_\mathrm{c}} \right)
\end{eqnarray}
so that they vary linearly between the origin and the maximum attained
tractions. As shown in :numref:`fig:smm:coh:linear_cohesive_law`,
in this law, the dissipated and reversible energies are:
.. math::
\begin{eqnarray}
E_\mathrm{diss} &= \frac{1}{2} \sigma_\mathrm{c}\, \delta_\mathrm{max}\\[1ex]
E_\mathrm{rev} &= \frac{1}{2} T\, \delta
\end{eqnarray}
Moreover, a damage parameter :math:`D` can be defined as:
.. math::
D = \min \left(
\frac{\delta_\mathrm{max}}{\delta_\mathrm{c}},1 \right)
which varies from 0 (undamaged condition) and 1 (fully
damaged condition). This variable can only increase because damage is
an irreversible process. A simple penalty contact model has been incorporated
in the cohesive law so that normal tractions can be returned in
case of compression:
.. math::
T_\mathrm{n} = \alpha \Delta_\mathrm{n} \quad\text{if}\quad
\Delta_\mathrm{n}\quad <\quad 0
where :math:`\alpha` is a stiffness parameter that defaults to zero. The
relative contact energy is equivalent to reversible energy but in
compression.
The material name of the linear decreasing cohesive law is
``material_cohesive_linear`` and its parameters with their respective default
values are:
- ``sigma_c = 0``
- ``delta_c = 0``
- ``beta = 0``
- ``G_c = 0``
- ``kappa = 1``
- ``penalty = 0``
where ``G_c`` corresponds to :math:`G_\mathrm{c, I}`. A random number
generator can be used to assign a random :math:`\sigma_\mathrm{c}` to each
facet following a given distribution (see
Section :ref:`sect-smm-cl`). Only one parameter between ``delta_c``
and ``G_c`` has to be specified. For random :math:`\sigma_\mathrm{c}`
distributions, the chosen parameter of these two is kept fixed and the
other one is varied.
The bi-linear constitutive law works exactly the same way as the linear
one, except for the additional parameter ``delta_0`` that by
default is zero. Two examples for the extrinsic and intrinsic cohesive
elements and also an example to assign different properties to
inter-granular and trans-granular cohesive elements can be found in
the folder ``examples/cohesive_element/``.
----
.. _ssect:smm:cl:coh-friction:
Linear Cohesive Law with Friction
-'''''''''''''''''''''''''''''''''
+"""""""""""""""""""""""""""""""""
Keyword: **cohesive_linear_friction**
Material description with input file:
.. code-block:: python
#input.dat
material cohesive_linear_friction [
name = interface
beta = 1 # weighting parameter for sliding and normal opening (default: 0)
G_c = 30e-3 # Mode I fracture energy
penalty = 1.0e6 # stiffness in compression to prevent penetration (α in the text)
sigma_c = 2.0 # critical stress sigma_c (default: 0)
contact_after_breaking = true # Activation of contact when the elements are fully damaged
mu = 0.5 # Maximum value of the friction coefficient
penalty_for_friction = 5.0e3 # Penalty parameter for the friction behavior
]
-----
This law represents a variation of the linear irreversible cohesive of
the previous section, which adds friction. The friction behavior is
approximated with an elasto-plastic law, which relates the friction
force to the relative sliding between the two faces of the cohesive
element. The slope of the elastic branch is called
``penalty_for_friction``, and is defined by the user, together
with the friction coefficient, as a material property. The friction
contribution evolves with the damage of the cohesive law: it is null
when the damage is zero, and it becomes maximum when the damage is
equal to one. This is done by defining a current value of the
friction coefficient (mu) that increases linearly with the damage, up
to the value of the friction coefficient defined by the user. The
yielding plateau of the friction law is given by the product of the
current friction coefficient and the local compression stress acting
in the cohesive element. Such an approach is equivalent to a
node-to-node contact friction. Its accuracy is acceptable only for
small displacements.
The material name of the linear cohesive law with friction is
``material_cohesive_linear_friction``. Its additional parameters
with respect to those of the linear cohesive law without friction,
with the respective default values, are:
- ``mu = 0``
- ``penalty_for_friction = 0``
.. _ssect:smm:cl:coh-fatigue:
Linear Cohesive Law with Fatigue
-''''''''''''''''''''''''''''''''
+""""""""""""""""""""""""""""""""
Keyword: **cohesive_linear_fatigue**
Material description with input file:
.. code-block:: python
#input.dat
material cohesive_linear_fatigue [
name = cohesive
sigma_c = 1 # critical stress sigma_c (default: 0)
beta = 1 # weighting parameter for sliding and normal opening (default: 0)
delta_c = 1 # Critical displacement
delta_f = 1 # delta_f (normalization of opening rate to alter reloading stiffness after fatigue)
count_switches = true # Count the opening/closing switches per element
]
-----
This law represents a variation of the linear irreversible cohesive
law of the previous section, that removes the hypothesis of elastic
unloading-reloading cycles. With this law, some energy is dissipated
also during unloading and reloading with hysteresis. The
implementation follows the work of :cite:`nguyen2001`. During the
unloading-reloading cycle, the traction increment is computed as
.. math::
\dot{T} =
\begin{cases}
K^- \, \dot{\delta} & \text{if $\dot{\delta} < 0$} \\
K^+ \, \dot{\delta} & \text{if $\dot{\delta} > 0$} \\
\end{cases}
where :math:`\dot{\delta}` and :math:`\dot{T}` are respectively the effective
opening displacement and the cohesive traction increments with respect
to time, while :math:`K^-` and :math:`K^+` are respectively the unloading and
reloading incremental stiffness. The unloading path is linear and
results in an unloading stiffness
.. math::
K^- = \frac{T_\mathrm{max}}{\delta_\mathrm{max}}
where :math:`T_\mathrm{max}` and :math:`\delta_\mathrm{max}` are the maximum
cohesive traction and the effective opening displacement reached
during the precedent loading phase. The unloading stiffness remains
constant during the unloading phase. On the other hand the reloading
stiffness increment :math:`\dot{K}^+` is calculated as
.. math::
\dot{K}^+ =
\begin{cases}
- K^+ \, \dot{\delta} / \delta_\mathrm{f} & \text{if $\dot{\delta}
> 0$} \\
\left( K^+ - K^- \right) \, \dot{\delta} / \delta_\mathrm{f} &
\text{if $\dot{\delta}$ < $0$}
\end{cases}
where :math:`\delta_\mathrm{f}` is a material parameter (refer
to :cite:`vocialta15` for more details). During unloading the stiffness
:math:`K^+` tends to :math:`K^-`, while during reloading :math:`K^+` gets decreased at
every time step. If the cohesive traction during reloading exceeds the
upper limit given by equation :eq:`eq-smm-coh-tractions`, it is
recomputed following the behavior of the linear decreasing cohesive
law for crack opening.
.. _ssect:smm:cl:coh-exponential:
Exponential Cohesive Law
-'''''''''''''''''''''''''
+""""""""""""""""""""""""
Keyword: **cohesive_exponential**
Material description with input file:
.. code-block:: python
#input.dat
material cohesive_exponential [
name = coh1
sigma_c = 1.5e6 # critical stress sigma_c (default: 0)
beta = 1 # weighting parameter for sliding and normal opening (default: 0)
delta_c = 1e-4 # Critical displacement
exponential_penalty = true # Is contact penalty following the exponential law?
contact_tangent = 1.0 # Ratio of contact tangent over the initial exponential tangent
]
-----
Ortiz and Pandolfi proposed this cohesive law in 1999 :cite:`ortiz1999`. The
traction-opening equation for this law is as follows:
.. math::
T = e \sigma_c \frac{\delta}{\delta_c}e^{-\delta/ \delta_c}
:label: eq:exponential_law
This equation is plotted in :numref:`fig:smm:cl:ecl`. The term
:math:`\partial{\vec{T}}/ \partial{\delta}` after the necessary derivation
can expressed as
.. math::
\frac{\partial{\vec{T}}} {\partial{\delta}} = \hat{\vec{T}} \otimes
\frac {\partial{(T/\delta)}}{\partial{\delta}}
\frac{\hat{\vec{T}}}{\delta}+ \frac{T}{\delta} \left[ \beta^2 \mat{I} +
\left(1-\beta^2\right) \left(\vec{n} \otimes \vec{n}\right)\right]
:label: eq:tangent_cohesive
where
.. math::
\frac{\partial{(T/ \delta)}}{\partial{\delta}} = \left\{\begin{array} {l l}
-e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad \text{if}
\delta \geq \delta_{max}\\
0 & \quad \text{if} \delta < \delta_{max}, \delta_n > 0
\end{array} \right.
As regards the behavior in compression, two options are available:
a contact penalty approach with stiffness following the formulation of
the exponential law and a contact penalty approach with constant
stiffness. In the second case, the stiffness is defined as a function
of the tangent of the exponential law at the origin.
.. figure:: figures/cl/cohesive_exponential.png
:alt: Exponential cohesive law
:name: fig:smm:cl:ecl
:align: center
Exponential cohesive law
diff --git a/doc/dev-doc/manual/contactmechanicsmodel.rst b/doc/dev-doc/manual/contactmechanicsmodel.rst
index 617895736..e937840e5 100644
--- a/doc/dev-doc/manual/contactmechanicsmodel.rst
+++ b/doc/dev-doc/manual/contactmechanicsmodel.rst
@@ -1,472 +1,472 @@
Contact Mechanics Model
-=======================
+-----------------------
The contact mechanics model is a specific implementation of
:cpp:class:`Model ` interface to handle contact between
bodies.
Theory
-------
+``````
.. _fig-contactmechanicsmodel-schematic:
.. figure:: figures/contact_mechanics_schematic.png
:align: center
:width: 90%
Basic notation for the contact between two bodies.
Let us consider two deformable bodies, represented as :math:`\Omega_\alpha`,
:math:`\alpha=1, 2`. The boundary :math:`\Gamma_\alpha` of a body is divided
into three non-intersecting regions : :math:`\Gamma^D_\alpha` with
prescribed displacements, :math:`\Gamma^N_\alpha` with prescribed tractions
and :math:`\Gamma^C_\alpha` where the two bodies :math:`\Omega_1` can potentially
come into contact such that:
.. math::
\Gamma^D_\alpha \cup \Gamma^N_\alpha \cup \Gamma^C_\alpha =
\Gamma_\alpha, \quad \Gamma^D_\alpha \cap \Gamma^N_\alpha \cap
\Gamma^C_\alpha = \emptyset
The motion of the two bodies is described in a fixed spatial frame
defined by orthonormal basis :math:`[\boldsymbol{e}_x,
\boldsymbol{e}_y, \boldsymbol{e}_z]` by mapping
:math:`\mathcal{M}_\alpha^t` in time interval :math:`t \in [0, T]`. In
*reference* configuration :math:`i.e.~ t=0`, the position vector for
an arbitrary point on body :math:`\Omega_\alpha` is represented as
:math:`\boldsymbol{X}_\alpha` and in actual configuration the points
are denoted by small letters for example :math:`\boldsymbol{x}`.
:numref:`fig-contactmechanicsmodel-schematic` shows the motion of two
bodies from the *reference* to actual configuration. During the
motion, the bodies can potentially come in contact along
:math:`\Gamma^C_\alpha` as shown in
:numref:`fig-contactmechanicsmodel-schematic`, termed as *potential
contact zone*. Upon contact, two physical conditions need to be
satisfied :math:`(i)` the two surfaces :math:`\Gamma^C_\alpha` cannot
interpenetrate at any time during motion :math:`(ii)` forces must be
exerted by bodies along the :math:`\Gamma^C_\alpha` to resist the
interpenetration as well as relative motion along the contacting
surfaces. Of the two conditions, the non-penetration of bodies can be
defined by a *gap function* :math:`g` which represents the separation
between the two bodies:
.. math::
g = (\boldsymbol{r} -\boldsymbol{\rho}).\boldsymbol{n}
where :math:`\boldsymbol{r} \in \mathcal{M}[\Gamma^C_1]` is the
position of point :math:`\boldsymbol{X}_1` at time :math:`t` ,
:math:`\boldsymbol{\rho} \in \mathcal{M}[\Gamma^C_2]` is the closest
point projection of :math:`\boldsymbol{r}` and :math:`\boldsymbol{n}`
is the outward normal at :math:`\boldsymbol{\rho}` (see
:numref:`fig-contactmechanicsmodel-schematic` b). To preclude the
interpenetration, the *gap function* is constrained with :math:`g \geq 0`, as it is shown in
\Cref{fig:body-contact}b. When the two bodies eventually come in
contact, the gap vanishes :math:`i.e.~ g=0` which leads to the
development of tractions :math:`\boldsymbol{T}_\alpha` along the
contact interface. Thus, the boundary value problem can be formulated
for the two bodies :math:`\Omega_\alpha` with contact constraints as an
extra boundary condition:
.. math::
\boldsymbol{T}\Big |_{{\Gamma^C}_{\alpha}} = \boldsymbol{T}_\alpha
Solution spaces :math:`\mathcal{U}_\alpha` and weighting spaces
:math:`\mathcal{W}_\alpha` are defined for each body:
.. math::
\mathcal{U}_\alpha = \{ \boldsymbol{u}_\alpha: \Omega_\alpha \to \mathbb{R}^d ~|~ \boldsymbol{u}=\boldsymbol{u}[\boldsymbol{x}_\alpha]~ \forall~\boldsymbol{x}_\alpha \in \Gamma^D_\alpha \}
.. math::
\mathcal{W}_\alpha = \{ \boldsymbol{w}_\alpha: \Omega_\alpha \to
\mathbb{R}^d ~|~ \boldsymbol{w}=0~ \forall~ \boldsymbol{x}_\alpha \in
\Gamma^D_\alpha \}
The variational form for each body :math:`\alpha` is given as
.. math::
\int_{\Omega_\alpha}\boldsymbol{\sigma}[\boldsymbol{u}_\alpha]:\boldsymbol{\epsilon}[\boldsymbol{w}_\alpha]~d\Omega_\alpha
=
\int_{\Omega_\alpha}\boldsymbol{b}_\alpha.\boldsymbol{w}_\alpha~d\Omega_\alpha +
\int_{\Gamma^C_\alpha}\boldsymbol{T}_\alpha.\boldsymbol{w}_\alpha~d\Gamma_\alpha +
\int_{\Gamma^N_\alpha}\boldsymbol{T}_\alpha^D.\boldsymbol{w}_\alpha~d\Gamma^N_\alpha
In equilibrium state, following Newton's :math:`3^{rd}` law, it can be
stated the forces exerted by two bodies along :math:`\Gamma^C_\alpha`
are equal and opposite:
.. math::
\boldsymbol{T}_1~d\Gamma^C_1 = - \boldsymbol{T}_2~d\Gamma^C_2
This allows to replace the two integrals over contact surfaces
$\Gamma^C_\alpha$ with a single integral over any one of the surfaces:
.. math::
\int_\alpha
\boldsymbol{T}_\alpha.\delta\boldsymbol{u}_\alpha~d\Gamma^C_\alpha &=
\int_{\Gamma^C_1}\boldsymbol{T}_1\delta\boldsymbol{u}_1~d\Gamma^C_1 +
\int_{\Gamma^C_2} \boldsymbol{T}_2\delta\boldsymbol{u}_2~d\Gamma^C_2 \\ &=
\int_{\Gamma^C_1} \boldsymbol{T}_1(\delta\boldsymbol{u}_1 -
\delta\boldsymbol{u}_2)~d\Gamma^C_1
The contact traction :math:`\boldsymbol{T}_1` can be decomposed into
its normal and tangential components as:
.. math::
\boldsymbol{T}_\alpha = \boldsymbol{T}^n + \boldsymbol{T}^t
= \sigma_n\boldsymbol{n}+ \boldsymbol{T}^t
where :math:`\boldsymbol{T}^n` is the component along the normal
vector :math:`\boldsymbol{n}` and :math:`\boldsymbol{T}^t` is the
component tangential to :math:`\boldsymbol{n}` developed due to the
friction along the surfaces. Upon contact , the normal pressure
:math:`\sigma_n` associated to :math:`\boldsymbol{T}^n` must be
compressive to preclude the interpenetration of bodies. In general, if
a point is not in contact :math:`g > 0`, then :math:`\sigma_n=0` and
in contact :math:`\sigma_n < 0`. This leads to the non-penetration
condition:
.. math::
\sigma_n g = 0
The above set of conditions care called *Hertz-Signorini-Moreau* or
*Karush-Kuhn-Tucker* condition given as:
.. math::
g \geq 0, \quad \sigma_n \leq 0 \quad g.\sigma_n = 0
The tangential component of the contact traction is defined as :
.. math::
\boldsymbol{T}_t = (\boldsymbol{I} - \boldsymbol{n} \otimes \boldsymbol{n})\boldsymbol{\sigma}
where :math:`\boldsymbol{\sigma}` is the Cauchy stress tensor and
$\boldsymbol{n}$ is the outward normal at
:math:`\boldsymbol{\rho}`. The direction of tangential traction,
:math:`\boldsymbol{s}` is determined by the relative sliding velocity
:math:`\boldsymbol{v}_t` of the point :math:`\boldsymbol{r}` and its
projection point :math:`\boldsymbol{\rho}` in contact and is given as:
.. math::
\boldsymbol{s} = \begin{cases} \dfrac{\boldsymbol{v}_t}{\|
\boldsymbol{v}_t \|}, & \text{ if } \| \boldsymbol{v}_t\| > 0 \\ &
\\ 0, & \text{ if } \| \boldsymbol{v}_t\| = 0 \end{cases}
According to the experimental observations of Amontons
and~\cite{coulomb}, in presence of friction, the interface develops a
frictional strength which governs the sliding between them and thus
constrains the tangential contact. In Akantu implementation, we
restrict to the classical non-associated Coulomb's friction
law~\citep{coulomb} which is widely used in many physical and
engineering applications. The Coulomb's friction law, as a first-order
approximation, states that the frictional strength is proportional to
the normal pressure:
.. math::
\sigma_{fric} = \mu |\sigma_n|
where :math:`\mu` is coefficient of friction between interfaces. If
the tangential traction :math:`||\boldsymbol{T}^t||` developed is
below the frictional strength, the relative tangential sliding is zero
:math:`i.e.~\boldsymbol{v}_t = 0` :
.. math::
||\boldsymbol{T}^t|| < \mu| \sigma_n |, \quad \boldsymbol{v}_t=0
The above equation denotes a *stick state*. As soon as the tangential
stress reaches the frictional strength, the two surfaces start
slipping relative to each other :math:`i.e.~\boldsymbol{v}_t > 0`. The
slipping of the surfaces ensures that the tangential stress does not
exceeds the frictional strength, :math:`||\boldsymbol{\sigma}_t||
-\mu|\sigma_n| = 0`. This definition of *slip state* is defined as:
.. math::
||\boldsymbol{\sigma}_t|| -\mu|\sigma_n| = 0, \quad
||\boldsymbol{v}_t|| > 0
Similar to *Karush-Kuhn-Tucker* condition for normal contact, the
above conditions formulate the necessary conditions for tangential
contact:
.. math::
||\boldsymbol{v}_t|| \geq 0, \quad ||\boldsymbol{T}^t||
-\mu|\sigma_n| = 0, \quad ||\boldsymbol{v}_t||
\Big(||\boldsymbol{T}^t|| -\mu|\sigma_n| \Big) = 0
The above contact and frictional constraints are unilateral in nature
*i.e.* they do not behave symmetrically with respect to *gap function,
g*. This renders the balance of work as a variational inequality:
.. math::
\sum_{\alpha=1}^2
\int_{\Omega_\alpha}\boldsymbol{\sigma}[\boldsymbol{u}_\alpha]:\boldsymbol{\epsilon}[\boldsymbol{w}_\alpha]~d\Omega_\alpha
\geq \sum_{\alpha=1}^2 \Big \{
\int_{\Omega_\alpha}\boldsymbol{b}_\alpha.\boldsymbol{w}_\alpha~d\Omega_\alpha +
\int_{\Gamma^Ct_\alpha}\boldsymbol{T}_\alpha.\boldsymbol{w}_\alpha~d\Gamma_\alpha +
\int_{\Gamma^N_\alpha}\boldsymbol{T}_\alpha^D.\boldsymbol{w}_\alpha~d\Gamma_\alpha\Big
\}
which makes it a non-linear optimization problem. The strategy
employed (optimization techniques) to find the solution of variational
inequality depends on the choice of numerical framework employed to
solve the B.V.P.
To solve the minimization problem, FEM introduces the concept of
active set strategy to overcome the problem. In active set strategy,
it is assumed that at current solution step, the part of potential
contact zone :math:`\Gamma^C_1` that are in contact are known,
:math:`\Gamma^{C\star}_1 \subseteq \Gamma^C_1`. This
is achieved by first allowing the two bodies to interpenetrate and
finding the interpenetrated part of potential contact zone which will
denote the active set. Knowing the active set, transforms the
optimization problem to a variational equality where constraints are
imposed along the active part of contacting interface:
.. math::
\sum_{\alpha=1}^2
\int_{\Omega_\alpha}\boldsymbol{\sigma}[\boldsymbol{u}_\alpha]:\boldsymbol{\epsilon}[\delta\boldsymbol{u}_\alpha]~d\Omega_\alpha
= \sum_{\alpha=1}^2 \Big \{
\int_{\Omega_\alpha}\boldsymbol{b}_\alpha.\delta\boldsymbol{u}_\alpha~d\Omega_\alpha +
\int_{\Gamma^N_\alpha}\boldsymbol{T}_\alpha^D.\delta\boldsymbol{u}_\alpha~d\Gamma_\alpha\Big
\} \\
+ \int_{\Gamma^{C\star}_1
}\boldsymbol{T}_1(\delta\boldsymbol{u}_1 -
\delta\boldsymbol{u}_2)~d\Gamma^{C\star}_1
Thus, the resolution of contact problem in FEM requires two steps:
finding the active set along the contacting interface and then
imposing the contact constraints along the active set only. In the
following section, we describe how to employ the contact detection strategies
implemented within Akantu in order to find the active set and to compute contact forces.
Using the Contact Mechanics Model
----------------------------------
+`````````````````````````````````
The :cpp:class:`ContactMechanicsModel `
object solves the contact problem. An instance of
the class can be created like this::
ContactMechanicsModel contact(mesh, spatial_dimension);
while an existing mesh has been used (see \ref{sect:common:mesh}). To
intialize the model object::
contact.initFull(_analysis_method = _explicit_lumped_mass);
The contact mechanics model contains :cpp:class:`Arrays `:
:cpp:func:`gaps `
contains the nodal interpenetrating value :math:`g` (positive for
interpenetration, zero by default after initialization)
:cpp:func:`normals `
contains the normal vector at the slave nodes (zero by default
after initialization).
In Akantu, the possible contact between surfaces is divided into 3
categories.
- Physical Surfaces - The contact occurs between two pyhsically
defined boundaries/surfaces of a body.
- Cohesive Surfaces - The contact occurs between fracturing surfaces
created using :cpp:class:`SolidMechanicsModelCohesive
`.
- All Surfaces - The contact can occur between physical as well as
cohesive surfaces.
To select the contacting surfaces, one must define a
:cpp:class:`SurfaceSelector ` of one of the
above defined types.
To define contact between Physical surfaces, an instance of
:cpp:class:`PhysicalSurfaceSelector `
is created where the mesh object (see \ref{sect:common:mesh}) is
passed as an argument::
auto && surface_selector = std::make_shared(mesh);
To define contact between cohesive surfaces, an instance of
:cpp:class:`CohesiveSurfaceSelector `
must be created. As the contact occurs between the cohesive facets,
therefore the mesh facet object is passed as an argument::
auto && surface_selector = std::make_shared(mesh.getMeshFacets());
To defind contact between physical and cohesive surfaces, an instance of
:cpp:class:`AllSurfaceSelector `
must be created. As the contact occurs between the cohesive facets,
therefore the mesh facet object is passed as an argument::
auto && surface_selector = std::make_shared(mesh.getMeshFacets());
Once a surface selector is created it must be assigned to the
:cpp:class:`ContactDetector ` class::
contact.getContactDetector().setSurfaceSelector(surface_selector);
Contact detection
'''''''''''''''''
The contact detection algorithm can receive the a few parameters. It is possible to
specify the master/slave surfaces with their string identifier. The geometrical projections
are performed with iterations which can be controlled as a classical optimization problem.
A typical detection configuration is given below:
.. code-block::
contact_detector [
type = explicit
master = contact_bottom
slave = contact_top
projection_tolerance = 1e-10
max_iterations = 100
extension_tolerance = 1e-5
]
Contact resolution
''''''''''''''''''
The contact resolution defines the surface tractions due to contact
both for normal and tangential contact. A typical configuration for a
penalization formulation is as follows:
.. code-block::
contact_resolution penalty_linear [
name = contact_top
mu = 0.0 # friction coefficient
epsilon_n = 4e5 # normal penalization
epsilon_t = 1e5 # friction penalization
is_master_deformable = false # when master is rigid => computational savings
]
Coupling with :cpp:class:`SolidMechanicsModel `
-''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
+````````````````````````````````````````````````````````````````````````````
To couple the
:cpp:class:`ContactMechancisModel`
contact mechanics model with
:cpp:class:`SolidMechanicsModel` a
dedicated coupler class :cpp:class:`CouplerSolidContact` is provided.
When an instance of a coupler class is created, it automatically
creates the instances of solid mechanics model and contact mechanics
model. The two objects can be retrived from the coupler class.
.. code-block:: c++
CouplerSolidContact coupler(mesh);
auto & solid = coupler.getSolidMechanicsModel();
auto & contact = coupler.getContactMechanicsModel();
Simply initializing the coupler initializes the two models.
.. code-block:: c++
coupler.initFull( _analysis_method = _explicit_lumped_mass);
However two set the material selector and the contact detector for the
two models, one must set them using directly the instance of the two
model classes.
.. code-block:: c++
auto && selector = std::make_shared>(
"physical_names",solid);
solid.setMaterialSelector(selector);
.. code-block:: c++
auto && surface_selector = std::make_shared(mesh);
contact.getContactDetector().setSurfaceSelector(surface_selector);
The dumping fields/vectors belonging to the solid mechanics model and
contact mechanics model can directly be set through the coupler
class.
.. code-block:: c++
coupler.setBaseName("contact-explicit-dynamic");
coupler.addDumpFieldVector("displacement");
coupler.addDumpFieldVector("normal_force");
coupler.addDumpFieldVector("external_force");
coupler.addDumpFieldVector("internal_force");
coupler.addDumpField("gaps");
coupler.addDumpField("areas");
coupler.addDumpField("stress");
Finally to solve the two models :cpp:func:`solveStep
` function of coupler class must be
invoked.
.. code-block:: c++
coupler.solveStep();
Coupling with :cpp:class:`SolidMechanicsModelCohesive `
-'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
+`````````````````````````````````````````````````````````````````````````````````````````````
To use the contact mechanics model with cohesive elements, one must use the
:cpp:class:`CouplerSolidCohesiveContact` to
coupler the :cpp:class:`ContactMechancisModel`
with
:cpp:class:`SolidMechanicsModelCohesive`.
The initialization and invocation of the functions are similar to
:cpp:class:`CouplerSolidContact` except a few
changes.
.. code-block:: c++
solid = coupler.getSolidMechanicsModelCohesive();
While initializing the coupler, the nature of cohesive elements
(extrinsic/intrinsic) should need to be passed.
.. code-block:: c++
coupler.initFull( _analysis_method = _explicit_lumped_mass, _is_extrinsic=true);
To ensure that cohesive elements break during an explicit insertion, one must
call the function :cpp:func:`checkCohesiveStress()
` after
:cpp:func:`solveStep() `.
.. code-block:: c++
coupler.solveStep();
solid.checkCohesiveStress();
diff --git a/doc/dev-doc/manual/heattransfermodel.rst b/doc/dev-doc/manual/diffusion.rst
similarity index 98%
copy from doc/dev-doc/manual/heattransfermodel.rst
copy to doc/dev-doc/manual/diffusion.rst
index d08f31f12..ca8341382 100644
--- a/doc/dev-doc/manual/heattransfermodel.rst
+++ b/doc/dev-doc/manual/diffusion.rst
@@ -1,169 +1,169 @@
-Heat Transfer Model
-===================
+Diffusion Model
+---------------
The heat transfer model is a specific implementation of the :cpp:class:`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.
Using the Heat Transfer Model
------------------------------
+`````````````````````````````
A material file name has to be provided during initialization.
Currently, the :cpp:class:`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 `
More precisely, the heat transfer model contains 4 :cpp:class:`Arrays `:
- **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 [
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 ` 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/doc/dev-doc/manual/extra.rst b/doc/dev-doc/manual/extra.rst
new file mode 100644
index 000000000..d46b839f6
--- /dev/null
+++ b/doc/dev-doc/manual/extra.rst
@@ -0,0 +1,6 @@
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/doc/dev-doc/manual/heattransfermodel.rst b/doc/dev-doc/manual/heattransfermodel.rst
index d08f31f12..ca8341382 100644
--- a/doc/dev-doc/manual/heattransfermodel.rst
+++ b/doc/dev-doc/manual/heattransfermodel.rst
@@ -1,169 +1,169 @@
-Heat Transfer Model
-===================
+Diffusion Model
+---------------
The heat transfer model is a specific implementation of the :cpp:class:`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.
Using the Heat Transfer Model
------------------------------
+`````````````````````````````
A material file name has to be provided during initialization.
Currently, the :cpp:class:`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 `
More precisely, the heat transfer model contains 4 :cpp:class:`Arrays `:
- **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 [
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 ` 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/doc/dev-doc/manual/models.rst b/doc/dev-doc/manual/models.rst
new file mode 100644
index 000000000..f9838c28b
--- /dev/null
+++ b/doc/dev-doc/manual/models.rst
@@ -0,0 +1,19 @@
+.. _sect-models:
+
+Models
+======
+
+.. include:: models_commons.rst
+
+.. include:: solidmechanicsmodel.rst
+.. include:: contactmechanicsmodel.rst
+
+Models [Unsable]
+================
+
+This models are implemented for some research purposes but might not have the
+same level of stability or tests than the models that are not marked unsable.
+
+.. include:: diffusion.rst
+.. include:: phasefieldmodel.rst
+.. include:: structuralmechanicsmodel.rst
diff --git a/doc/dev-doc/manual/models_commons.rst b/doc/dev-doc/manual/models_commons.rst
new file mode 100644
index 000000000..e69de29bb
diff --git a/doc/dev-doc/manual/new-constitutive-laws.rst b/doc/dev-doc/manual/new-constitutive-laws.rst
index 300b9f4e3..d7aee1eb3 100644
--- a/doc/dev-doc/manual/new-constitutive-laws.rst
+++ b/doc/dev-doc/manual/new-constitutive-laws.rst
@@ -1,373 +1,373 @@
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 :cpp:class:`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::
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 1);
solver.set("threshold", 1e-4);
solver.set("convergence_type", SolveConvergenceCriteria::_residual);
model.solveStep();
/* ---------------------------------------------------------------------- */
#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 (Int 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(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*.
diff --git a/doc/dev-doc/manual/phasefieldmodel.rst b/doc/dev-doc/manual/phasefieldmodel.rst
index 1935beb16..6dc660ef8 100644
--- a/doc/dev-doc/manual/phasefieldmodel.rst
+++ b/doc/dev-doc/manual/phasefieldmodel.rst
@@ -1,265 +1,265 @@
PhaseField Model
-================
+----------------
The phasefield model is a specific implementation of
:cpp:class:`Model ` interface to handle brittle fracture
for infinitesimal strains.
Theory
-------
+``````
The variational formulation of brittle fracture was first proposed by
Francfort and Marigo:cite:`franc` so as to overcome the shortcomings
of Griffith's criteria. The approach used by Francfort and Margio for
variational formulation is based on the principle of global minimality
of total energy. Similar to Griffith's work, they define a surface
energy corresponding to the discontinuity :math:`\Gamma` as
.. math::
\Phi_{s}(\Gamma) = \int\limits_{\Gamma}\mathcal{G}_cd\mathcal{H}^{d-1}(\Gamma)
where :math:`\mathcal{H}^{d-1}` is d-1 dimensional Hausdroff measure
which is surface measure for smooth hypersurfaces. the material
behavior is considered to be linearly elastic with small strains
throughout the body. Based on these assumptions the elastic energy of
the body is defined as
.. math::
\Phi_e(\Gamma,u) = \int\limits_{\Omega_{p}}\psi(\boldsymbol{\epsilon}(u))d\Omega
where :math:`\psi` is elastic energy density and it is given as a
function of strain as
.. math::
\psi_e(\boldsymbol{\epsilon}) = \dfrac{1}{2}\lambda(\text{tr}(\boldsymbol{\epsilon}))^{2}+\mu\boldsymbol{\epsilon}:\boldsymbol{\epsilon}
Exponential phasefield law
''''''''''''''''''''''''''
Bourdin et. al.:cite:`bourdin` in their work proposed a regularized
version of variational formulation where a scalar field variable
:math:`d(\vec{x},t)` is used to represent crack. This scalar field
variable approximates the sharp crack topology by taking value 1 at
crack location and smoothly diffusing into value 0 away from the
crack. \figref{fig:diffusive_topology} shows both sharp crack topology
and an approximated regularized crack topology for one-dimensional
case.
The strong form of the phasefield can be expressed as
.. math::
2(1-d)\mathcal{H} - \dfrac{G_c}{l_0}(d-l_0^2d) = 0
\nabla d . \boldsymbol n = 0
Using the PhaseField Model
---------------------------
+``````````````````````````
The :cpp:class:`PhaseFieldModel ` object
solves the strong form using an inplicit solver. An instance of the
class can be created like this::
PhaseFieldModel phase(mesh, spatial_dimension);
while ans existing mesh has been used (see \ref{sect:common:mesh}). To
intialize the model object::
phase.initFull();
Currently, implicit solver is defined for the phasefield model and no
explicit solver is implemented in Akantu. Furthermore, By default, the
implicit solver defined is of linear type. One can change the linear
solver to a non-linear by initiating a new solver type.
The phasefield model contains :cpp:class:`Arrays `:
:cpp:func:`damage `:
contains the nodal damage :math:`d` (zero by default after the
initialization)
:cpp:func:`blocked_dofs `
contains a Boolean value specifying whether the damage at a node
is to be blocked or not. A Dirichlet boundary condition can be
prescribed by setting the **blocked_dofs** value of damage to
``true``. The **damage** ais computed for all nodes where the
**blocked_dofs** value is set to ``false``. For the remaining
nodes, the imposed values (zero by default after initialization)
are kept.
:cpp:func:`internal_force `
contains the driving force responsible for the crack to nucelate
or propagate.
Currently, the phasefield model uses a exponential shaped scalar
field variable :math:`d(\vec{x}, t)` to approximate the sharp crack
topology. A Phasefield variable thus requires a length scale
parameter :math:`l_0` to control the width of the diffusive crack
topology along with the material parameters such as elastic modulus,
poisson's ratio, critical energy release rate.
The data input file provides the parameters for the exponential
phasefield law as follows:
.. code-block::
phasefield exponential [
name = plate
E = 210.0
nu = 0.3
gc = 5e-3
l0 = 0.1
]
:cpp:class:`PhaseFieldModel` can handle
phasefield laws for multiple materials. To define so:
.. code-block::
phasefield exponential [
name = hard
E = 210.0
nu = 0.3
gc = 5e-3
l0 = 0.1
]
phasefield exponential [
name = soft
E = 21.0
nu = 0.3
gc = 5e-5
l0 = 0.1
]
In order to assign correct phasefield variable properties based on the
names of region as defined in mesh file,
:cpp:class:`MeshDataPhaseFieldSelector
` must be set as phasefield
selector::
auto && selector = std::make_shared>(
"physical_names", phase);
phase.setPhaseFieldSelector(selector);
For the crack to nucleate or propagate, phasefield model requires
strain measure at each quadrature point. The strain measure computed
from solid mechanics model is provided to the phasefield
model. Similarly, damage value computed at each quadrature point by
the phasefield model is provided to the solid mechanics model. This
damage is thus used to degrade the the elastic strain energy. To do
so, a new material which is a specific implementation of
:cpp:class:`MaterialDamage ` is defined for
solid mechanics model. The governing equation is given as
.. math::
\boldsymbol \sigma = ((1-d)^2 + \eta)\dfrac{\psi}{\boldsymbol \epsilon}
where :math:`\psi` is the elastic energy and :math:`\eta` is a
numerical parmaeter avoid numerical difficulties due to full
degradation of elastic energy for fully broken state. The material
properties are thus provided in the input data filled as:
.. code-block::
material phasefield [
name = hard
rho = 1.
E = 210.0
nu = 0.3
eta = 0.0
Plane_Stress = false
]
material phasefield [
name = soft
rho = 1.
E = 21.0
nu = 0.3
eta = 0.0
Plane_Stress = false
]
To simplify the execution of phasefield model coupled with
solidmechanis model, a special class
:cpp:class:`CouplerSolidPhaseField` is
provided.
Coupling Phase Field Model and Solid Mechanics Model
''''''''''''''''''''''''''''''''''''''''''''''''''''
A dedicated coupler class :cpp:class:`CouplerSolidPhaseField
` is defined in Akantu to ease the
coupling of the two models.
When an instance of a Coupler class is created, it automatically
creates the instances of solid mechanics model and phasefield
model. The two objects can be retrived from the coupler class.
.. code-block:: c++
CouplerSolidPhaseField coupler(mesh);
auto & phase = coupler.getPhaseFieldModel();
auto & solid = coupler.getSolidMechanicsModel();
The two objects must be used to define the solver type and apply
boundary conditions.
.. code-block:: c++
solid.initFull(_analysis_method = _explicit_lumped_mass);
phase.initFull(_analysis_method = _static);
The whole process of coupling the two models at a given time step is
made easy by the :cpp:func:`solve `
function of coupler class.
.. code-block:: c++
coupler.solve(, );
Staggered scheme
-----------------
+````````````````
To solve the solid mechanics model and phasefield model, staggered
scheme is implemented. In staggered scheme, at current time step first
solid mechanics model is solved assuming the damage values from the
previous time step. The strain thus computed are passed to the
phasefield model and now, phasefield model is solved for the damage
variable. At each iteration step, convergence in displacement and
damage is checked for. If convergence is not reached, the combined
newton-raphson iteration continues.
To illustrate the staggered scheme solution of phasefield model, one
can refere to the examples provided for the phasefield model. Two
examples provided are for static `phasefield-static.py` and dynamic
crack propgation `phasefield-dynamic.py`.
For the static problem, both the solid mechanics model and the
phsefield model are solved using alinear implicit solver. Convergence
in value of both displacement as well as damage is checked at each
loading step.
In case of the dynamic problem, solid mechanics model is solved
dynamically using an explicit solver and the phasefield model is
solved using a linear implcit solver. In this sceanrio, the staggered
scheme doesnot check for any convergence. Below is the crack
propagation observed for the dyanmic problem.
.. _fig-phasefieldmodel-dynamic:
.. figure:: figures/phasefield-dynamic.png
:align: center
Dynamic crack propagation using phasefield model
.. LocalWords: phasefield SolidMechanics PhaseField akantu cpp
diff --git a/doc/dev-doc/manual/solidmechanicsmodel.rst b/doc/dev-doc/manual/solidmechanicsmodel.rst
index 6b214a866..7480b162e 100644
--- a/doc/dev-doc/manual/solidmechanicsmodel.rst
+++ b/doc/dev-doc/manual/solidmechanicsmodel.rst
@@ -1,886 +1,886 @@
.. _sect-smm:
Solid Mechanics Model
-=====================
+---------------------
The solid mechanics model is a specific implementation of the :cpp:class:`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 :cpp:class:`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``:
:cpp:func:`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.
:cpp:func:`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).
:cpp:func:`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).
:cpp:func:`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).
:cpp:func:`external_force `
contains the external forces applied on the nodes
(:math:`\vec{f}_{\st{ext}}` in the following).
:cpp:func:`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. In addition to vector quantities, the solid
mechanics model can be queried for energies with the :cpp:func:`getEnergy
` which accepts an energy type as an
arguement (e.g. ``kinetic`` or ``potential``).
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 (Int node = 0; node < mesh.getNbNodes(); ++node) {
disp(node, 0) = 0.1;
velo(node, 1) = 1.;
}
.. _sect-smm-boundary:
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
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
.. code-block:: c++
auto & blocked = model.getBlockedDOFs();
const auto & pos = mesh.getNodes();
UInt nb_nodes = mesh.getNbNodes();
for (Int 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``
.. code-block:: c++
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 :cpp:func:`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:
:cpp:class:`FlagOnly `, :cpp:class:`FixedValue
` and :cpp:class:`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
.. code-block:: c++
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
:cpp:class:`FromTraction ` of Neumann
boundary conditions is called. Otherwise, the functor :cpp:class:`FromStress
` should be called which gets the stress tensor
as an input parameter
.. code-block:: c++
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
.. code-block:: c++
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 :numref:`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.
:cpp:class:`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``.
:cpp:class:`MaterialSelector ` and
:cpp:class:`DefaultMaterialSelector ` is used
to assign a material to regular elements by default. For the regular elements,
as in the example above, :cpp:class:`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.
For cohesive material, ``Akantu`` has a pre-defined material selector to assign
the first cohesive material by default to the cohesive elements which is called
:cpp:class:`DefaultMaterialCohesiveSelector
` and it inherits its properties from
:cpp:class:`DefaultMaterialSelector `. Multiple
cohesive materials can be assigned using mesh data information (for more
details, see :ref:`sect-smm-intrinsic-insertion`).
Insertion of Cohesive Elements
-``````````````````````````````
+''''''''''''''''''''''''''''''
Cohesive elements are currently compatible only with static simulation
and dynamic simulation with an explicit time integration scheme (see
section :ref:`ssect-smm-expl-time-integration`). They do not have to be
inserted when the mesh is generated (intrinsic) but can be added
during the simulation (extrinsic). At any time during the simulation,
it is possible to access the following energies with the relative
function:
.. code-block:: c++
Real Ed = model.getEnergy("dissipated");
Real Er = model.getEnergy("reversible");
Real Ec = model.getEnergy("contact");
A new model have to be call in a very similar way that the solid
mechanics model:
.. code-block:: c++
SolidMechanicsModelCohesive model(mesh);
model.initFull(_analysis_method = _explicit_lumped_mass,
_is_extrinsic = true);
Cohesive element insertion can be either realized at the beginning of
the simulation or it can be carried out dynamically during the
simulation. The first approach is called ``intrinsic``, the second
one ``extrinsic``. When an element is present from the beginning, a
bi-linear or exponential cohesive law should be used instead of a
linear one. A bi-linear law works exactly like a linear one except for
an additional parameter :math:`\delta_0` separating an initial linear
elastic part from the linear irreversible one. For additional details
concerning cohesive laws see Section~\ref{sec:cohesive-laws}.
.. _fig-smm-coh-insertion:
.. figure:: figures/insertion.svg
:align: center
Insertion of a cohesive element.
Extrinsic cohesive elements are dynamically inserted between two
standard elements when
.. math::
\sigma_\st{eff} > \sigma_\st{c} \quad\text {with} \quad \sigma_\st{eff} = \sqrt{\sigma_\st{n} ^ 2 + \frac{\tau ^ 2} {\beta ^ 2 }}
in which :math:`\sigma_\st{n}` is the tensile normal traction and $\tau$ the
resulting tangential one (:numref:`fig-smm-coh-insertion`).
Extrinsic approach
-''''''''''''''''''
+""""""""""""""""""
During the simulation, stress has to be checked along each facet in order to
insert cohesive elements where the stress criterion is reached. This check is
performed by calling the method :cpp:func:`checkCohesiveStress
`, as example before
each step resolution:
.. code-block:: c++
model.checkCohesiveStress();
model.solveStep();
The area where stresses are checked and cohesive elements inserted can be
limited using the method :cpp:func:`setLimit
` on the :cpp:class:`CohesiveInserter
` during initialization. As example, to limit
insertion in the range :math:`[-1.5, 1.5]` in the $x$ direction:
.. code-block:: c++
auto & inserter = model.getElementInserter();
inserter.setLimit(_x, -1.5, 1.5);
Additional restrictions with respect to :math:`_y` and :math:`_z` directions can
be added as well.
.. _sect-smm-intrinsic-insertion:
Intrinsic approach
-''''''''''''''''''
+""""""""""""""""""
Intrinsic cohesive elements are inserted in the mesh with the method
:cpp:func:`initFull `.
Similarly, the range of insertion can be limited with :cpp:func:`setLimit
` before the :cpp:func:`initFull
` call.
In both cases extrinsic and intrinsic the insertion can be restricted to
surfaces or element groups. To do so the list of groups should be specified in
the input file.
.. code-block::
model solid_mechanics_model_cohesive [
cohesive_inserter [
cohesive_surfaces = [surface1, surface2, ...]
cohesive_zones = [group1, group2, ...]
]
material cohesive_linear [
name = insertion
beta = 1
G_c = 10
sigma_c = 1e6
]
]
``cohesive_surfaces`` defines are the physical surfaces defined in the mesh
(``Curves`` in 1D and ``Surfaces`` in 2D) on which the insertion is allowed,
``cohesive_zones`` are the physical volumes in which the insertion is allowed
(``Surfaces`` in 2D and ``Volumes`` in 3D).
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}}
:label: eqn-smm-static
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 :cpp:func:`initFull
` method by giving the
:cpp:class:`SolidMechanicsModelOptions `.
These options specify the type of analysis to be performed and whether the
materials should be initialized with :cpp:func:`initMaterials
` or not
.. code-block:: c++
SolidMechanicsModel model(mesh);
model.initFull(_analysis_method = _static);
Here, a static analysis is chosen by passing the argument
:cpp:enumerator:`_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
:cpp:func:`solveStep ` can be called
.. code-block:: c++
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 1);
solver.set("threshold", 1e-4);
solver.set("convergence_type", SolveConvergenceCriteria::_residual);
model.solveStep();
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 (:eq:`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 (:eq:`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
(:cpp:enumerator:`SolveConvergenceCriteria::_residual
`), i.e. :math:`||\vec{r}|| <`
:cpp:member:`threshold
`. One can also
iterate until the increment of displacement is zero
(:cpp:enumerator:`SolveConvergenceCriteria::_solution
`) 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`
.. code-block:: c++
solver.set("max_iterations", 100);
model.solveStep();
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`.
.. _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).
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:
.. math:: \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = \vec{f}_{\mathrm{ext}}
:label: eqn-equation-motion
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, (:eq:`eqn-equation-motion`) becomes a
system of three equations (see :cite:`curnier92a,hughes-83a` for more details):
.. math::
\begin{eqnarray}
\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}\\
\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}\\
\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}
\end{eqnarray}
:label: eqn-equation-motion-discret
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~(:eq:`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`.
.. csv-table::
:header: ":math:`\\alpha`", "Method (:math:`\\beta = 1/2`)", "Type"
":math:`0`", "central difference", "explicit"
":math:`\frac{1}{6}`", "Fox-Goodwin(royal road)", "implicit"
":math:`\frac{1}{3}`", "Linear acceleration", "implicit"
":math:`\frac{1}{2}`", "Average acceleration (trapeziodal rule)", "implicit"
The solution of this system of equations,
(:eq:`eqn-equation-motion-discret`) 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
.. csv-table::
:header: "", ":math:`\\vec{w}`", ":math:`e`", ":math:`d`", ":math:`c`"
"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}`"
.. 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 (Int 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 (Int 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.
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: :math:`7800\mathrm{kg/m}^3`, Young's modulus:
:math:`210\mathrm{GPa}` 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`.
.. include:: ./constitutive-laws.rst
.. include:: ./new-constitutive-laws.rst
diff --git a/doc/dev-doc/manual/structuralmechanicsmodel.rst b/doc/dev-doc/manual/structuralmechanicsmodel.rst
index 1da5dd206..d1087c70d 100644
--- a/doc/dev-doc/manual/structuralmechanicsmodel.rst
+++ b/doc/dev-doc/manual/structuralmechanicsmodel.rst
@@ -1,214 +1,214 @@
Structural Mechanics Model
-==========================
+--------------------------
Static structural mechanics problems can be handled using the class
:cpp:class:`StructuralMechanicsModel `. So
far, ``Akantu`` provides 2D and 3D Bernoulli beam elements :cite:`frey2009`.
This model is instantiated for a given :cpp:class:`Mesh `, as for
the :cpp:class:`SolidMechanicsModel `. The model
will create its own :cpp:class:`FEEngine ` object to compute
the interpolation, gradient, integration and assembly operations. The
:cpp:class:`StructuralMechanicsModel `
constructor is called in the following way:
.. code-block:: c++
StructuralMechanicsModel model(mesh, spatial_dimension);
where ``mesh`` is a :cpp:class:`Mesh ` object defining the structure for
which the equations of statics are to be solved, and
``spatial_dimension`` is the dimensionality of the problem. If
``spatial_dimension`` is omitted, the problem is assumed to have
the same dimensionality as the one specified by the mesh.
.. warning::
Dynamic computations are not supported to date.
.. note::
Structural meshes are created and loaded
with ``_miot_gmsh_struct`` instead of ``_miot_gmsh`` (cf. :ref:`loading_mesh`)
.. code-block:: c++
Mesh mesh;
mesh.read("structural_mesh.msh", _miot_gmsh_struct);
This model contains at least the following :cpp:class:`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``. The **displacement** is 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_rotation** contains the generalized displacements (*i.e.* displacements and rotations) 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).
- **external_force** contains the generalized external forces (forces and moments) applied to the nodes (:math:`\vec{f_{\st{ext}}}` in the following).
- **internal_force** contains the generalized internal forces (forces and moments) applied to the nodes (:math:`\vec{f_{\st{int}}}` in the following).
An example to help understand how to use this model will be presented in the
next section.
.. _sec:structMechMod:setup:
Model Setup
------------
+```````````
Initialization
-``````````````
+''''''''''''''
The easiest way to initialize the structural mechanics model is:
.. code-block:: c++
model.initFull();
The method :cpp:class:`initFull ` computes the shape
functions, initializes the internal vectors mentioned above and allocates the
memory for the stiffness matrix, unlike the solid mechanics model, its default
argument is ``_static``.
Material properties are defined using the :cpp:class:`StructuralMaterial
` structure described in
:numref:`tab-structmechmod-strucmaterial`. Such a definition could, for
instance, look like
.. code-block:: c++
StructuralMaterial mat1;
mat.E=3e10;
mat.I=0.0025;
mat.A=0.01;
.. _tab-structmechmod-strucmaterial:
.. table:: Material properties for structural elements defined in the class :cpp:class:`StructuralMaterial `.
:align: center
====== ======
Field Description
====== ======
``E`` Young's modulus
``A`` Cross section area
``I`` Second cross sectional moment of inertia (for 2D elements)
``Iy`` ``I`` around beam :math:`y`--axis (for 3D elements)
``Iz`` ``I`` around beam :math:`z`--axis (for 3D elements)
``GJ`` Polar moment of inertia of beam cross section (for 3D elements)
====== ======
Materials can be added to the model's ``element_material`` vector using
.. code-block:: c++
model.addMaterial(mat1);
They are successively numbered and then assigned to specific elements.
.. code-block:: c++
for (Int i = 0; i < nb_element_mat_1; ++i) {
model.getElementMaterial(_bernoulli_beam_2)(i,0) = 1;
}
.. _sect:structMechMod:boundary:
Setting Boundary Conditions
-```````````````````````````
+'''''''''''''''''''''''''''
As explained before, the Dirichlet boundary conditions are applied through the
array **blocked_dofs**. Two options exist to define Neumann conditions.
If a nodal force is applied, it has to be directly set in the array
**force_momentum**. For loads distributed along the beam length, the
method :cpp:class:`computeForcesFromFunction ` integrates them into nodal forces. The
method takes as input a function describing the distribution of loads along the
beam and a functor :cpp:class:`BoundaryFunctionType ` specifing if the function is expressed in the local coordinates (``_bft_traction_local``) or in the
global system of coordinates (``_bft_traction``).
.. code-block:: c++
static void lin_load(double * position, double * load,
Real * normal, UInt surface_id){
memset(load,0,sizeof(Real)*3);
load[1] = position[0]*position[0]-250;
}
int main(){
...
model.computeForcesFromFunction<_bernoulli_beam_2>(lin_load,
_bft_traction_local);
...
}
.. _sect:structMechMod:static:
Static Analysis
----------------
+```````````````
The :cpp:class:`StructuralMechanicsModel ` class can perform static analyses of structures. In this case, the equation to solve is the same as for the :cpp:class:`SolidMechanicsModel ` used for static analyses
.. math:: \mat{K} \vec{u} = \vec{f_{\st{ext}}}~,
:label: eqn-structmechmod-static
where :math:`\mat{K}` is the global stiffness matrix, :math:`\vec{u}` the
generalized displacement vector and :math:`\vec{f_{\st{ext}}}` the vector of
generalized external forces applied to the system.
To solve such a problem, the static solver of the
:cpp:class:`StructuralMechanicsModel ` object
is used. First a model has to be created and initialized.
.. code-block:: c++
StructuralMechanicsModel model(mesh);
model.initFull();
- :cpp:func:`model.initFull ` initializes all
internal vectors to zero.
Once the model is created and initialized, the boundary conditions can be set as explained in Section :ref:`sect:structMechMod:boundary`. Boundary conditions will prescribe the external forces or moments for the free degrees of freedom :math:`\vec{f_{\st{ext}}}` and displacements or rotations for the others. To completely define the system represented by equation (:eq:`eqn-structmechmod-static`), the global stiffness matrix :math:`\mat{K}` must be assembled.
.. code-block:: c++
model.assembleStiffnessMatrix();
The computation of the static equilibrium is performed using the same
Newton-Raphson algorithm as described in
Section~\ref{sect:smm:static}.
\note{To date, :cpp:class:`StructuralMechanicsModel
` handles only constitutively and
geometrically linear problems, the algorithm is therefore guaranteed to converge
in two iterations.}
.. code-block:: c++
model.solveStep();
- :cpp:func:`model.solveStep ` solves the :eq:`eqn-structmechmod-static`.
The **increment** vector of the model will contain the new
increment of displacements, and the **displacement_rotation**
vector is also updated to the new displacements.
At the end of the analysis, the final solution is stored in the
**displacement_rotation** vector. A full example of how to solve a structural
mechanics problem is presented in the code
``example/structural_mechanics/bernoulli_beam_2_example.cc``. This example is
composed of a 2D beam, clamped at the left end and supported by two rollers as
shown in :numref:`fig-structmechmod-exam1-1`. The problem is defined by the
applied load :math:`q=6 \text{\kN/m}`, moment :math:`\bar{M} = 3.6 \text{kN m}`,
moments of inertia :math:`I_1 = 250\,000 \text{cm}^4` and :math:`I_2 = 128\,000
\text{cm}^4` and lengths :math:`L_1 = 10\text{m}` and :math:`L_2 = 8\text{m}`.
The resulting rotations at node two and three are :math:`\varphi_2 = 0.001\,167`
and :math:`\varphi_3 = -0.000\,771`.
.. _fig-structmechmod-exam1-1:
.. figure:: figures/beam_example.svg
:align: center
2D beam example
diff --git a/doc/dev-doc/requirements.txt b/doc/dev-doc/requirements.txt
index dd722034f..da3a21c59 100644
--- a/doc/dev-doc/requirements.txt
+++ b/doc/dev-doc/requirements.txt
@@ -1,8 +1,9 @@
urllib3<2
sphinx==5.3.0
sphinx-rtd-theme==1.2.0
sphinxcontrib-bibtex==2.5.0
breathe==4.34.0
jinja2==3.1.2
gitpython==3.1.30
myst-parser==0.18.1
+sphinx-copybutton==0.4.0
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index cc1bdad9c..cd3756c11 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,43 +1,31 @@
#===============================================================================
# Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This file is part of Akantu
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with Akantu. If not, see .
#
#===============================================================================
-add_example(new_material "Example on how to add a new material in Akantu" PACKAGE core)
-add_example(boundary_conditions "Example on hoy to set boundary conditions" PACKAGE core)
-add_example(explicit "Example on how to run an explicit simulation" PACKAGE core)
-add_example(io "Example on how to perform Input/Output operations" PACKAGE core)
-add_example(implicit "Example on how to run an implicit simulation" PACKAGE implicit)
-add_example(static "Example on how to run a static simulation" PACKAGE implicit)
-add_example(parallel "Example of how to write a parallel code with Akantu" PACKAGE parallel)
-add_example(cohesive_element "Cohesive element examples" PACKAGE cohesive_element)
-add_example(structural_mechanics "Structural mechanics model examples" PACKAGE structural_mechanics)
-add_example(heat_transfer "Example on how to run heat transfer simulation" PACKAGE heat_transfer)
-add_example(python "Example on how to use the python interface" PACKAGE python_interface)
-add_example(embedded "Example on how to run embedded model simulation" PACKAGE embedded)
-add_example(contact_mechanics "Example on how to run contact mechanics model simulation" PACKAGE contact_mechanics)
-add_example(phase_field "Example on how to run phase field model simulation" PACKAGE phase_field)
+add_example(c++ "C++ examples" PACKAGE core)
+add_example(python "Python example" PACKAGE python_interface)
package_add_files_to_package(
examples/README.rst
cmake/AkantuExampleMacros.cmake
)
diff --git a/examples/README.rst b/examples/README.rst
index e317e4b13..64b5ab67d 100644
--- a/examples/README.rst
+++ b/examples/README.rst
@@ -1,59 +1,4 @@
-========
Examples
========
-This examples demonstrate some of the principal features of **Akantu**. Each
-sub-folder contain concentrate on one particular features.
-
-Common features
----------------
-
-*boundary_conditions*
- This shows how to impose boundary conditions by using the group of elements
- defining the boundaries. Boundaries are applied mainly by using functors, some
- are predefined in akantu asshown in the sub-example 'predefined_bc'. If this
- is not sufficient the user can define it's one functor, example
- 'user_defined_bc'
-
-*io*
- This examples show some advanced features of the dumpers, to dump extra
- fields than the default ones
-
-*parallel*
- This example shows how to initialize a parallel computation
-
-*python*
- This examples show some basic usage of the python interface
-
-Solid Mechanics
----------------
-
-*explicit*
- This examples shows a dynamic wave propagation in a beam. It uses a central
- difference explicit time integration scheme.
-
-*implicit*
- This example shows how to use the implicit dynamic solver.
-
-*static*
- This example shows how to do a static computation
-
-*cohesive_element*
- This examples show some usage of the cohesive element in explicit or implicit
- and with extrinsic or intrinsic elements
-
-*new_material*
- This example shows how to write a custom material outside of **Akantu** and
- how to use it
-
-Heat transfer
--------------
-
-*heat_transfer*
- This example shows how to use the HeatTransferModel
-
-Structural mechanics
---------------------
-
-*structural_mechanics*
- This example shows how to use the StructuralMechanicsModel
+Akantu's example are separated in 2 types, the C++ and the python ones.
diff --git a/examples/c++/CMakeLists.txt b/examples/c++/CMakeLists.txt
new file mode 100644
index 000000000..b643bbad3
--- /dev/null
+++ b/examples/c++/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_example(contact_mechanics_model "Example on how to run contact mechanics model simulation" PACKAGE contact_mechanics)
+add_example(diffusion_model "Example on how to run heat transfer simulation" PACKAGE diffusion)
+add_example(embedded_model "Example on how to run embedded model simulation" PACKAGE embedded)
+add_example(phase_field_model "Example on how to run phase field model simulation" PACKAGE phase_field)
+add_example(solid_mechanics_cohesive_model "Cohesive element examples" PACKAGE cohesive_element)
+add_example(solid_mechanics_model "Solid mechanics examples" PACKAGE solid_mechanics)
+add_example(structural_mechanics_model "Structural mechanics model examples" PACKAGE structural_mechanics)
diff --git a/examples/c++/README.rst b/examples/c++/README.rst
new file mode 100644
index 000000000..db3ca35e5
--- /dev/null
+++ b/examples/c++/README.rst
@@ -0,0 +1,2 @@
+C++ Examples
+------------
diff --git a/examples/contact_mechanics/CMakeLists.txt b/examples/c++/contact_mechanics_model/CMakeLists.txt
similarity index 100%
rename from examples/contact_mechanics/CMakeLists.txt
rename to examples/c++/contact_mechanics_model/CMakeLists.txt
diff --git a/examples/c++/contact_mechanics_model/README.rst b/examples/c++/contact_mechanics_model/README.rst
new file mode 100644
index 000000000..35200d7c3
--- /dev/null
+++ b/examples/c++/contact_mechanics_model/README.rst
@@ -0,0 +1,9 @@
+Examples for the contact mechanics model
+````````````````````````````````````````
+
+`contact_explicit_static` and `contact_explicit_dynamic` are solving a 2D Hertz contact patch test.
+
+.. image:: images/hertz.svg
+
+
+.. image:: images/hertz.png
diff --git a/examples/contact_mechanics/cohesive-contact.geo b/examples/c++/contact_mechanics_model/cohesive-contact.geo
similarity index 100%
rename from examples/contact_mechanics/cohesive-contact.geo
rename to examples/c++/contact_mechanics_model/cohesive-contact.geo
diff --git a/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc b/examples/c++/contact_mechanics_model/cohesive_contact_explicit_dynamic.cc
similarity index 100%
rename from examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc
rename to examples/c++/contact_mechanics_model/cohesive_contact_explicit_dynamic.cc
diff --git a/examples/contact_mechanics/contact_explicit_dynamic.cc b/examples/c++/contact_mechanics_model/contact_explicit_dynamic.cc
similarity index 100%
rename from examples/contact_mechanics/contact_explicit_dynamic.cc
rename to examples/c++/contact_mechanics_model/contact_explicit_dynamic.cc
diff --git a/examples/contact_mechanics/contact_explicit_static.cc b/examples/c++/contact_mechanics_model/contact_explicit_static.cc
similarity index 100%
rename from examples/contact_mechanics/contact_explicit_static.cc
rename to examples/c++/contact_mechanics_model/contact_explicit_static.cc
diff --git a/examples/contact_mechanics/hertz.geo b/examples/c++/contact_mechanics_model/hertz.geo
similarity index 100%
rename from examples/contact_mechanics/hertz.geo
rename to examples/c++/contact_mechanics_model/hertz.geo
diff --git a/examples/c++/contact_mechanics_model/hertz.png b/examples/c++/contact_mechanics_model/hertz.png
new file mode 100644
index 000000000..97d9963fe
Binary files /dev/null and b/examples/c++/contact_mechanics_model/hertz.png differ
diff --git a/examples/c++/contact_mechanics_model/hertz.svg b/examples/c++/contact_mechanics_model/hertz.svg
new file mode 100644
index 000000000..a4b76fc53
--- /dev/null
+++ b/examples/c++/contact_mechanics_model/hertz.svg
@@ -0,0 +1,921 @@
+
+
diff --git a/examples/contact_mechanics/material-cohesive.dat b/examples/c++/contact_mechanics_model/material-cohesive.dat
similarity index 100%
rename from examples/contact_mechanics/material-cohesive.dat
rename to examples/c++/contact_mechanics_model/material-cohesive.dat
diff --git a/examples/contact_mechanics/material.dat b/examples/c++/contact_mechanics_model/material.dat
similarity index 100%
rename from examples/contact_mechanics/material.dat
rename to examples/c++/contact_mechanics_model/material.dat
diff --git a/examples/heat_transfer/CMakeLists.txt b/examples/c++/diffusion_model/CMakeLists.txt
similarity index 100%
rename from examples/heat_transfer/CMakeLists.txt
rename to examples/c++/diffusion_model/CMakeLists.txt
diff --git a/examples/heat_transfer/cube.geo b/examples/c++/diffusion_model/cube.geo
similarity index 100%
rename from examples/heat_transfer/cube.geo
rename to examples/c++/diffusion_model/cube.geo
diff --git a/examples/heat_transfer/cube_tet.geo b/examples/c++/diffusion_model/cube_tet.geo
similarity index 100%
rename from examples/heat_transfer/cube_tet.geo
rename to examples/c++/diffusion_model/cube_tet.geo
diff --git a/examples/heat_transfer/heat_transfer_dynamics_2d.cc b/examples/c++/diffusion_model/heat_transfer_dynamics_2d.cc
similarity index 100%
rename from examples/heat_transfer/heat_transfer_dynamics_2d.cc
rename to examples/c++/diffusion_model/heat_transfer_dynamics_2d.cc
diff --git a/examples/heat_transfer/heat_transfer_dynamics_3d.cc b/examples/c++/diffusion_model/heat_transfer_dynamics_3d.cc
similarity index 100%
rename from examples/heat_transfer/heat_transfer_dynamics_3d.cc
rename to examples/c++/diffusion_model/heat_transfer_dynamics_3d.cc
diff --git a/examples/heat_transfer/heat_transfer_static_2d.cc b/examples/c++/diffusion_model/heat_transfer_static_2d.cc
similarity index 100%
rename from examples/heat_transfer/heat_transfer_static_2d.cc
rename to examples/c++/diffusion_model/heat_transfer_static_2d.cc
diff --git a/examples/heat_transfer/material.dat b/examples/c++/diffusion_model/material.dat
similarity index 100%
rename from examples/heat_transfer/material.dat
rename to examples/c++/diffusion_model/material.dat
diff --git a/examples/heat_transfer/square.geo b/examples/c++/diffusion_model/square.geo
similarity index 100%
rename from examples/heat_transfer/square.geo
rename to examples/c++/diffusion_model/square.geo
diff --git a/examples/embedded/CMakeLists.txt b/examples/c++/embedded_model/CMakeLists.txt
similarity index 100%
rename from examples/embedded/CMakeLists.txt
rename to examples/c++/embedded_model/CMakeLists.txt
diff --git a/examples/embedded/concrete.geo b/examples/c++/embedded_model/concrete.geo
similarity index 100%
rename from examples/embedded/concrete.geo
rename to examples/c++/embedded_model/concrete.geo
diff --git a/examples/embedded/embedded.cc b/examples/c++/embedded_model/embedded.cc
similarity index 100%
rename from examples/embedded/embedded.cc
rename to examples/c++/embedded_model/embedded.cc
diff --git a/examples/embedded/material.dat b/examples/c++/embedded_model/material.dat
similarity index 100%
rename from examples/embedded/material.dat
rename to examples/c++/embedded_model/material.dat
diff --git a/examples/embedded/reinforcement.geo b/examples/c++/embedded_model/reinforcement.geo
similarity index 100%
rename from examples/embedded/reinforcement.geo
rename to examples/c++/embedded_model/reinforcement.geo
diff --git a/examples/phase_field/CMakeLists.txt b/examples/c++/phase_field_model/CMakeLists.txt
similarity index 100%
rename from examples/phase_field/CMakeLists.txt
rename to examples/c++/phase_field_model/CMakeLists.txt
diff --git a/examples/phase_field/material_notch.dat b/examples/c++/phase_field_model/material_notch.dat
similarity index 100%
rename from examples/phase_field/material_notch.dat
rename to examples/c++/phase_field_model/material_notch.dat
diff --git a/examples/phase_field/phase_field_notch.cc b/examples/c++/phase_field_model/phase_field_notch.cc
similarity index 100%
rename from examples/phase_field/phase_field_notch.cc
rename to examples/c++/phase_field_model/phase_field_notch.cc
diff --git a/examples/phase_field/phase_field_parallel.cc b/examples/c++/phase_field_model/phase_field_parallel.cc
similarity index 100%
rename from examples/phase_field/phase_field_parallel.cc
rename to examples/c++/phase_field_model/phase_field_parallel.cc
diff --git a/examples/phase_field/square_notch.geo b/examples/c++/phase_field_model/square_notch.geo
similarity index 100%
rename from examples/phase_field/square_notch.geo
rename to examples/c++/phase_field_model/square_notch.geo
diff --git a/examples/cohesive_element/CMakeLists.txt b/examples/c++/solid_mechanics_cohesive_model/CMakeLists.txt
similarity index 100%
rename from examples/cohesive_element/CMakeLists.txt
rename to examples/c++/solid_mechanics_cohesive_model/CMakeLists.txt
diff --git a/examples/c++/solid_mechanics_cohesive_model/README.rst b/examples/c++/solid_mechanics_cohesive_model/README.rst
new file mode 100644
index 000000000..2e94052b0
--- /dev/null
+++ b/examples/c++/solid_mechanics_cohesive_model/README.rst
@@ -0,0 +1,2 @@
+Solid Mechanics Cohesive Model
+``````````````````````````````
diff --git a/examples/cohesive_element/cohesive_extrinsic/CMakeLists.txt b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/CMakeLists.txt
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic/CMakeLists.txt
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/CMakeLists.txt
diff --git a/examples/cohesive_element/cohesive_extrinsic/cohesive_extrinsic.cc b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/cohesive_extrinsic.cc
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic/cohesive_extrinsic.cc
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/cohesive_extrinsic.cc
diff --git a/examples/cohesive_element/cohesive_extrinsic/material.dat b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/material.dat
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic/material.dat
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/material.dat
diff --git a/examples/cohesive_element/cohesive_extrinsic/triangle.geo b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/triangle.geo
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic/triangle.geo
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/triangle.geo
diff --git a/examples/cohesive_element/cohesive_extrinsic/triangle.msh b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/triangle.msh
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic/triangle.msh
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic/triangle.msh
diff --git a/examples/cohesive_element/cohesive_extrinsic_ig_tg/CMakeLists.txt b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/CMakeLists.txt
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic_ig_tg/CMakeLists.txt
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/CMakeLists.txt
diff --git a/examples/cohesive_element/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc
diff --git a/examples/cohesive_element/cohesive_extrinsic_ig_tg/material.dat b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/material.dat
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic_ig_tg/material.dat
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/material.dat
diff --git a/examples/cohesive_element/cohesive_extrinsic_ig_tg/triangle.geo b/examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/triangle.geo
similarity index 100%
rename from examples/cohesive_element/cohesive_extrinsic_ig_tg/triangle.geo
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_extrinsic_ig_tg/triangle.geo
diff --git a/examples/cohesive_element/cohesive_intrinsic/CMakeLists.txt b/examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/CMakeLists.txt
similarity index 100%
rename from examples/cohesive_element/cohesive_intrinsic/CMakeLists.txt
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/CMakeLists.txt
diff --git a/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc b/examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/cohesive_intrinsic.cc
similarity index 100%
rename from examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/cohesive_intrinsic.cc
diff --git a/examples/cohesive_element/cohesive_intrinsic/material.dat b/examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/material.dat
similarity index 100%
rename from examples/cohesive_element/cohesive_intrinsic/material.dat
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/material.dat
diff --git a/examples/cohesive_element/cohesive_intrinsic/triangle.geo b/examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/triangle.geo
similarity index 100%
rename from examples/cohesive_element/cohesive_intrinsic/triangle.geo
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/triangle.geo
diff --git a/examples/cohesive_element/cohesive_intrinsic/triangle.msh b/examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/triangle.msh
similarity index 100%
rename from examples/cohesive_element/cohesive_intrinsic/triangle.msh
rename to examples/c++/solid_mechanics_cohesive_model/cohesive_intrinsic/triangle.msh
diff --git a/examples/c++/solid_mechanics_model/CMakeLists.txt b/examples/c++/solid_mechanics_model/CMakeLists.txt
new file mode 100644
index 000000000..a44ebf2a1
--- /dev/null
+++ b/examples/c++/solid_mechanics_model/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_example(new_material "Example on how to add a new material in Akantu" PACKAGE core)
+add_example(boundary_conditions "Example on hoy to set boundary conditions" PACKAGE core)
+add_example(explicit "Example on how to run an explicit simulation" PACKAGE core)
+add_example(io "Example on how to perform Input/Output operations" PACKAGE core)
+add_example(implicit "Example on how to run an implicit simulation" PACKAGE implicit)
+add_example(static "Example on how to run a static simulation" PACKAGE implicit)
+add_example(parallel "Example of how to write a parallel code with Akantu" PACKAGE parallel)
diff --git a/examples/c++/solid_mechanics_model/README.rst b/examples/c++/solid_mechanics_model/README.rst
new file mode 100644
index 000000000..4eff9735f
--- /dev/null
+++ b/examples/c++/solid_mechanics_model/README.rst
@@ -0,0 +1,2 @@
+Solid Mechanics Model
+`````````````````````
diff --git a/examples/boundary_conditions/CMakeLists.txt b/examples/c++/solid_mechanics_model/boundary_conditions/CMakeLists.txt
similarity index 100%
rename from examples/boundary_conditions/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/boundary_conditions/CMakeLists.txt
diff --git a/examples/boundary_conditions/predefined_bc/CMakeLists.txt b/examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/CMakeLists.txt
similarity index 100%
rename from examples/boundary_conditions/predefined_bc/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/CMakeLists.txt
diff --git a/examples/boundary_conditions/predefined_bc/material.dat b/examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/material.dat
similarity index 100%
rename from examples/boundary_conditions/predefined_bc/material.dat
rename to examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/material.dat
diff --git a/examples/boundary_conditions/predefined_bc/predefined_bc.cc b/examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/predefined_bc.cc
similarity index 100%
rename from examples/boundary_conditions/predefined_bc/predefined_bc.cc
rename to examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/predefined_bc.cc
diff --git a/examples/boundary_conditions/predefined_bc/square.geo b/examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/square.geo
similarity index 100%
rename from examples/boundary_conditions/predefined_bc/square.geo
rename to examples/c++/solid_mechanics_model/boundary_conditions/predefined_bc/square.geo
diff --git a/examples/boundary_conditions/python_user_defined_bc/CMakeLists.txt b/examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/CMakeLists.txt
similarity index 100%
rename from examples/boundary_conditions/python_user_defined_bc/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/CMakeLists.txt
diff --git a/examples/boundary_conditions/python_user_defined_bc/boundary_condition.py b/examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/boundary_condition.py
similarity index 100%
rename from examples/boundary_conditions/python_user_defined_bc/boundary_condition.py
rename to examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/boundary_condition.py
diff --git a/examples/boundary_conditions/python_user_defined_bc/fine_mesh.geo b/examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/fine_mesh.geo
similarity index 100%
rename from examples/boundary_conditions/python_user_defined_bc/fine_mesh.geo
rename to examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/fine_mesh.geo
diff --git a/examples/boundary_conditions/python_user_defined_bc/material.dat b/examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/material.dat
similarity index 100%
rename from examples/boundary_conditions/python_user_defined_bc/material.dat
rename to examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/material.dat
diff --git a/examples/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc b/examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc
similarity index 100%
rename from examples/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc
rename to examples/c++/solid_mechanics_model/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc
diff --git a/examples/boundary_conditions/user_defined_bc/CMakeLists.txt b/examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/CMakeLists.txt
similarity index 100%
rename from examples/boundary_conditions/user_defined_bc/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/CMakeLists.txt
diff --git a/examples/boundary_conditions/user_defined_bc/fine_mesh.geo b/examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/fine_mesh.geo
similarity index 100%
rename from examples/boundary_conditions/user_defined_bc/fine_mesh.geo
rename to examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/fine_mesh.geo
diff --git a/examples/boundary_conditions/user_defined_bc/material.dat b/examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/material.dat
similarity index 100%
rename from examples/boundary_conditions/user_defined_bc/material.dat
rename to examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/material.dat
diff --git a/examples/boundary_conditions/user_defined_bc/user_defined_bc.cc b/examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/user_defined_bc.cc
similarity index 100%
rename from examples/boundary_conditions/user_defined_bc/user_defined_bc.cc
rename to examples/c++/solid_mechanics_model/boundary_conditions/user_defined_bc/user_defined_bc.cc
diff --git a/examples/explicit/CMakeLists.txt b/examples/c++/solid_mechanics_model/explicit/CMakeLists.txt
similarity index 100%
rename from examples/explicit/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/explicit/CMakeLists.txt
diff --git a/examples/c++/solid_mechanics_model/explicit/README.rst b/examples/c++/solid_mechanics_model/explicit/README.rst
new file mode 100644
index 000000000..6e6bf735d
--- /dev/null
+++ b/examples/c++/solid_mechanics_model/explicit/README.rst
@@ -0,0 +1,7 @@
+This examples simulates a pulse in a 3D bar, the pulse is of the form::
+
+.. math:: u_x = A \mathrm{sin}(\frac{2 \pi k}{L} X_x) e^{- \frac{\frac{2 \pi
+ k}{L} X_x}/{L}}^2
+
+with :math:`u_x` and :math:`X_x` being the displacement and position in the
+:math:`x` direction.
diff --git a/examples/explicit/bar.geo b/examples/c++/solid_mechanics_model/explicit/bar.geo
similarity index 100%
rename from examples/explicit/bar.geo
rename to examples/c++/solid_mechanics_model/explicit/bar.geo
diff --git a/examples/explicit/explicit_dynamic.cc b/examples/c++/solid_mechanics_model/explicit/explicit_dynamic.cc
similarity index 99%
rename from examples/explicit/explicit_dynamic.cc
rename to examples/c++/solid_mechanics_model/explicit/explicit_dynamic.cc
index cbc98f7df..9b6a44f1a 100644
--- a/examples/explicit/explicit_dynamic.cc
+++ b/examples/c++/solid_mechanics_model/explicit/explicit_dynamic.cc
@@ -1,98 +1,97 @@
/**
* Copyright (©) 2014-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 3;
const Real pulse_width = 2.;
const Real A = 0.01;
Real time_step;
Real time_factor = 0.8;
Int max_steps = 1000;
-
Mesh mesh(spatial_dimension);
if (Communicator::getStaticCommunicator().whoAmI() == 0)
mesh.read("bar.msh");
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass);
time_step = model.getStableTimeStep();
std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step
<< "s)" << std::endl;
time_step *= time_factor;
model.setTimeStep(time_step);
/// boundary and initial conditions
Array & displacement = model.getDisplacement();
const Array & nodes = mesh.getNodes();
for (Int n = 0; n < mesh.getNbNodes(); ++n) {
Real x = nodes(n) - 2;
// Sinus * Gaussian
Real L = pulse_width;
Real k = 0.1 * 2 * M_PI * 3 / L;
displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L));
}
std::ofstream energy;
energy.open("energy.csv");
energy << "id,rtime,epot,ekin,tot" << std::endl;
model.setBaseName("explicit_dynamic");
model.addDumpField("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("stress");
model.dump();
for (Int s = 1; s <= max_steps; ++s) {
model.solveStep();
Real epot = model.getEnergy("potential");
Real ekin = model.getEnergy("kinetic");
energy << s << "," << s * time_step << "," << epot << "," << ekin << ","
<< epot + ekin << "," << std::endl;
if (s % 10 == 0)
std::cout << "passing step " << s << "/" << max_steps << std::endl;
model.dump();
}
energy.close();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/c++/solid_mechanics_model/explicit/images/sphx_glr_explicit.png b/examples/c++/solid_mechanics_model/explicit/images/sphx_glr_explicit.png
new file mode 100644
index 000000000..f3028eb06
Binary files /dev/null and b/examples/c++/solid_mechanics_model/explicit/images/sphx_glr_explicit.png differ
diff --git a/examples/explicit/material.dat b/examples/c++/solid_mechanics_model/explicit/material.dat
similarity index 100%
rename from examples/explicit/material.dat
rename to examples/c++/solid_mechanics_model/explicit/material.dat
diff --git a/examples/implicit/CMakeLists.txt b/examples/c++/solid_mechanics_model/implicit/CMakeLists.txt
similarity index 100%
rename from examples/implicit/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/implicit/CMakeLists.txt
diff --git a/examples/implicit/beam.geo b/examples/c++/solid_mechanics_model/implicit/beam.geo
similarity index 100%
rename from examples/implicit/beam.geo
rename to examples/c++/solid_mechanics_model/implicit/beam.geo
diff --git a/examples/implicit/implicit_dynamic.cc b/examples/c++/solid_mechanics_model/implicit/implicit_dynamic.cc
similarity index 100%
rename from examples/implicit/implicit_dynamic.cc
rename to examples/c++/solid_mechanics_model/implicit/implicit_dynamic.cc
diff --git a/examples/implicit/material_dynamic.dat b/examples/c++/solid_mechanics_model/implicit/material_dynamic.dat
similarity index 100%
rename from examples/implicit/material_dynamic.dat
rename to examples/c++/solid_mechanics_model/implicit/material_dynamic.dat
diff --git a/examples/io/CMakeLists.txt b/examples/c++/solid_mechanics_model/io/CMakeLists.txt
similarity index 100%
rename from examples/io/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/io/CMakeLists.txt
diff --git a/examples/io/dumper/CMakeLists.txt b/examples/c++/solid_mechanics_model/io/dumper/CMakeLists.txt
similarity index 100%
rename from examples/io/dumper/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/io/dumper/CMakeLists.txt
diff --git a/examples/io/dumper/dumpable_interface.cc b/examples/c++/solid_mechanics_model/io/dumper/dumpable_interface.cc
similarity index 100%
rename from examples/io/dumper/dumpable_interface.cc
rename to examples/c++/solid_mechanics_model/io/dumper/dumpable_interface.cc
diff --git a/examples/io/dumper/dumper_low_level.cc b/examples/c++/solid_mechanics_model/io/dumper/dumper_low_level.cc
similarity index 100%
rename from examples/io/dumper/dumper_low_level.cc
rename to examples/c++/solid_mechanics_model/io/dumper/dumper_low_level.cc
diff --git a/examples/io/dumper/locomotive_tools.cc b/examples/c++/solid_mechanics_model/io/dumper/locomotive_tools.cc
similarity index 100%
rename from examples/io/dumper/locomotive_tools.cc
rename to examples/c++/solid_mechanics_model/io/dumper/locomotive_tools.cc
diff --git a/examples/io/dumper/locomotive_tools.hh b/examples/c++/solid_mechanics_model/io/dumper/locomotive_tools.hh
similarity index 100%
rename from examples/io/dumper/locomotive_tools.hh
rename to examples/c++/solid_mechanics_model/io/dumper/locomotive_tools.hh
diff --git a/examples/io/dumper/swiss_train.geo b/examples/c++/solid_mechanics_model/io/dumper/swiss_train.geo
similarity index 100%
rename from examples/io/dumper/swiss_train.geo
rename to examples/c++/solid_mechanics_model/io/dumper/swiss_train.geo
diff --git a/examples/io/parser/CMakeLists.txt b/examples/c++/solid_mechanics_model/io/parser/CMakeLists.txt
similarity index 100%
rename from examples/io/parser/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/io/parser/CMakeLists.txt
diff --git a/examples/io/parser/example_parser.cc b/examples/c++/solid_mechanics_model/io/parser/example_parser.cc
similarity index 100%
rename from examples/io/parser/example_parser.cc
rename to examples/c++/solid_mechanics_model/io/parser/example_parser.cc
diff --git a/examples/io/parser/input_file.dat b/examples/c++/solid_mechanics_model/io/parser/input_file.dat
similarity index 100%
rename from examples/io/parser/input_file.dat
rename to examples/c++/solid_mechanics_model/io/parser/input_file.dat
diff --git a/examples/io/parser/swiss_cheese.geo b/examples/c++/solid_mechanics_model/io/parser/swiss_cheese.geo
similarity index 100%
rename from examples/io/parser/swiss_cheese.geo
rename to examples/c++/solid_mechanics_model/io/parser/swiss_cheese.geo
diff --git a/examples/new_material/CMakeLists.txt b/examples/c++/solid_mechanics_model/new_material/CMakeLists.txt
similarity index 100%
rename from examples/new_material/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/new_material/CMakeLists.txt
diff --git a/examples/new_material/barre_trou.geo b/examples/c++/solid_mechanics_model/new_material/barre_trou.geo
similarity index 100%
rename from examples/new_material/barre_trou.geo
rename to examples/c++/solid_mechanics_model/new_material/barre_trou.geo
diff --git a/examples/new_material/local_material_damage.cc b/examples/c++/solid_mechanics_model/new_material/local_material_damage.cc
similarity index 100%
rename from examples/new_material/local_material_damage.cc
rename to examples/c++/solid_mechanics_model/new_material/local_material_damage.cc
diff --git a/examples/new_material/local_material_damage.hh b/examples/c++/solid_mechanics_model/new_material/local_material_damage.hh
similarity index 100%
rename from examples/new_material/local_material_damage.hh
rename to examples/c++/solid_mechanics_model/new_material/local_material_damage.hh
diff --git a/examples/new_material/local_material_damage_inline_impl.hh b/examples/c++/solid_mechanics_model/new_material/local_material_damage_inline_impl.hh
similarity index 100%
rename from examples/new_material/local_material_damage_inline_impl.hh
rename to examples/c++/solid_mechanics_model/new_material/local_material_damage_inline_impl.hh
diff --git a/examples/new_material/material.dat b/examples/c++/solid_mechanics_model/new_material/material.dat
similarity index 100%
rename from examples/new_material/material.dat
rename to examples/c++/solid_mechanics_model/new_material/material.dat
diff --git a/examples/new_material/new_local_material.cc b/examples/c++/solid_mechanics_model/new_material/new_local_material.cc
similarity index 100%
rename from examples/new_material/new_local_material.cc
rename to examples/c++/solid_mechanics_model/new_material/new_local_material.cc
diff --git a/examples/new_material/viscoelastic_maxwell/CMakeLists.txt b/examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/CMakeLists.txt
similarity index 100%
rename from examples/new_material/viscoelastic_maxwell/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/CMakeLists.txt
diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat b/examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat
similarity index 100%
rename from examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat
rename to examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat
diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc b/examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc
similarity index 100%
rename from examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc
rename to examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc
diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo b/examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo
similarity index 100%
rename from examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo
rename to examples/c++/solid_mechanics_model/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo
diff --git a/examples/parallel/CMakeLists.txt b/examples/c++/solid_mechanics_model/parallel/CMakeLists.txt
similarity index 100%
rename from examples/parallel/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/parallel/CMakeLists.txt
diff --git a/examples/parallel/material.dat b/examples/c++/solid_mechanics_model/parallel/material.dat
similarity index 100%
rename from examples/parallel/material.dat
rename to examples/c++/solid_mechanics_model/parallel/material.dat
diff --git a/examples/parallel/parallel_2d.cc b/examples/c++/solid_mechanics_model/parallel/parallel_2d.cc
similarity index 100%
rename from examples/parallel/parallel_2d.cc
rename to examples/c++/solid_mechanics_model/parallel/parallel_2d.cc
diff --git a/examples/parallel/square_2d.geo b/examples/c++/solid_mechanics_model/parallel/square_2d.geo
similarity index 100%
rename from examples/parallel/square_2d.geo
rename to examples/c++/solid_mechanics_model/parallel/square_2d.geo
diff --git a/examples/static/CMakeLists.txt b/examples/c++/solid_mechanics_model/static/CMakeLists.txt
similarity index 100%
rename from examples/static/CMakeLists.txt
rename to examples/c++/solid_mechanics_model/static/CMakeLists.txt
diff --git a/examples/static/material.dat b/examples/c++/solid_mechanics_model/static/material.dat
similarity index 100%
rename from examples/static/material.dat
rename to examples/c++/solid_mechanics_model/static/material.dat
diff --git a/examples/static/square.geo b/examples/c++/solid_mechanics_model/static/square.geo
similarity index 100%
rename from examples/static/square.geo
rename to examples/c++/solid_mechanics_model/static/square.geo
diff --git a/examples/static/static.cc b/examples/c++/solid_mechanics_model/static/static.cc
similarity index 100%
rename from examples/static/static.cc
rename to examples/c++/solid_mechanics_model/static/static.cc
diff --git a/examples/structural_mechanics/CMakeLists.txt b/examples/c++/structural_mechanics_model/CMakeLists.txt
similarity index 100%
rename from examples/structural_mechanics/CMakeLists.txt
rename to examples/c++/structural_mechanics_model/CMakeLists.txt
diff --git a/examples/structural_mechanics/bernoulli_beam_2_example.cc b/examples/c++/structural_mechanics_model/bernoulli_beam_2_example.cc
similarity index 100%
rename from examples/structural_mechanics/bernoulli_beam_2_example.cc
rename to examples/c++/structural_mechanics_model/bernoulli_beam_2_example.cc
diff --git a/examples/python/CMakeLists.txt b/examples/python/CMakeLists.txt
index a8e77e062..98e5e8a58 100644
--- a/examples/python/CMakeLists.txt
+++ b/examples/python/CMakeLists.txt
@@ -1,41 +1,35 @@
#===============================================================================
# Copyright (©) 2016-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This file is part of Akantu
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with Akantu. If not, see .
#
#===============================================================================
-add_subdirectory(custom-material)
-add_subdirectory(dynamics)
-add_subdirectory(eigen_modes)
-add_subdirectory(plate-hole)
-add_subdirectory(stiffness_matrix)
-
-add_example(structural_mechanics "structural mechanics example in python"
+add_example(structural_mechanics_model "structural mechanics example in python"
PACKAGE structural_mechanics)
-add_example(fragmentation "example of fragmentation in python"
- PACKAGE cohesive_element)
-add_example(phase-field "phase-field example in python"
+add_example(phase_field_model "phase-field example in python"
PACKAGE phase_field)
-add_example(cohesive "cohesive element examples in python"
+add_example(solid_mechanics_cohesive_model "cohesive element examples in python"
PACKAGE cohesive_element)
-add_example(contact-mechanics "contact mechanics example in python"
+add_example(solid_mechanics_model "solid_mechanics_model in python"
+ PACKAGE solid_mechanics)
+add_example(contact_mechanics_model "contact mechanics example in python"
PACKAGE contact_mechanics)
package_add_files_to_package(
examples/python/README.rst
)
diff --git a/examples/python/README.rst b/examples/python/README.rst
deleted file mode 100644
index 57361e88d..000000000
--- a/examples/python/README.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-In order to use the Akantu examples using the python mode, one has to source
-the file ``akantu_environement.sh``:
-
-| > source /akantu_environement.sh
-
-or for an installed akantu version:
-
-| > source /share/akantu/akantu_environement.sh
diff --git a/examples/python/contact-mechanics/CMakeLists.txt b/examples/python/contact_mechanics_model/CMakeLists.txt
similarity index 100%
rename from examples/python/contact-mechanics/CMakeLists.txt
rename to examples/python/contact_mechanics_model/CMakeLists.txt
diff --git a/examples/python/contact-mechanics/compression.dat b/examples/python/contact_mechanics_model/compression.dat
similarity index 100%
rename from examples/python/contact-mechanics/compression.dat
rename to examples/python/contact_mechanics_model/compression.dat
diff --git a/examples/python/contact-mechanics/compression.geo b/examples/python/contact_mechanics_model/compression.geo
similarity index 100%
rename from examples/python/contact-mechanics/compression.geo
rename to examples/python/contact_mechanics_model/compression.geo
diff --git a/examples/python/contact-mechanics/compression.py b/examples/python/contact_mechanics_model/compression.py
similarity index 100%
rename from examples/python/contact-mechanics/compression.py
rename to examples/python/contact_mechanics_model/compression.py
diff --git a/examples/python/contact-mechanics/detection-explicit.dat b/examples/python/contact_mechanics_model/detection-explicit.dat
similarity index 100%
rename from examples/python/contact-mechanics/detection-explicit.dat
rename to examples/python/contact_mechanics_model/detection-explicit.dat
diff --git a/examples/python/contact-mechanics/detection-explicit.msh b/examples/python/contact_mechanics_model/detection-explicit.msh
similarity index 100%
rename from examples/python/contact-mechanics/detection-explicit.msh
rename to examples/python/contact_mechanics_model/detection-explicit.msh
diff --git a/examples/python/contact-mechanics/detection-explicit.py b/examples/python/contact_mechanics_model/detection-explicit.py
similarity index 100%
rename from examples/python/contact-mechanics/detection-explicit.py
rename to examples/python/contact_mechanics_model/detection-explicit.py
diff --git a/examples/python/phase-field/CMakeLists.txt b/examples/python/phase_field_model/CMakeLists.txt
similarity index 100%
rename from examples/python/phase-field/CMakeLists.txt
rename to examples/python/phase_field_model/CMakeLists.txt
diff --git a/examples/python/phase-field/material.dat b/examples/python/phase_field_model/material.dat
similarity index 100%
rename from examples/python/phase-field/material.dat
rename to examples/python/phase_field_model/material.dat
diff --git a/examples/python/phase-field/material_static.dat b/examples/python/phase_field_model/material_static.dat
similarity index 100%
rename from examples/python/phase-field/material_static.dat
rename to examples/python/phase_field_model/material_static.dat
diff --git a/examples/python/phase-field/phasefield-dynamic.py b/examples/python/phase_field_model/phasefield-dynamic.py
similarity index 100%
rename from examples/python/phase-field/phasefield-dynamic.py
rename to examples/python/phase_field_model/phasefield-dynamic.py
diff --git a/examples/python/phase-field/phasefield-static.py b/examples/python/phase_field_model/phasefield-static.py
similarity index 100%
rename from examples/python/phase-field/phasefield-static.py
rename to examples/python/phase_field_model/phasefield-static.py
diff --git a/examples/python/phase-field/plate.geo b/examples/python/phase_field_model/plate.geo
similarity index 100%
rename from examples/python/phase-field/plate.geo
rename to examples/python/phase_field_model/plate.geo
diff --git a/examples/python/phase-field/plate_static.geo b/examples/python/phase_field_model/plate_static.geo
similarity index 100%
rename from examples/python/phase-field/plate_static.geo
rename to examples/python/phase_field_model/plate_static.geo
diff --git a/examples/python/solid_mechanics_cohesive_model/CMakeLists.txt b/examples/python/solid_mechanics_cohesive_model/CMakeLists.txt
new file mode 100644
index 000000000..9ed940edd
--- /dev/null
+++ b/examples/python/solid_mechanics_cohesive_model/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(fragmentation)
+add_subdirectory(cohesive)
diff --git a/examples/python/cohesive/CMakeLists.txt b/examples/python/solid_mechanics_cohesive_model/cohesive/CMakeLists.txt
similarity index 100%
rename from examples/python/cohesive/CMakeLists.txt
rename to examples/python/solid_mechanics_cohesive_model/cohesive/CMakeLists.txt
diff --git a/examples/python/cohesive/custom_material.py b/examples/python/solid_mechanics_cohesive_model/cohesive/custom_material.py
similarity index 100%
rename from examples/python/cohesive/custom_material.py
rename to examples/python/solid_mechanics_cohesive_model/cohesive/custom_material.py
diff --git a/examples/python/cohesive/local_material.dat b/examples/python/solid_mechanics_cohesive_model/cohesive/local_material.dat
similarity index 100%
rename from examples/python/cohesive/local_material.dat
rename to examples/python/solid_mechanics_cohesive_model/cohesive/local_material.dat
diff --git a/examples/python/cohesive/material.dat b/examples/python/solid_mechanics_cohesive_model/cohesive/material.dat
similarity index 100%
rename from examples/python/cohesive/material.dat
rename to examples/python/solid_mechanics_cohesive_model/cohesive/material.dat
diff --git a/examples/python/cohesive/plate.geo b/examples/python/solid_mechanics_cohesive_model/cohesive/plate.geo
similarity index 100%
rename from examples/python/cohesive/plate.geo
rename to examples/python/solid_mechanics_cohesive_model/cohesive/plate.geo
diff --git a/examples/python/cohesive/plate.py b/examples/python/solid_mechanics_cohesive_model/cohesive/plate.py
similarity index 100%
rename from examples/python/cohesive/plate.py
rename to examples/python/solid_mechanics_cohesive_model/cohesive/plate.py
diff --git a/examples/python/fragmentation/CMakeLists.txt b/examples/python/solid_mechanics_cohesive_model/fragmentation/CMakeLists.txt
similarity index 100%
rename from examples/python/fragmentation/CMakeLists.txt
rename to examples/python/solid_mechanics_cohesive_model/fragmentation/CMakeLists.txt
diff --git a/examples/python/fragmentation/fragmentation.py b/examples/python/solid_mechanics_cohesive_model/fragmentation/fragmentation.py
similarity index 100%
rename from examples/python/fragmentation/fragmentation.py
rename to examples/python/solid_mechanics_cohesive_model/fragmentation/fragmentation.py
diff --git a/examples/python/fragmentation/fragmentation_mesh.geo b/examples/python/solid_mechanics_cohesive_model/fragmentation/fragmentation_mesh.geo
similarity index 100%
rename from examples/python/fragmentation/fragmentation_mesh.geo
rename to examples/python/solid_mechanics_cohesive_model/fragmentation/fragmentation_mesh.geo
diff --git a/examples/python/fragmentation/material.dat b/examples/python/solid_mechanics_cohesive_model/fragmentation/material.dat
similarity index 100%
rename from examples/python/fragmentation/material.dat
rename to examples/python/solid_mechanics_cohesive_model/fragmentation/material.dat
diff --git a/examples/python/solid_mechanics_model/CMakeLists.txt b/examples/python/solid_mechanics_model/CMakeLists.txt
new file mode 100644
index 000000000..1e0a74775
--- /dev/null
+++ b/examples/python/solid_mechanics_model/CMakeLists.txt
@@ -0,0 +1,5 @@
+add_subdirectory(custom-material)
+add_subdirectory(dynamics)
+add_subdirectory(eigen_modes)
+add_subdirectory(plate-hole)
+add_subdirectory(stiffness_matrix)
diff --git a/examples/python/custom-material/CMakeLists.txt b/examples/python/solid_mechanics_model/custom-material/CMakeLists.txt
similarity index 100%
rename from examples/python/custom-material/CMakeLists.txt
rename to examples/python/solid_mechanics_model/custom-material/CMakeLists.txt
diff --git a/examples/python/custom-material/bar.geo b/examples/python/solid_mechanics_model/custom-material/bar.geo
similarity index 100%
rename from examples/python/custom-material/bar.geo
rename to examples/python/solid_mechanics_model/custom-material/bar.geo
diff --git a/examples/python/custom-material/bi-material.py b/examples/python/solid_mechanics_model/custom-material/bi-material.py
similarity index 100%
rename from examples/python/custom-material/bi-material.py
rename to examples/python/solid_mechanics_model/custom-material/bi-material.py
diff --git a/examples/python/custom-material/custom-material.py b/examples/python/solid_mechanics_model/custom-material/custom-material.py
similarity index 100%
rename from examples/python/custom-material/custom-material.py
rename to examples/python/solid_mechanics_model/custom-material/custom-material.py
diff --git a/examples/python/custom-material/material.dat b/examples/python/solid_mechanics_model/custom-material/material.dat
similarity index 100%
rename from examples/python/custom-material/material.dat
rename to examples/python/solid_mechanics_model/custom-material/material.dat
diff --git a/examples/python/custom-material/square.geo b/examples/python/solid_mechanics_model/custom-material/square.geo
similarity index 100%
rename from examples/python/custom-material/square.geo
rename to examples/python/solid_mechanics_model/custom-material/square.geo
diff --git a/examples/python/dynamics/CMakeLists.txt b/examples/python/solid_mechanics_model/dynamics/CMakeLists.txt
similarity index 100%
rename from examples/python/dynamics/CMakeLists.txt
rename to examples/python/solid_mechanics_model/dynamics/CMakeLists.txt
diff --git a/examples/python/dynamics/bar.geo b/examples/python/solid_mechanics_model/dynamics/bar.geo
similarity index 100%
rename from examples/python/dynamics/bar.geo
rename to examples/python/solid_mechanics_model/dynamics/bar.geo
diff --git a/examples/python/dynamics/dynamics.py b/examples/python/solid_mechanics_model/dynamics/dynamics.py
similarity index 100%
rename from examples/python/dynamics/dynamics.py
rename to examples/python/solid_mechanics_model/dynamics/dynamics.py
diff --git a/examples/python/dynamics/material.dat b/examples/python/solid_mechanics_model/dynamics/material.dat
similarity index 100%
rename from examples/python/dynamics/material.dat
rename to examples/python/solid_mechanics_model/dynamics/material.dat
diff --git a/examples/python/eigen_modes/CMakeLists.txt b/examples/python/solid_mechanics_model/eigen_modes/CMakeLists.txt
similarity index 100%
rename from examples/python/eigen_modes/CMakeLists.txt
rename to examples/python/solid_mechanics_model/eigen_modes/CMakeLists.txt
diff --git a/examples/python/eigen_modes/eigen_modes.py b/examples/python/solid_mechanics_model/eigen_modes/eigen_modes.py
similarity index 100%
rename from examples/python/eigen_modes/eigen_modes.py
rename to examples/python/solid_mechanics_model/eigen_modes/eigen_modes.py
diff --git a/examples/python/eigen_modes/image_saver.py b/examples/python/solid_mechanics_model/eigen_modes/image_saver.py
similarity index 100%
rename from examples/python/eigen_modes/image_saver.py
rename to examples/python/solid_mechanics_model/eigen_modes/image_saver.py
diff --git a/examples/python/eigen_modes/material.dat b/examples/python/solid_mechanics_model/eigen_modes/material.dat
similarity index 100%
rename from examples/python/eigen_modes/material.dat
rename to examples/python/solid_mechanics_model/eigen_modes/material.dat
diff --git a/examples/python/plate-hole/CMakeLists.txt b/examples/python/solid_mechanics_model/plate-hole/CMakeLists.txt
similarity index 100%
rename from examples/python/plate-hole/CMakeLists.txt
rename to examples/python/solid_mechanics_model/plate-hole/CMakeLists.txt
diff --git a/examples/python/plate-hole/material.dat b/examples/python/solid_mechanics_model/plate-hole/material.dat
similarity index 100%
rename from examples/python/plate-hole/material.dat
rename to examples/python/solid_mechanics_model/plate-hole/material.dat
diff --git a/examples/python/plate-hole/plate.geo b/examples/python/solid_mechanics_model/plate-hole/plate.geo
similarity index 100%
rename from examples/python/plate-hole/plate.geo
rename to examples/python/solid_mechanics_model/plate-hole/plate.geo
diff --git a/examples/python/plate-hole/plate.py b/examples/python/solid_mechanics_model/plate-hole/plate.py
similarity index 100%
rename from examples/python/plate-hole/plate.py
rename to examples/python/solid_mechanics_model/plate-hole/plate.py
diff --git a/examples/python/solver_callback/CMakeLists.txt b/examples/python/solid_mechanics_model/solver_callback/CMakeLists.txt
similarity index 100%
rename from examples/python/solver_callback/CMakeLists.txt
rename to examples/python/solid_mechanics_model/solver_callback/CMakeLists.txt
diff --git a/examples/python/solver_callback/bar.geo b/examples/python/solid_mechanics_model/solver_callback/bar.geo
similarity index 100%
rename from examples/python/solver_callback/bar.geo
rename to examples/python/solid_mechanics_model/solver_callback/bar.geo
diff --git a/examples/python/solver_callback/material.dat b/examples/python/solid_mechanics_model/solver_callback/material.dat
similarity index 100%
rename from examples/python/solver_callback/material.dat
rename to examples/python/solid_mechanics_model/solver_callback/material.dat
diff --git a/examples/python/solver_callback/solver_callback.py b/examples/python/solid_mechanics_model/solver_callback/solver_callback.py
similarity index 100%
rename from examples/python/solver_callback/solver_callback.py
rename to examples/python/solid_mechanics_model/solver_callback/solver_callback.py
diff --git a/examples/python/stiffness_matrix/CMakeLists.txt b/examples/python/solid_mechanics_model/stiffness_matrix/CMakeLists.txt
similarity index 100%
rename from examples/python/stiffness_matrix/CMakeLists.txt
rename to examples/python/solid_mechanics_model/stiffness_matrix/CMakeLists.txt
diff --git a/examples/python/stiffness_matrix/material.dat b/examples/python/solid_mechanics_model/stiffness_matrix/material.dat
similarity index 100%
rename from examples/python/stiffness_matrix/material.dat
rename to examples/python/solid_mechanics_model/stiffness_matrix/material.dat
diff --git a/examples/python/stiffness_matrix/plate.geo b/examples/python/solid_mechanics_model/stiffness_matrix/plate.geo
similarity index 100%
rename from examples/python/stiffness_matrix/plate.geo
rename to examples/python/solid_mechanics_model/stiffness_matrix/plate.geo
diff --git a/examples/python/stiffness_matrix/plate.msh b/examples/python/solid_mechanics_model/stiffness_matrix/plate.msh
similarity index 100%
rename from examples/python/stiffness_matrix/plate.msh
rename to examples/python/solid_mechanics_model/stiffness_matrix/plate.msh
diff --git a/examples/python/stiffness_matrix/stiffness_matrix.py b/examples/python/solid_mechanics_model/stiffness_matrix/stiffness_matrix.py
similarity index 100%
rename from examples/python/stiffness_matrix/stiffness_matrix.py
rename to examples/python/solid_mechanics_model/stiffness_matrix/stiffness_matrix.py
diff --git a/examples/python/structural_mechanics/CMakeLists.txt b/examples/python/structural_mechanics_model/CMakeLists.txt
similarity index 100%
rename from examples/python/structural_mechanics/CMakeLists.txt
rename to examples/python/structural_mechanics_model/CMakeLists.txt
diff --git a/examples/python/structural_mechanics/structural_mechanics.py b/examples/python/structural_mechanics_model/structural_mechanics.py
similarity index 100%
rename from examples/python/structural_mechanics/structural_mechanics.py
rename to examples/python/structural_mechanics_model/structural_mechanics.py
diff --git a/examples/python/structural_mechanics/structural_mechanics_dynamics.py b/examples/python/structural_mechanics_model/structural_mechanics_dynamics.py
similarity index 100%
rename from examples/python/structural_mechanics/structural_mechanics_dynamics.py
rename to examples/python/structural_mechanics_model/structural_mechanics_dynamics.py
diff --git a/examples/python/structural_mechanics/structural_mechanics_softening.py b/examples/python/structural_mechanics_model/structural_mechanics_softening.py
similarity index 100%
rename from examples/python/structural_mechanics/structural_mechanics_softening.py
rename to examples/python/structural_mechanics_model/structural_mechanics_softening.py