diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 044a41d8c..0062cf10f 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -1,92 +1,93 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 # @date last modification: Fri Jan 30 2015 # # @brief Build the documentation # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # 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 . # #=============================================================================== if(AKANTU_DOCUMENTATION_DEVELOPER_MANUAL) add_subdirectory(dev-doc) endif() file(GLOB_RECURSE __all_files "*.*") set(_all_files) foreach(_file ${__all_files}) if("${_file}" MATCHES "${PROJECT_SOURCE_DIR}/doc/manual") file(RELATIVE_PATH __file ${PROJECT_SOURCE_DIR} ${_file}) list(APPEND _all_files ${__file}) endif() endforeach() set(AKANTU_MANUAL_FILES) set(_akantu_manual_not_found_files) foreach(_pkg ${${_project}_PACKAGE_SYSTEM_PACKAGES_ON}) string(TOUPPER "${_pkg}" _pkg) foreach(_f ${AKANTU_${_pkg}_MANUAL_FILES}) #check if some file are registered but not present list(FIND _all_files doc/manual/${_f} _ret) if(_ret EQUAL -1) list(APPEND _akantu_manual_not_found_files "${_pkg}: doc/manual/${_f} ") else() list(APPEND AKANTU_MANUAL_FILES doc/manual/${_f}) endif() endforeach() endforeach() if(_akantu_manual_not_found_files) message("") message("******************************************************") message("There are files registered in packages but not present") message("******************************************************") foreach(_file ${_akantu_manual_not_found_files}) message(${_file}) endforeach() message("******************************************************") message(FATAL_ERROR "abort") endif() ################################################################ #construct the list of files to exclude because not registered ################################################################ set(_akantu_doc_exclude_files) foreach(_file ${_all_files}) list(FIND AKANTU_MANUAL_FILES ${_file} _ret) if(_ret EQUAL -1) list(APPEND _akantu_doc_exclude_files /${_file}) endif() endforeach() list(REMOVE_ITEM _akantu_doc_exclude_files /doc/manual/CMakeLists.txt) -set(AKANTU_DOC_EXCLUDE_FILES ${_akantu_doc_exclude_files} CACHE INTERNAL "Documentation files to excluse from Akantu Package" FORCE) +set(AKANTU_DOC_EXCLUDE_FILES ${_akantu_doc_exclude_files} + CACHE INTERNAL "Documentation files to excluse from Akantu Package" FORCE) if (AKANTU_DOCUMENTATION_MANUAL) add_subdirectory(manual) endif() diff --git a/doc/dev-doc/CMakeLists.txt b/doc/dev-doc/CMakeLists.txt index 35f27a4ea..18b4160ba 100644 --- a/doc/dev-doc/CMakeLists.txt +++ b/doc/dev-doc/CMakeLists.txt @@ -1,196 +1,199 @@ set(DOXYGEN_INPUT_DOX ${CMAKE_CURRENT_BINARY_DIR}/akantu.dox) set(DOXYGEN_XML_DIR ${CMAKE_CURRENT_BINARY_DIR}/xml) set(DOXYGEN_OUTPUT ${DOXYGEN_XML_DIR}/index.xml) # 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") # ---------------------------------------------------------------------------- # # Doxygen # # ---------------------------------------------------------------------------- # find_package(Doxygen REQUIRED) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) package_get_all_source_files( AKANTU_LIBRARY_SRCS AKANTU_LIBRARY_PUBLIC_HDRS AKANTU_LIBRARY_PRIVATE_HDRS ) package_get_all_include_directories( _akantu_include_dirs ) package_get_all_external_informations( PRIVATE_INCLUDE AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR INTERFACE_INCLUDE AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR LIBRARIES AKANTU_EXTERNAL_LIBRARIES ) list(APPEND _akantu_include_dirs ${AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR} ${AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR} ${PROJECT_BINARY_DIR}/src) file(STRINGS ${PROJECT_SOURCE_DIR}/.clang-format AKANTU_TAB_SIZE REGEX "^TabWidth: *([0-9]*)" ) string(REGEX REPLACE ".*([0-9]+)" "\\1" AKANTU_TAB_SIZE "${AKANTU_TAB_SIZE}") if (CMAKE_VERSION VERSION_GREATER 3.9.5) #set(DOXYGEN_WARNINGS YES) #set(DOXYGEN_QUIET NO) set(DOXYGEN_STRIP_FROM_PATH ${PROJECT_SOURCE_DIR}) set(DOXYGEN_STRIP_FROM_INC_PATH ${PROJECT_SOURCE_DIR}) set(DOXYGEN_TAB_SIZE ${AKANTU_TAB_SIZE}) set(DOXYGEN_ALIASES "rst=\\verbatim embed:rst" "endrst=\\endverbatim" ) set(DOXYGEN_WARN_IF_UNDOCUMENTED NO) set(DOXYGEN_WARN_IF_DOC_ERROR NO) set(DOXYGEN_WARN_AS_ERROR NO) set(DOXYGEN_EXCLUDE "${PROJECT_SOURCE_DIR}/src/common/aka_fwd.hh") set(DOXYGEN_RECURSIVE YES) set(DOXYGEN_EXCLUDE "aka_named_argument.hh" ) set(DOXYGEN_EXAMPLE_PATH "${PROJECT_SOURCE_DIR}/examples") set(DOXYGEN_EXAMPLE_RECURSIVE YES) set(DOXYGEN_SOURCE_BROWSER NO) set(DOXYGEN_CLANG_ASSISTED_PARSING NO) #set(DOXYGEN_CLANG_OPTIONS ) set(DOXYGEN_CLANG_DATABASE_PATH ${CMAKE_BINARY_DIR}) set(DOXYGEN_USE_MATHJAX YES) set(DOXYGEN_GENERATE_HTML NO) set(DOXYGEN_GENERATE_HTMLHELP NO) set(DOXYGEN_GENERATE_LATEX NO) set(DOXYGEN_GENERATE_XML YES) set(DOXYGEN_XML_OUTPUT xml) set(DOXYGEN_ENABLE_PREPROCESSING YES) set(DOXYGEN_MACRO_EXPANSION YES) set(DOXYGEN_INCLUDE_PATH ${_akantu_include_dirs}) -# set(DOXYGEN_PREDEFINED -# ${AKANTU_DEFINITIONS} -# "DOXYGEN" + set(DOXYGEN_PREDEFINED + ${AKANTU_DEFINITIONS} + "DOXYGEN" # "AKANTU_TO_IMPLEMENT()=" # "DECLARE_NAMED_ARGUMENT()=" # "OPTIONAL_NAMED_ARGUMENT(n, v)=v" # "REQUIRED_NAMED_ARGUMENT(n)=" -# ) + ) set(DOXYGEN_COLLABORATION_GRAPH NO) set(DOXYGEN_UML_LOOK YES) set(DOXYGEN_TEMPLATE_RELATIONS YES) set(DOXYGEN_CALL_GRAPH YES) set(DOXYGEN_CALLER_GRAPH YES) set(DOXYGEN_DOT_GRAPH_MAX_NODES 500) set(DOXYGEN_SHOW_FILES NO) set(DOXYGEN_LOOKUP_CACHE_SIZE 9) set(_SRCS + ${PROJECT_BINARY_DIR}/src/aka_config.hh + ${PROJECT_BINARY_DIR}/src/aka_element_classes_info.hh ${AKANTU_LIBRARY_SRCS} ${AKANTU_LIBRARY_PUBLIC_HDRS} ${AKANTU_LIBRARY_PRIVATE_HDRS} ) list(REMOVE_ITEM _SRCS "${PROJECT_SOURCE_DIR}/src/common/aka_named_argument.hh") doxygen_add_docs(doxygen-doc ${_SRCS} USE_STAMP_FILE COMMENT "Building XML documentation with Doxygen in ${DOXYGEN_XML_DIR}" ) else() string(REGEX REPLACE ";" " " AKANTU_DOXYGEN_DEFINTIONS "${AKANTU_DEFINITIONS};DOXYGEN") string(REGEX REPLACE ";" " " AKANTU_DOXYGEN_INCLUDE_DIRS "${_akantu_include_dirs}") make_directory(${DOXYGEN_XML_DIR}) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/akantu.dox.in ${DOXYGEN_INPUT_DOX} ) add_custom_command( OUTPUT ${DOXYGEN_OUTPUT} COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_INPUT_DOX} DEPENDS ${DOXYGEN_INPUT_DOX} COMMENT "Building XML documentation with Doxygen in ${DOXYGEN_XML_DIR}" ) add_custom_target(doxygen-doc ALL DEPENDS ${DOXYGEN_OUTPUT} ) add_custom_target(doxygen-doc-forced COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_INPUT} COMMENT "Building XML documentation with Doxygen (forced) in ${DOXYGEN_XML_DIR}" ) endif() # ---------------------------------------------------------------------------- # # Sphinx # # ---------------------------------------------------------------------------- # find_package(Sphinx REQUIRED) set(SPHINX_VERBOSE_FLAG "-q") if(CMAKE_VERBOSE_MAKEFILE) set(SPHINX_VERBOSE_FLAG) endif(CMAKE_VERBOSE_MAKEFILE) # set(AKANTU_IN_READTHEDOC TRUE) # configure_file( # "${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in" # "${CMAKE_CURRENT_SOURCE_DIR}/conf_rtd.py" # @ONLY) set(AKANTU_IN_READTHEDOC FALSE) configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in" - "${CMAKE_CURRENT_BINARY_DIR}/conf.py" + "${SPHINX_INPUT}" @ONLY) 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 -c "${CMAKE_CURRENT_BINARY_DIR}" -d "${SPHINX_CACHE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}" "${SPHINX_HTML_DIR}" ) file(GLOB_RECURSE _SPHINX_SRCS "*.rst") add_custom_command( OUTPUT ${SPHINX_OUTPUT} COMMAND ${_sphinx_command} DEPENDS doxygen-doc ${SPHINX_INPUT} ${_SPHINX_SRCS} COMMENT "Building HTML documentation with Sphinx in ${SPHINX_HTML_DIR}" ) add_custom_target(sphinx-doc ALL DEPENDS ${SPHINX_OUTPUT}) add_custom_target(sphinx-doc-forced COMMAND ${_sphinx_command} DEPENDS doxygen-doc COMMENT "Building HTML documentation with Sphinx (forced) in ${SPHINX_HTML_DIR}" ) diff --git a/doc/dev-doc/conf.py.in b/doc/dev-doc/conf.py.in index 8f53ac5b9..4cc768f29 100644 --- a/doc/dev-doc/conf.py.in +++ b/doc/dev-doc/conf.py.in @@ -1,225 +1,236 @@ # -*- coding: utf-8 -*- # # Configuration file for the Sphinx documentation builder. # # This file does only contain a selection of the most common options. For a # full list see the documentation: # http://www.sphinx-doc.org/en/master/config # -- Path setup -------------------------------------------------------------- # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # import os # import subprocess # import sys # sys.path.insert(0, os.path.abspath('.')) # -- Project information ----------------------------------------------------- project = '@PROJECT_NAME@' copyright = '@AKANTU_COPYRING@' author = '@AKANTU_MAINTAINER@' version = '@AKANTU_VERSION_MAJOR@.@AKANTU_VERSION_MINOR@' # The full version, including alpha/beta/rc tags release = '@AKANTU_VERSION@' # -- General configuration --------------------------------------------------- # If your documentation needs a minimal Sphinx version, state it here. # # needs_sphinx = '1.0' # Number figures numfig = True # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ 'sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.ifconfig', 'sphinx.ext.viewcode', + 'sphinxcontrib.bibtex', 'breathe', ] read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True' # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: # # source_suffix = ['.rst', '.md'] source_suffix = '.rst' # The master toctree document. master_doc = 'index' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. language = None # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. exclude_patterns = ["CMakeLists.txt"] # The name of the Pygments (syntax highlighting) style to use. pygments_style = 'sphinx' primary_domain = 'cpp' highlight_language = 'cpp' # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # if read_the_docs_build: html_theme = 'default' else: html_theme = 'sphinx_rtd_theme' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. # # html_theme_options = {} # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] # Custom sidebar templates, must be a dictionary that maps document names # to template names. # # The default sidebars (for documents that don't match any pattern) are # defined by theme itself. Builtin themes are using these templates by # default: ``['localtoc.html', 'relations.html', 'sourcelink.html', # 'searchbox.html']``. # # html_sidebars = {} html_sidebars = { '**': [ 'relations.html', # needs 'show_related': True theme option to display 'searchbox.html', ]} +math_eqref_format = "Eq. {number}" + # MathJax configuration mathjax_config = { + 'extensions': [ + "tex2jax.js", + "siunitx.js" + ], 'TeX': { 'Macros': { 'st': [r'\mathrm{#1}', 1], 'mat': [r'\mathbf{#1}', 1], }, + 'extensions': ["AMSmath.js", "AMSsymbols.js", "sinuitx.js"], }, } # -- Options for HTMLHelp output --------------------------------------------- # Output file base name for HTML help builder. htmlhelp_basename = 'Akantudoc' # -- Options for LaTeX output ------------------------------------------------ latex_elements = { # The paper size ('letterpaper' or 'a4paper'). # # 'papersize': 'letterpaper', # The font size ('10pt', '11pt' or '12pt'). # # 'pointsize': '10pt', # Additional stuff for the LaTeX preamble. # 'preamble': r'''\usepackage{amsmath}''', # Latex figure (float) alignment # # 'figure_align': 'htbp', } # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, # author, documentclass [howto, manual, or own class]). latex_documents = [ (master_doc, 'Akantu.tex', 'Akantu Documentation', 'Nicolas Richart', 'manual'), ] # -- Options for manual page output ------------------------------------------ # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). man_pages = [ (master_doc, 'akantu', 'Akantu Documentation', [author], 1) ] # -- Options for Texinfo output ---------------------------------------------- # Grouping the document tree into Texinfo files. List of tuples # (source start file, target name, title, author, # dir menu entry, description, category) texinfo_documents = [ (master_doc, 'Akantu', 'Akantu Documentation', author, 'Akantu', 'One line description of project.', 'Miscellaneous'), ] # -- Options for Epub output ------------------------------------------------- # Bibliographic Dublin Core info. epub_title = project epub_author = author epub_publisher = author epub_copyright = copyright # The unique identifier of the text. This can be a ISBN number # or the project homepage. # epub_identifier = '' # A unique identification for the text. # # epub_uid = '' # A list of files that should not be packed into the epub file. epub_exclude_files = ['search.html'] # -- Extension configuration ------------------------------------------------- if read_the_docs_build: akantu_path = "." else: akantu_path = "@CMAKE_CURRENT_BINARY_DIR@" os.makedirs("@CMAKE_CURRENT_BINARY_DIR@/_static", exist_ok=True) # print("akantu_path = '{}'".format(akantu_path)) breathe_projects = {"Akantu": os.path.join(akantu_path, "xml")} breathe_default_project = "Akantu" +breathe_default_members = ('members', 'undoc-members') +breathe_implementation_filename_extensions = ['.c', '.cc', '.cpp'] +breathe_show_enumvalue_initializer = True # -- Options for intersphinx extension --------------------------------------- intersphinx_mapping = { 'numpy': ('https://docs.scipy.org/doc/numpy/', None), 'scipy': ('https://docs.scipy.org/doc/scipy/reference', None), } diff --git a/doc/dev-doc/index.rst b/doc/dev-doc/index.rst index 2d0b0164d..8affc93b7 100644 --- a/doc/dev-doc/index.rst +++ b/doc/dev-doc/index.rst @@ -1,29 +1,34 @@ .. Akantu documentation master file, created by sphinx-quickstart on Fri Apr 17 16:35:46 2020. You can adapt this file completely to your liking, but it should at least contain the root `toctree` directive. Akantu: a FEM library ===================== .. toctree:: :maxdepth: 2 :caption: User Manual ./manual/getting_started.rst ./manual/fe_engine.rst ./manual/solidmechanicsmodel.rst .. toctree:: :maxdepth: 2 :caption: API Reference ./reference.rst +.. toctree:: + :maxdepth: 2 + :caption: Bibliography + + ./manual/bibliography.rst Indices and tables ================== * :ref:`genindex` * :ref:`modindex` * :ref:`search` diff --git a/doc/dev-doc/manual/bibliography.rst b/doc/dev-doc/manual/bibliography.rst new file mode 100644 index 000000000..7a9e4db7d --- /dev/null +++ b/doc/dev-doc/manual/bibliography.rst @@ -0,0 +1,7 @@ +.. _bibliography: + +Bibliography +------------ + +.. bibliography:: manual-bibliography.bib + :cited: diff --git a/doc/dev-doc/manual/constitutive-laws.rst b/doc/dev-doc/manual/constitutive-laws.rst new file mode 100644 index 000000000..019a908bd --- /dev/null +++ b/doc/dev-doc/manual/constitutive-laws.rst @@ -0,0 +1,534 @@ +.. _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 + 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 + `[sect:io:material] <#sect:io:material>`__). +| In order to conveniently store values at each quadrature in a material + point provides a special data structure, the . The internal fields are + inheriting from the . Furthermore, it provides several functions for + initialization, auto-resizing and auto removal of quadrature points. + +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: + +seed = 1.0 + +where the value 1.0 can be substituted with any number. Currently 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: + +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 . +In Appendix `[app:material-parameters] <#app:material-parameters>`__ a +summary of the parameters for all materials of is provided. + +Elasticity +`````````` + +The elastic law is a commonly used constitutive relationship that can be +used for a wide range of engineering materials (metals, concrete, rock, +wood, glass, rubber, etc.) provided that the strains remain small (small +deformation and stress lower than yield strength). + +The elastic laws are often expressed as :math:`\mat{\sigma} = +\mat{C}:\mat{\varepsilon}` with where :math:`\mat{\sigma}` is the Cauchy +stress tensor, :math:`\mat{\varepsilon}` represents the infinitesimal +strain tensor and :math:`\mat{C}` is the elastic modulus tensor. + +Linear isotropic +'''''''''''''''' + +.. figure:: figures/cl/stress_strain_el.svg + :align: center + :name: fig-smm-cl-elastic + + Stress-strain curve for elastic material and (b) schematic + representation of Hooke’slaw, denoted as a spring. + +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-elastic`. + +The equation that relates the strains to the displacements is: point) +from the displacements as follows: + +.. math:: \mat{\varepsilon} = + \frac{1}{2} \left[ \nabla_0 \vec{u}+\nabla_0 \vec{u}^T \right] + :label: eqn-smm-strain_inf + +where :math:`\mat{\varepsilon}` represents the infinitesimal strain +tensor, :math:`\nabla_{0}\vec{u}` the displacement gradient tensor +according to the initial configuration. The constitutive equation for +isotropic homogeneous media can be expressed as: + +.. math:: \mat{\sigma } = \lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu\mat{\varepsilon} + :label: eqn-smm-material-constitutive_elastic + +where :math:`\mat{\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} + +Linear anisotropic +'''''''''''''''''' + +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. For this we define the elastic modulus +tensor as follow: + +.. 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 tensor is expressed in the +material basis. And this basis as to be given too. This basis +:math:`\Omega_{\st{mat}} = \{\vec{n_1}, \vec{n_2}, \vec{n_3}\}` is used to define the rotation +:math:`R_{ij} = \vec{n_j} . \vec{e_i}`. And :math:`\mat{C}` can be rotated in the global +basis :math:`\Omega = \{\vec{e_1}, \vec{e_2}, \vec{e_3}\}` as follow: + +.. math:: + + \begin{aligned} + \mat{C}_{\Omega} &= \mat{R}_1 \mat{C}_{\Omega_{\st{mat}}} \mat{R}_2\\ + \mat{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]\\ + \mat{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} + + +Linear orthotropic +'''''''''''''''''' + +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 \\ + & \multicolumn{2}{l}{\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`. + +Neo-Hookean +``````````` + +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 + :align: center + :name: fig-smm-sl-neo-hookean + + Neo-hookean Stress-strain curve. + +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(\mat{C}) = \frac{1}{2}\lambda_0\left(\ln J\right)^2-\mu_0\ln J+\frac{1}{2} + \mu_0\left(\mathrm{tr}(\mat{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:`\mat{F}=\nabla_{\!\!\vec{X}}\vec{x}`): +:math:`J=\text{det}(\mat{F})`. Finally :math:`\mat{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:`\mat{\sigma}`) should be computed through the second +Piola-Kirchhoff stress tensor :math:`\mat{S}`: + +.. math:: \mat{\sigma } = \frac{1}{J}\mat{F}\mat{S}\mat{F}^T + +Finally the second Piola-Kirchhoff stress tensor is given by: + +.. math:: + + \mat{S} = 2\frac{\partial\Psi}{\partial\mat{C}} = \lambda_0\ln J + \mat{C}^{-1}+\mu_0\left(\mat{I}-\mat{C}^{-1}\right) + +The parameters to indicate in the material file are the same as those +for the elastic case: (Young’s modulus), (Poisson’s ratio). + +Visco-Elasticity +```````````````` + +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 +Figure `[fig:smm:cl:visco-elastic:hyst] <#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). + +The standard rheological linear solid model (see Sections 10.2 and 10.3 +of :cite:`simo92`) has been implemented in . This model +results from the combination of a spring mounted in parallel with a +spring and a dashpot connected in series, as illustrated in +Figure `[fig:smm:cl:visco-elastic:model] <#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 + E_V \right ) ^ {-1} \cdot \left [ \frac{d\sigma(t)}{dt} + \frac{E_V}{\eta}\sigma(t) - \frac{EE_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: +(density), (Young’s modulus), (Poisson’s ratio), (if set to zero plane +strain, otherwise plane stress), (dashpot viscosity) and (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. + +Small-Deformation Plasticity +```````````````````````````` + +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:: + + \begin{aligned} + \mat{\varepsilon} &= \mat{\varepsilon}^e +\mat{\varepsilon}^p\\ + {\mat{\sigma}} &= 2G(\mat{\varepsilon}^e) + \lambda \mathrm{tr}(\mat{\varepsilon}^e)\mat{I} + \end{aligned} + +.. figure:: figures/cl/isotropic_hardening_plasticity.svg + :name: fig:smm:cl:Lin-strain-hard + + 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:: f = \sigma_{\st{eff}} - \sigma_y = \left(\frac{3}{2} {\mat{\sigma}}^{\st{tr}} : {\mat{\sigma}}^{\st{tr}}\right)^\frac{1}{2}-\sigma_y (\mat{\varepsilon}^p) + :label: eqn-smm-constitutive-von_Mises + +.. math:: f < 0 \quad \textrm{Elastic deformation,} \qquad f = 0 \quad \textrm{Plastic deformation} + :label: eqn-smm-constitutive-yielding + +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:`{\mat{\sigma}}^{\st{tr}}` is the deviatoric part of stress given +by + +.. math:: + + \label{eqn:smm:constitutive:deviatoric stress} + {\mat{\sigma}}^{\st{tr}}=\mat{\sigma} - \frac{1}{3} \mathrm{tr}(\mat{\sigma}) \mat {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:: \Delta {\mat{\varepsilon}}^p = \Delta p \frac {\partial{f}}{\partial{\mat \sigma}}=\frac{3}{2} \Delta p \frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}} + :label: eqn-smm-constitutive-plastic_contitutive_equation + +In these expressions, the direction of the plastic strain increment (or +equivalently, plastic strain rate) is given by +:math:`\frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{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:: {\mat{\sigma}}^{\st{tr}} = {\mat{\sigma}}_t + 2G\Delta \mat{\varepsilon} + \lambda \mathrm{tr}(\Delta \mat{\varepsilon})\mat{I} + +#. Check the Yielding criteria: + + .. math:: f = (\frac{3}{2} {\mat{\sigma}}^{\st{tr}} : {\mat{\sigma}}^{\st{tr}})^{1/2}-\sigma_y (\mat{\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 {\mat{\varepsilon}}^p = \frac{3}{2} \Delta p \frac{{\mat{\sigma}}^{\st{tr}}}{\sigma_{\st{eff}}} + +#. Compute the stress increment: + + .. math:: {\Delta \mat{\sigma}} = 2G(\Delta \mat{\varepsilon}-\Delta \mat{\varepsilon}^p) + \lambda \mathrm{tr}(\Delta \mat{\varepsilon}-\Delta \mat{\varepsilon}^p)\mat{I} + +#. Update the variables: + + .. math:: \begin{aligned} + {\mat{\varepsilon^p}} &= {\mat{\varepsilon}}^p_t+{\Delta {\mat{\varepsilon}}^p}\\ + {\mat{\sigma}} &= {\mat{\sigma}}_t+{\Delta \mat{\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: (Yield stress) and (Hardening modulus). In addition, the +elastic parameters need to be defined as previously mentioned: (Young’s +modulus), (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:: \mat{\sigma} = (1-d)\, \mat{C}:\mat{\varepsilon} + +where :math:`\mat{\sigma}`, :math:`\mat{\varepsilon}` are the Cauchy +stress and strain tensors, and :math:`\mat{C}` is the elastic stiffness +tensor. This formulation relies on the definition of an evolution law +for the damage variable. In , many possibilities exist and they are +listed below. + +Marigo +'''''' + +This damage evolution law is energy based as defined by Marigo +:cite:`marigo81a,lemaitre96a`. It is an isotropic damage law. + +.. math:: + + \begin{aligned} + Y &= \frac{1}{2}\mat{\varepsilon}:\mat{C}:\mat{\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`. + +Mazars +'''''' + +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_{\st{eq}} = +\sqrt{<\mat{\varepsilon}>_+:<\mat{\varepsilon}>_+}`, with :math:`<.>_+` +the positive part of the tensor. This part is defined in the principal +coordinates (I, II, III) as :math:`\varepsilon_{\st{eq}} = +\sqrt{<\mat{\varepsilon_I}>_+^2 + <\mat{\varepsilon_{II}}>_+^2 + <\mat{\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_{\st{eq}}} - A_t \exp^{-B_t(\varepsilon_{\st{eq}}-\kappa_0)}\\ + D_c &= 1 - \frac{\kappa_0 (1- A_c)}{\varepsilon_{\st{eq}}} - A_c + \exp^{-B_c(\varepsilon_{\st{eq}}-\kappa_0)}\\ + \alpha_t &= \frac{\sum_{i=1}^3<\varepsilon_i>_+\varepsilon_{\st{nd}\;i}}{\varepsilon_{\st{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_{\st{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. diff --git a/doc/dev-doc/manual/fe_engine.rst b/doc/dev-doc/manual/fe_engine.rst index 5e693213f..b0a87ee43 100644 --- a/doc/dev-doc/manual/fe_engine.rst +++ b/doc/dev-doc/manual/fe_engine.rst @@ -1,183 +1,308 @@ ``FEEngine`` ============ -The ``FEEngine`` interface is dedicated to handle the -finite-element approximations and the numerical integration of the -weak form. As we will see in Chapter sect:smm, ``Model`` -creates its own ``FEEngine`` object so the explicit creation of the +The :cpp:class:`FEEngine` interface is dedicated to handle the +finite-element approximations and the numerical integration of the weak form. As +we will see in Chapter :doc:`./solidmechanicsmodel`, +:cpp:class:`Model` creates its own +:cpp:class:`FEEngine` object so the explicit creation of the object is not required. Mathematical Operations ----------------------- -Using the ``FEEngine`` object, one can compute a interpolation, an -integration or a gradient. A simple example is given below:: - - // having a FEEngine object - std::unique_ptr fem = - std::make_unique>( - my_mesh, dim, "my_fem"); - // instead of this, a FEEngine object can be get using the model: - // model.getFEEngine() - - //compute the gradient - Array u; //append the values you want - Array nablauq; //gradient array to be computed - // compute the gradient - fem->gradientOnIntegrationPoints(const Array &u, - Array &nablauq, - const UInt nb_degree_of_freedom, - ElementType type); - - // interpolate - Array uq; //interpolated array to be computed - // compute the interpolation - fem->interpolateOnIntegrationPoints(const Array &u, - Array &uq, - UInt nb_degree_of_freedom, - ElementType type); - - // interpolated function can be integrated over the elements - Array int_val_on_elem; - // integrate - fem->integrate(const Array &uq, - Array &int_uq, - UInt nb_degree_of_freedom, - ElementType type); - -Another example below shows how to integrate stress and strain fields -over elements assigned to a particular material:: - - UInt sp_dim = 3; //spatial dimension - UInt m = 1; //material index of interest - const ElementType type = _tetrahedron_4; //element type - - // get the stress and strain arrays associated to the material index m - const Array & strain_vec = model.getMaterial(m).getGradU(type); - const Array & stress_vec = model.getMaterial(m).getStress(type); - - // get the element filter for the material index - const Array & elem_filter = model.getMaterial(m).getElementFilter(type); - - // initialize the integrated stress and strain arrays - Array int_strain_vec(elem_filter.getSize(), - sp_dim*sp_dim, "int_of_strain"); - Array int_stress_vec(elem_filter.getSize(), - sp_dim*sp_dim, "int_of_stress"); - - // integrate the fields - model.getFEEngine().integrate(strain_vec, int_strain_vec, - sp_dim*sp_dim, type, _not_ghost, elem_filter); - model.getFEEngine().integrate(stress_vec, int_stress_vec, - sp_dim*sp_dim, type, _not_ghost, elem_filter); +Using the :cpp:class:`FEEngine` object, one can compute a interpolation, +an integration or a gradient.A simple example is given below: + +.. code-block:: c++ + + // having a FEEngine object + auto fem = std::make_unique>(my_mesh, dim, "my_fem"); + // instead of this, a FEEngine object can be get using the model: + // model.getFEEngine() + + // compute the gradient + Array u; // append the values you want + Array nablauq; // gradient array to be computed + // compute the gradient + fem->gradientOnIntegrationPoints(const Array & u, Array & nablauq, + const UInt nb_degree_of_freedom, + ElementType type); + + // interpolate + Array uq; // interpolated array to be computed + // compute the interpolation + fem->interpolateOnIntegrationPoints(const Array & u, Array & uq, + UInt nb_degree_of_freedom, + ElementType type); + + // interpolated function can be integrated over the elements + Array int_val_on_elem; + // integrate + fem->integrate(const Array & uq, Array & int_uq, + UInt nb_degree_of_freedom, ElementType type); + + +Another example below shows how to integrate stress and strain fields over + elements assigned to a particular material: + +.. code-block:: c++ + + UInt sp_dim{3}; // spatial dimension + UInt m{1}; // material index of interest + const auto type{_tetrahedron_4}; // element type + + // get the stress and strain arrays associated to the material index m + const auto & strain_vec = model.getMaterial(m).getGradU(type); + const auto & stress_vec = model.getMaterial(m).getStress(type); + + // get the element filter for the material index + const auto & elem_filter = model.getMaterial(m).getElementFilter(type); + + // initialize the integrated stress and strain arrays + Array int_strain_vec(elem_filter.getSize(), sp_dim * sp_dim, + "int_of_strain"); + Array int_stress_vec(elem_filter.getSize(), sp_dim * sp_dim, + "int_of_stress"); + + // integrate the fields + model.getFEEngine().integrate(strain_vec, int_strain_vec, sp_dim * sp_dim, type, + _not_ghost, elem_filter); + model.getFEEngine().integrate(stress_vec, int_stress_vec, sp_dim * sp_dim, type, + _not_ghost, elem_filter); + Elements -------- The base for every Finite-Elements computation is its mesh and the elements that are used within that mesh. The element types that can be used depend on the -mesh, but also on the dimensionality of the problem (1D, 2D or 3D). In Akantu, -several isoparametric Lagrangian element types are supported (and one -serendipity element). Each of these types is discussed in some detail below, +mesh, but also on the dimensionality of the problem (1D, 2D or 3D). In +``Akantu``, several iso-parametric Lagrangian element types are supported (and +one serendipity element). Each of these types is discussed in some detail below, starting with the 1D-elements all the way to the 3D-elements. More detailed information (shape function, location of Gaussian quadrature points, and so on) can be found in Appendix app:elements. -Isoparametric Elements -...................... +Iso-parametric Elements +....................... 1D ```` -In Akantu, there are two types of isoparametric elements defined in 1D. These -element types are called ``_segment_2`` and ``_segment_3``, and are -depicted schematically in :numref:`fig:elements:1D`. Some of the basic -properties of these elements are listed in :numref:`tab:elements:1D`. +There are two types of iso-parametric elements defined in 1D. These element +types are called :cpp:enumerator:`_segment_2 ` and +:cpp:enumerator:`_segment_3 `, and are depicted +schematically in :numref:`fig-elements-1D`. Some of the basic properties of +these elements are listed in :numref:`tab-elements-1D`. -.. _fig:elements:1D: +.. _fig-elements-1D: .. figure:: figures/elements/segments.svg :align: center - Schematic overview of the two 1D element types in Akantu. In each - element, the node numbering as used in Akantu is indicated and also the + Schematic overview of the two 1D element types in ``Akantu``. In each + element, the node numbering as used in ``Akantu`` is indicated and also the quadrature points are highlighted (gray circles). -.. _tab:elements:1D: -.. table:: Some basic properties of the two 1D isoparametric elements in Akantu +.. _tab-elements-1D: +.. csv-table:: Some basic properties of the two 1D iso-parametric elements in ``Akantu`` + :header: "Element type", "Order", "#nodes", "#quad points" - +--------------+---------+------+------+ - |Element |Order |#nodes|#quad | - |type | | |points| - +--------------+---------+------+------+ - |``_segment_2``|linear |2 |1 | - +--------------+---------+------+------+ - |``_segment_3``|quadratic|3 |2 | - +--------------+---------+------+------+ + ":cpp:enumerator:`_segment_2 `", "linear", 2, 1 + ":cpp:enumerator:`_segment_3 `", "quadratic", 3, 2 2D ```` -In Akantu, there are four types of isoparametric elements defined in 2D. These -element types are called ``_triangle_3``, ``_triangle_6``, -``_quadrangle_4`` and ``_quadrangle_8``, and all of them are depicted -in :numref:`fig:elements:2D`. As with the 1D elements, some of the most basic -properties of these elements are listed in :numref:`tab:elements:2D`. It is -important to note that the first element is linear, the next two quadratic and -the last one cubic. Furthermore, the last element type (``_quadrangle_8``) + +There are four types of iso-parametric elements defined in 2D. These element +types are called :cpp:enumerator:`_triangle_3 `, +:cpp:enumerator:`_triangle_6 `, +:cpp:enumerator:`_quadrangle_4 ` and +:cpp:enumerator:`_quadrangle_8 `, and all of them are +depicted in :numref:`fig-elements-2D`. As with the 1D elements, some of the most +basic properties of these elements are listed in :numref:`tab-elements-2D`. It +is important to note that the first element is linear, the next two quadratic +and the last one cubic. Furthermore, the last element type (``_quadrangle_8``) is not a Lagrangian but a serendipity element. -.. _fig:elements:2D: +.. _fig-elements-2D: .. figure:: figures/elements/elements_2d.svg :align: center - Schematic overview of the four 2D element types in Akantu. In each - element, the node numbering as used in Akantu is indicated and also the + Schematic overview of the four 2D element types in ``Akantu``. In each + element, the node numbering as used in ``Akantu`` is indicated and also the quadrature points are highlighted (gray circles). -.. _tab:elements:2D: -.. table:: Some basic properties of the 2D isoparametric elements in Akantu +.. _tab-elements-2D: +.. csv-table:: Some basic properties of the 2D iso-parametric elements in ``Akantu`` + :header: "Element type", "Order", "#nodes", "#quad points" - +--------------------+----------+------+------+ - |Element |Order |#nodes|#quad | - |type | | |points| - +--------------------+----------+------+------+ - |``_triangle_3`` |linear |3 |1 | - +--------------------+----------+------+------+ - |``_triangle_6`` |quadratic |6 |3 | - +--------------------+----------+------+------+ - |``_quadrangle_4`` |linear |4 |4 | - +--------------------+----------+------+------+ - |``_quadrangle_8`` |quadratic |8 |9 | - +--------------------+----------+------+------+ + ":cpp:enumerator:`_triangle_3 `", "linear", 3, 1 + ":cpp:enumerator:`_triangle_6 `", "quadratic", 6, 3 + ":cpp:enumerator:`_quadrangle_4 `", "linear", 4, 4 + ":cpp:enumerator:`_quadrangle_8 `", "quadratic", 8, 9 3D ```` -In Akantu, there are three types of isoparametric elements defined in 3D. These -element types are called ``_tetrahedron_4``, ``_tetrahedron_10`` and -``_hexahedron_8``, and all of them are depicted schematically in -:numref:`fig:elements:3D`. As with the 1D and 2D elements some of the most -basic properties of these elements are listed in :numref:`tab:elements:3D`. +In ``Akantu``, there are three types of iso-parametric elements defined in 3D. +These element types are called :cpp:enumerator:`_tetrahedron_4 +`, :cpp:enumerator:`_tetrahedron_10 +` and :cpp:enumerator:`_hexadedron_8 +`, and all of them are depicted schematically in +:numref:`fig-elements-3D`. As with the 1D and 2D elements some of the most basic +properties of these elements are listed in :numref:`tab-elements-3D`. -.. _fig:elements:3D: +.. _fig-elements-3D: .. figure:: figures/elements/elements_3d.svg :align: center - Schematic overview of the three 3D element types in Akantu. In each - element, the node numbering as used in Akantu is indicated and also the + Schematic overview of the three 3D element types in ``Akantu``. In each + element, the node numbering as used in ``Akantu`` is indicated and also the quadrature points are highlighted (gray circles). -.. _tab:elements:3D: -.. table:: Some basic properties of the 3D isoparametric elements in Akantu - - +--------------------+----------+------+------+ - |Element |Order |#nodes|#quad | - |type | | |points| - +--------------------+----------+------+------+ - |``_tetrahedron_4`` |linear |4 |1 | - +--------------------+----------+------+------+ - |``_tetrahedron_10`` |quadratic |10 |4 | - +--------------------+----------+------+------+ - |``_hexadedron_8`` |cubic |8 |8 | - +--------------------+----------+------+------+ +.. _tab-elements-3D: +.. csv-table:: Some basic properties of the 3D iso-parametric elements in ``Akantu`` + :header: "Element type", "Order", "#nodes", "#quad points" + + ":cpp:enumerator:`_tetrahedron_4 `", "linear", 4, 1 + ":cpp:enumerator:`_tetrahedron_10 `", "quadratic", 10, 4 + ":cpp:enumerator:`_hexadedron_8 `", "cubic", 8, 8 + +Cohesive Elements +................. + +The cohesive elements that have been implemented in ``Akantu`` are based +on the work of Ortiz and Pandolfi :cite:`ortiz1999`. Their main +properties are reported in :numref:`tab-coh-cohesive_elements`. + +.. _fig-smm-coh-cohesive2d: +.. figure:: figures/elements/cohesive_2d_6.svg + :align: center + + Cohesive element in 2D for quadratic triangular elements T6. + +.. _tab-coh-cohesive_elements: +.. csv-table:: Some basic properties of the cohesive elements in ``Akantu``. + :header: "Element type", "Facet type", "Order", "#nodes", "#quad points" + + ":cpp:enumerator:`_cohesive_1d_2 <_cohesive_1d_2>`", ":cpp:enumerator:`_point_1 `", "linear", 2, 1 + ":cpp:enumerator:`_cohesive_2d_4 `", ":cpp:enumerator:`_segment_2 `", "linear", 4, 1 + ":cpp:enumerator:`_cohesive_2d_6 `", ":cpp:enumerator:`_segment_3 `", "quadratic", 6, 2 + ":cpp:enumerator:`_cohesive_3d_6 `", ":cpp:enumerator:`_triangle_3 `","linear", 6, 1 + ":cpp:enumerator:`_cohesive_3d_12 `", ":cpp:enumerator:`_triangle_6 `", "quadratic", 12, 3 + +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_\mathrm{eff} > \sigma_\mathrm{c} \quad\text {with} \quad \sigma_\mathrm{eff} = \sqrt{\sigma_\mathrm{n} ^ 2 + \frac{\tau ^ 2} {\beta ^ 2 }} + +in which :math:`\sigma_\mathrm { n } +` is the tensile normal traction and $\tau$ the resulting tangential one( :numref:`fig-smm-coh-insertion`). + +For the static analysis of the structures containing cohesive elements, the +stiffness of the cohesive elements should also be added to the total stiffness +of the structure.Considering a 2D quadratic cohesive element as that in +:numref:`fig-smm-coh-cohesive2d`, the opening displacement along the mid-surface +can be written as: + +.. _eq-coh-opening: +.. math:: + \begin{align} + \vec{\Delta}(s) &= \left[\!\!\left[ \mat{u}\right]\!\!\right] \,\mat{N}(s)\\ + &= \begin{bmatrix} + u_3 - u_0 & u_4 - u_1 & u_5 - u_2\\ + v_3 - v_0 & v_4 - v_1 & v_5 - v_2 + \end{bmatrix} + \begin{bmatrix} + N_0(s)\\ + N_1(s)\\ + N_2(s) + \end{bmatrix}\\ + &= \mat{N}^\mathrm{k} \mat{A U} = \mat{PU} + \end{align} + +The :math:`\mat{U}`, :math:`\mat{A}` and :math:`\mat{N}^\mathrm{k}` are as following : + +.. math:: + \begin{align} + \mat{U} &= \left[\begin{array}{c c c c c c c c c c c c} + u_0 & v_0 & u_1 & v_1 & u_2 & v_2 & u_3 & v_3 & u_4 & v_4 & u_5 & v_5 + \end{array}\right]\\ + \mat{A} &= \left[\begin{array}{c c c c c c c c c c c c} + 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0\\ + 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0\\ + 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0\\ + 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0\\ + 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0\\ + 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 + \end{array}\right]\\ + \mat{N}^\mathrm{k} &= \begin{bmatrix} + N_0(s) & 0 & N_1(s) & 0 & N_2(s) & 0\ + 0 & N_0(s) & 0 & N_1(s) & 0 & N_2(s) + \end{bmatrix} + \end{align} + +The consistent stiffness matrix for the element is obtained as + +.. _eq-cohesive_stiffness: +.. math:: + \mat{K}=\int_{S_0}\mat{P}^\mathrm{T} \, \frac{\partial{\vec{T}}}{\partial{\delta}}\mat{P}\,\mathrm{d} S_0 + +where :math:`\vec{T}` is the cohesive traction and :math:`\delta` the opening +displacement (for more details check :numref:`tab-coh-cohesive_elements`). + + +Structural Elements +................... + +Bernoulli Beam Elements +``````````````````````` + +These elements allow to compute the displacements and rotations of +structures constituted by Bernoulli beams. ``Akantu`` defines them for +both 2D and 3D problems respectively in the element types +:cpp:enumerator:`_bernoulli_beam_2 ` and :cpp:enumerator:`_bernoulli_beam_3 `. A +schematic depiction of a beam element is shown in +:numref:`fig-elements-bernoulli` and some of its properties are +listed in :numref:`tab-elements-bernoulli`. + +.. note:: + Beam elements are of mixed order: the axial displacement is + linearly interpolated while transverse displacements and rotations + use cubic shape functions. + +.. _fig-elements-bernoulli: +.. figure:: figures/elements/bernoulli_2.svg + :align: center + + Schematic depiction of a Bernoulli beam element (applied to 2D and + 3D) in ``Akantu``. The node numbering as used in ``Akantu`` is + indicated, and also the quadrature points are highlighted (gray + circles). + +.. _tab-elements-bernoulli: +.. csv-table:: Some basic properties of the beam elements in ``Akantu`` + :header: "Element type", "Dimension", "# nodes", "# quad. points", "# d.o.f." + + ":cpp:enumerator:`_bernoulli_beam_2 `", "2D", 2, 3, 6 + ":cpp:enumerator:`_bernoulli_beam_3 `", "3D", 2, 3, 12 diff --git a/doc/dev-doc/manual/figures/cl/hooke_law.svg b/doc/dev-doc/manual/figures/cl/hooke_law.svg new file mode 100644 index 000000000..94e4b0062 --- /dev/null +++ b/doc/dev-doc/manual/figures/cl/hooke_law.svg @@ -0,0 +1,78 @@ + +image/svg+xmlE diff --git a/doc/dev-doc/manual/figures/cl/isotropic_hardening_plasticity.svg b/doc/dev-doc/manual/figures/cl/isotropic_hardening_plasticity.svg new file mode 100644 index 000000000..631e1c525 --- /dev/null +++ b/doc/dev-doc/manual/figures/cl/isotropic_hardening_plasticity.svg @@ -0,0 +1,1186 @@ + +image/svg+xmlσεEhσyLinear Isotropic HardeningPerfect PlasticityEεpεeσ diff --git a/doc/dev-doc/manual/figures/cl/stress_strain_el.svg b/doc/dev-doc/manual/figures/cl/stress_strain_el.svg new file mode 100644 index 000000000..e7973b480 --- /dev/null +++ b/doc/dev-doc/manual/figures/cl/stress_strain_el.svg @@ -0,0 +1,166 @@ + +image/svg+xmlεσEE diff --git a/doc/dev-doc/manual/figures/cl/stress_strain_hook.tex b/doc/dev-doc/manual/figures/cl/stress_strain_hook.tex new file mode 100644 index 000000000..71c671036 --- /dev/null +++ b/doc/dev-doc/manual/figures/cl/stress_strain_hook.tex @@ -0,0 +1,10 @@ +\documentclass[tikz,convert={outfile=\jobname.svg}]{standalone} +%\usetikzlibrary{...}% tikz package already loaded by 'tikz' option +\begin{document} +\begin{tikzpicture}% Example: + \draw[thick,latex-latex] (0,5) node[left] {$\sigma$} |- (5,0) node (x) [right, below] {$\varepsilon$}; + \draw[thin] (1.5,1.5) -- (2.5,1.5) -- (2.5,2.5) node [midway, right] {E}; + \draw[very thick,color=red] (0,0) -- (4,4); + \draw[very thick,latex-latex,color=red] (1,1) -- (3,3); +\end{tikzpicture} +\end{document} diff --git a/doc/dev-doc/manual/figures/cl/stress_strain_neo.svg b/doc/dev-doc/manual/figures/cl/stress_strain_neo.svg new file mode 100644 index 000000000..a183cf930 --- /dev/null +++ b/doc/dev-doc/manual/figures/cl/stress_strain_neo.svg @@ -0,0 +1,135 @@ + +image/svg+xml diff --git a/doc/dev-doc/manual/figures/cl/stress_strain_visco.svg b/doc/dev-doc/manual/figures/cl/stress_strain_visco.svg new file mode 100644 index 000000000..0de647977 --- /dev/null +++ b/doc/dev-doc/manual/figures/cl/stress_strain_visco.svg @@ -0,0 +1,616 @@ + +image/svg+xmlεσ diff --git a/doc/dev-doc/manual/figures/elements/bernoulli_2.svg b/doc/dev-doc/manual/figures/elements/bernoulli_2.svg new file mode 100644 index 000000000..9e96665f2 --- /dev/null +++ b/doc/dev-doc/manual/figures/elements/bernoulli_2.svg @@ -0,0 +1,878 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/dev-doc/manual/figures/elements/cohesive_2d_6.svg b/doc/dev-doc/manual/figures/elements/cohesive_2d_6.svg new file mode 100644 index 000000000..7bf0eca89 --- /dev/null +++ b/doc/dev-doc/manual/figures/elements/cohesive_2d_6.svg @@ -0,0 +1,192 @@ + +image/svg+xml132465 diff --git a/doc/dev-doc/manual/figures/insertion.svg b/doc/dev-doc/manual/figures/insertion.svg new file mode 100644 index 000000000..455503347 --- /dev/null +++ b/doc/dev-doc/manual/figures/insertion.svg @@ -0,0 +1,258 @@ + +image/svg+xmlboundaryloopdoubled nodes diff --git a/doc/dev-doc/manual/getting_started.rst b/doc/dev-doc/manual/getting_started.rst index dbbe5ad53..3b3957148 100644 --- a/doc/dev-doc/manual/getting_started.rst +++ b/doc/dev-doc/manual/getting_started.rst @@ -1,302 +1,306 @@ Getting Started =============== -Compiling Akantu ----------------- +Compiling ``Akantu`` +-------------------- -Akantu is a `CMake `_ project, so to configure it, you can either +``Akantu`` is a `CMake `_ project, so to configure it, you can either follow the usual way:: > cd akantu > mkdir build > cd build > ccmake .. [ Set the options that you need ] > make > make install -All the Akantu options are documented in Appendix app:package-dependencies. +All the ``Akantu`` options are documented in Appendix app:package-dependencies. Writing a ``main`` function --------------------------- -Akantu first needs to be initialized. The memory management included in the core -library handles the correct allocation and de-allocation of vectors, structures -and/or objects. Moreover, in parallel computations, the initialization procedure -performs the communication setup. This is achieved by the function -:cpp:func:`initialize ` that is used as follows:: +``Akantu`` first needs to be initialized. The memory management included in the +core library handles the correct allocation and de-allocation of vectors, +structures and/or objects. Moreover, in parallel computations, the +initialization procedure performs the communication setup. This is achieved by +the function :cpp:func:`initialize ` that is used as +follows:: #include "aka_common.hh" #include "..." using namespace akantu; int main(int argc, char *argv[]) { initialize("input_file.dat", argc, argv); // your code ... } The :cpp:func:`initialize ` function takes the text inpute -file and the program parameters which can be parsed by Akantu in due form (see +file and the program parameters which can be parsed by ``Akantu`` in due form (see sect:parser). Obviously it is necessary to include all files needed in main. In this manual all provided code implies the usage of ``akantu`` as namespace. Creating and Loading a Mesh --------------------------- -In its current state, Akantu supports three types of meshes: Gmsh, Abaqus and +In its current state, ``Akantu`` supports three types of meshes: Gmsh, Abaqus and Diana. Once a :cpp:class:`akantu::Mesh` object is created with a given spatial dimension, it can be filled by reading a mesh input file. The method -:cpp:member:`read ` of the class :cpp:class:`Mesh +:cpp:func:`read ` of the class :cpp:class:`Mesh ` infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified. :: UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); // Reading Gmsh files mesh.read("my_gmsh_mesh.msh"); mesh.read("my_gmsh_mesh", _miot_gmsh); The Gmsh reader adds the geometrical and physical tags as mesh data. The physical values are stored as a :cpp:type:`UInt ` data called ``tag_0``, if a string name is provided it is stored as a ``std::string`` data named ``physical_names``. The geometrical tag is stored as a :cpp:type:`UInt ` data named ``tag_1``. Using Arrays ------------ -Data in Akantu can be stored in data containers implemented by the +Data in ``Akantu`` can be stored in data containers implemented by the :cpp:class:`akantu::Array` class. In its most basic usage, the :cpp:class:`Array ` class implemented in \akantu is similar to the ``std::vector`` class of the Standard Template Library (STL) for C++. A simple :cpp:class:`Array ` containing a sequence of ``nb_element`` values (of a given type) can be generated with:: Array example_array(nb_element); -where ``type`` usually is ``Real``, ``Int``, ``UInt`` or ``bool``. -Each value is associated to an index, so that data can be accessed by -typing:: +where ``type`` usually is :cpp:type:`Real `, :cpp:type:`Int +`, :cpp:type:`UInt ` or ``bool``. Each value is +associated to an index, so that data can be accessed by typing:: auto & val = example_array(index); ``Arrays`` can also contain tuples of values for each index. In that case, the number of components per tuple must be specified at the :cpp:class:`Array ` creation. For example, if we want to create an :cpp:class:`Array ` to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following:: UInt nb_nodes = 10; UInt spatial_dimension = 3; Array position(nb_nodes, spatial_dimension); In this case the :math:`x` position of the eighth node number will be given by ``position(7, 0)`` (in C++, numbering starts at 0 and not 1). If the number of components for the sequences is not specified, the default value of 1 is used. Here is a list of some basic operations that can be performed on :cpp:class:`Array `: - - ``resize(size)`` change the size of the :cpp:class:`Array `. - - ``clear()`` set all entries of the :cpp:class:`Array ` to - zero. - - ``set(t)`` set all entries of the :cpp:class:`Array ` to - ``t``. - - ``copy(const Array & other)`` copy another :cpp:class:`Array - ` into the current one. The two :cpp:class:`Array - ` should have the same number of components. - - ``push_back(tuple)`` append a tuple with the correct number of components at - the end of the :cpp:class:`Array `. - - ``erase(i)`` erase the value at the i-th position. - - ``find(value)`` search ``value`` in the current :cpp:class:`Array - `. Return position index of the first occurence or -1 if not - found. - - ``storage()`` Return the address of the allocated memory of the - :cpp:class:`Array `. + - :cpp:func:`resize(size) ` change the size of + the :cpp:class:`Array `. + - :cpp:func:`clear ` reset the size of the + :cpp:class:`Array ` to zero. (*warning* this changed in > + v4.0) + - :cpp:func:`set(t) ` set all entries of the + :cpp:class:`Array ` to ``t``. + - :cpp:func:`copy(const Array & other) ` copy another + :cpp:class:`Array ` into the current one. The two + :cpp:class:`Arrays ` should have the same number of + components. + - :cpp:func:`push_back(tuple) ` append a tuple with + the correct number of components at the end of the :cpp:class:`Array `. + - :cpp:func:`erase(i) ` erase the value at the i-th position. + - :cpp:func:`find(value) ` search ``value`` in the + current :cpp:class:`Array `. Return position index of the + first occurence or -1 if not found. + - :cpp:func:`storage() ` Return the address of the + allocated memory of the :cpp:class:`Array `. Array iterators ------------------- -It is very common in Akantu to loop over arrays to perform a specific treatment. +It is very common in ``Akantu`` to loop over arrays to perform a specific treatment. This ranges from geometric calculation on nodal quantities to tensor algebra (in constitutive laws for example). The :cpp:class:`Array ` object has the possibility to request iterators in order to make the writing of loops easier and enhance readability. For instance, a loop over the nodal coordinates can be performed like:: // accessing the nodal coordinates Array // with spatial_dimension components const auto & nodes = mesh.getNodes(); for (const auto & coords : make_view(nodes, spatial_dimension)) { // do what you need .... } -In that example, each ``coords`` is a ``Vector`` containing -geometrical array of size ``spatial_dimension`` and the iteration is +In that example, each ``coords`` is a :cpp:class:`Vector\ ` +containing geometrical array of size ``spatial_dimension`` and the iteration is conveniently performed by the :cpp:class:`Array ` iterator. The :cpp:class:`Array ` object is intensively used to store second order tensor values. In that case, it should be specified that the returned object type is a matrix when constructing the iterator. This is done when calling the :cpp:func:`make_view `. For instance, assuming that we have a :cpp:class:`Array ` storing stresses, we can loop over the stored tensors by:: for (const auto & stress : make_view(stresses, spatial_dimension, spatial_dimension)) { // stress is of type `const Matrix&` } -In that last example, the :cpp:class:`Matrix ` objects are +In that last example, the :cpp:class:`Matrix\ ` objects are ``spatial_dimension`` :math:`\times` ``spatial_dimension`` matrices. The light -objects :cpp:class:`Matrix ` and :cpp:class:`Vector -` can be used and combined to do most common linear algebra. If -the number of component is 1, it is possible to use :cpp:func:`make_view -` to this effect. +objects :cpp:class:`Matrix\ ` and +:cpp:class:`Vector\ ` can be used and combined to do most +common linear algebra. If the number of component is 1, it is possible to use +:cpp:func:`make_view ` to this effect. In general, a mesh consists of several kinds of elements. Consequently, the amount of data to be stored can differ for each element type. The straightforward example is the connectivity array, namely the sequences of nodes belonging to each element (linear triangular elements have fewer nodes than, say, rectangular quadratic elements etc.). A particular data structure called -:cpp:class:`ElementTypeMapArray ` is provided to -easily manage this kind of data. It consists of a group of ``Arrays``, each +:cpp:class:`ElementTypeMapArray\ ` is provided +to easily manage this kind of data. It consists of a group of ``Arrays``, each associated to an element type. The following code can retrieve the -``ElementTypeMapArray`` which stores the connectivity arrays for a mesh:: +:cpp:class:`ElementTypeMapArray\ ` which +stores the connectivity arrays for a mesh:: const ElementTypeMapArray & connectivities = mesh.getConnectivities(); Then, the specific array associated to a given element type can be obtained by:: const Array & connectivity_triangle = connectivities(_triangle_3); where the first order 3-node triangular element was used in the presented piece of code. Vector & Matrix ``````````````` -The :cpp:class:`Array ` iterators as presented in the previous -section can be shaped as :cpp:class:`Vector ` or -:cpp:class:`Matrix `. This objects represent 1st and 2nd order +The :cpp:class:`Array\ ` iterators as presented in the previous +section can be shaped as :cpp:class:`Vector\ ` or +:cpp:class:`Matrix\ `. This objects represent 1st and 2nd order tensors. As such they come with some functionalities that we will present a bit more into detail in this here. ``Vector`` ''''''''''''' - Accessors: - ``v(i)`` gives the ``i`` -th component of the vector ``v`` - ``v[i]`` gives the ``i`` -th component of the vector ``v`` - ``v.size()`` gives the number of component - Level 1: (results are scalars) - ``v.norm()`` returns the geometrical norm (:math:`L_2`) - ``v.norm()`` returns the :math:`L_N` norm defined as :math:`\left(\sum_i |v(i)|^N\right)^{1/N}`. N can take any positive integer value. There are also some particular values for the most commonly used norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical norm and ``L_inf`` for the norm infinity. - ``v.dot(x)`` return the dot product of ``v`` and ``x`` - ``v.distance(x)`` return the geometrical norm of :math:`v - x` - Level 2: (results are vectors) - ``v += s``, ``v -= s``, ``v *= s``, ``v /= s`` those are element-wise operators that sum, substract, multiply or divide all the component of ``v`` by the scalar ``s`` - ``v += x``, ``v -= x`` sums or substracts the vector ``x`` to/from ``v`` - - ``v.mul(A, x, alpha)`` stores the result of :math:`\alpha \boldsymbol{A} - \vec{x}` in ``v``, :math:`\alpha` is equal to 1 by default + - ``v.mul(A, x, alpha)`` stores the result of :math:`\alpha \boldsymbol{A} \vec{x}` in ``v``, :math:`\alpha` is equal to 1 by default - ``v.solve(A, b)`` stores the result of the resolution of the system :math:`\boldsymbol{A} \vec{x} = \vec{b}` in ``v`` - ``v.crossProduct(v1, v2)`` computes the cross product of ``v1`` and ``v2`` and stores the result in ``v`` ``Matrix`` ''''''''''''' - Accessors: - ``A(i, j)`` gives the component :math:`A_{ij}` of the matrix ``A`` - ``A(i)`` gives the :math:`i^{th}` column of the matrix as a ``Vector`` - ``A[k]`` gives the :math:`k^{th}` component of the matrix, matrices are stored in a column major way, which means that to access :math:`A_{ij}`, :math:`k = i + j M` - ``A.rows()`` gives the number of rows of ``A`` (:math:`M`) - ``A.cols()`` gives the number of columns of ``A`` (:math:`N`) - ``A.size()`` gives the number of component in the matrix (:math:`M \times N`) - Level 1: (results are scalars) - ``A.norm()`` is equivalent to ``A.norm()`` - ``A.norm()`` returns the :math:`L_N` norm defined as :math:`\left(\sum_i\sum_j |A(i,j)|^N\right)^{1/N}`. N can take any positive integer value. There are also some particular values for the most commonly used norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical norm and ``L_inf`` for the norm infinity. - ``A.trace()`` return the trace of ``A`` - ``A.det()`` return the determinant of ``A`` - ``A.doubleDot(B)`` return the double dot product of ``A`` and ``B``, :math:`\mat{A}:\mat{B}` - Level 3: (results are matrices) - ``A.eye(s)``, ``Matrix::eye(s)`` fills/creates a matrix with the :math:`s\mat{I}` with :math:`\mat{I}` the identity matrix - ``A.inverse(B)`` stores :math:`\mat{B}^{-1}` in ``A`` - ``A.transpose()`` returns :math:`\mat{A}^{t}` - ``A.outerProduct(v1, v2)`` stores :math:`\vec{v_1} \vec{v_2}^{t}` in ``A`` - ``C.mul(A, B, alpha)``: stores the result of the product of ``A`` and code{B} time the scalar ``alpha`` in ``C``. ``t_A`` and ``t_B`` are boolean defining if ``A`` and ``B`` should be transposed or not. +----------+----------+--------------+ |``t_A`` |``t_B`` |result | | | | | +----------+----------+--------------+ |false |false |:math:`\mat{C}| | | |= \alpha | | | |\mat{A} | | | |\mat{B}` | | | | | +----------+----------+--------------+ |false |true |:math:`\mat{C}| | | |= \alpha | | | |\mat{A} | | | |\mat{B}^t` | | | | | +----------+----------+--------------+ |true |false |:math:`\mat{C}| | | |= \alpha | | | |\mat{A}^t | | | |\mat{B}` | | | | | +----------+----------+--------------+ |true |true |:math:`\mat{C}| | | |= \alpha | | | |\mat{A}^t | | | |\mat{B}^t` | +----------+----------+--------------+ - ``A.eigs(d, V)`` this method computes the eigenvalues and eigenvectors of ``A`` and store the results in ``d`` and ``V`` such that :math:`d(i) = \lambda_i` and :math:`V(i) = \vec{v_i}` with :math:`\mat{A}\vec{v_i} = \lambda_i\vec{v_i}` and :math:`\lambda_1 > ... > \lambda_i > ... > \lambda_N` diff --git a/doc/dev-doc/manual/manual-bibliography.bib b/doc/dev-doc/manual/manual-bibliography.bib new file mode 100644 index 000000000..4ebe2f8ff --- /dev/null +++ b/doc/dev-doc/manual/manual-bibliography.bib @@ -0,0 +1,582 @@ +% This file was created with JabRef 2.10b2. +% Encoding: UTF-8 + + +@Article{aifantis84a, + Title = {On the microstructural origin of certain inelastic models}, + Author = {E. C. Aifantis}, + Journal = {Journal of Engineering Materials and Technology}, + Year = {1984}, + Pages = {326 - 330}, + Volume = {106} +} + +@Article{Aragon:2013d, + Title = {A hierarchical detection framework for computational contact mechanics}, + Author = {Alejandro M. Arag{\'o}n and Jean-Fran{\c c}ois Molinari}, + Journal = {Computer Methods in Applied Mechanics and Engineering}, + Year = {2014}, + Number = {0}, + Pages = {574 - 588}, + Volume = {268}, + + Bdsk-url-1 = {http://www.sciencedirect.com/science/article/pii/S0045782513002533}, + Bdsk-url-2 = {http://dx.doi.org/10.1016/j.cma.2013.10.001}, + Date-added = {2014-05-22 09:14:48 +0000}, + Date-modified = {2014-05-22 09:14:48 +0000}, + Doi = {10.1016/j.cma.2013.10.001}, + ISSN = {0045-7825}, + Url = {http://dx.doi.org/10.1016/j.cma.2013.10.001} +} + +@Article{bazant86a, + Title = {Mechanics of distributed cracking}, + Author = {Zden\v{e}k P. Ba\v{z}ant}, + Journal = {Applied Mechanics Reviews}, + Year = {1986}, + Number = {5}, + Pages = {675 - 705}, + Volume = {39} +} + +@Article{bazant76a, + Title = {Instability, Ductility, and Size Effect in Strain-Softening Concrete}, + Author = {Zden\v{e}k P. Ba\v{z}ant}, + Journal = {Journal of the Engineering Mechanics Division}, + Year = {1976}, + Number = {5}, + Pages = {331 - 344}, + Volume = {102} +} + +@Article{bazant85a, + Title = {Wave Propagation in a Strain-Softening Bar: Exact Solution}, + Author = {Zden\v{e}k P. Ba\v{z}ant and Ted Belytschko}, + Journal = {Journal of Engineering Mechanics}, + Year = {1985}, + Number = {3}, + Pages = {381 - 389}, + Volume = {111} +} + +@Article{bazant88a, + Title = {{Nonlocal Continuum Damage, Localization Instability and Convergence}}, + Author = {Zden\v{e}k P. Ba\v{z}ant and Gilles Pijaudier-Cabot}, + Journal = {Journal of Applied Mechanics}, + Year = {1988}, + Number = {2}, + Pages = {287-293}, + Volume = {55}, + + Optdoi = {10.1115/1.3173674}, + Publisher = {ASME} +} + +@Article{belytschko03a, + Title = {Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment}, + Author = {T Belytschko and H Chen and J Xu and G Zi}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {2003}, + Pages = {1875-1905}, + Volume = {58} +} + +@Article{camacho_computational_1996, + Title = {Computational modelling of impact damage in brittle materials}, + Author = {Camacho, {G.T.} and Ortiz, M.}, + Journal = {International Journal of Solids and Structures}, + Year = {1996}, + + Month = aug, + Number = {20{\textendash}22}, + Pages = {2899--2938}, + Volume = {33}, + + Urldate = {2012-01-17} +} + +@Article{caughey1960, + Title = {Classical normal modes in damped linear dynamic systems}, + Author = {Caughey, TK}, + Journal = {Journal of Applied Mechanics}, + Year = {1960}, + Number = {2}, + Pages = {269--271}, + Volume = {27}, + + Publisher = {American Society of Mechanical Engineers} +} + +@Book{courant56a, + Title = {On the partial difference equations of mathematical physics}, + Author = {Courant, Richard and Friedrichs, Kurt Otto and Lewy, H.}, + Publisher = {Courant Institute of Mathematical Sciences, New York University}, + Year = {1956}, + + Address = {New York} +} + +@Book{curnier92a, + Title = {M{\'e}thodes num{\'e}riques en m{\'e}canique des solides}, + Author = {Curnier, A.}, + Publisher = {EPFL}, + Year = {1992}, + + Bdsk-url-1 = {http://books.google.ch/books?id=-Z2zygAACAAJ}, + Url = {http://books.google.ch/books?id=-Z2zygAACAAJ} +} + +@Article{youn10a, + Title = {Studies of dynamic crack propagation and crack branching with peridynamics}, + Author = {Youn Doh-Ha and Florin Bobaru}, + Journal = {International Journal of Fracture}, + Year = {2010}, + Pages = {229-1244}, + Volume = {162} +} + +@Book{frey2009, + Title = {M{\'e}thodes des {\'e}l{\'e}ments finis: Analyse des structures et milieux continus}, + Author = {Frey, F. and Jirousek, J.}, + Publisher = {PPUR}, + Year = {2009}, + Series = {Trait{\'e} de g{\'e}nie civil de l'Ecole Polytechnique F{\'e}d{\'e}rale de Lausanne}, + + Bdsk-url-1 = {http://books.google.fr/books?id=bCBtQgAACAAJ}, + ISBN = {9782880748524}, + Url = {http://books.google.fr/books?id=bCBtQgAACAAJ} +} + +@Article{gmsh, + Title = {Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities}, + Author = {Geuzaine, Christophe and Remacle, Jean-Fran{\c c}ois}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {2009}, + Number = {11}, + Pages = {1309--1331}, + Volume = {79}, + + Bdsk-url-1 = {http://dx.doi.org/10.1002/nme.2579}, + Doi = {10.1002/nme.2579}, + ISSN = {1097-0207}, + Keywords = {computer-aided design, mesh generation, post-processing, finite element method, open-source software}, + Publisher = {John Wiley \& Sons, Ltd.}, + Url = {http://dx.doi.org/10.1002/nme.2579} +} + +@Article{giry11a, + Title = {Stress-based nonlocal damage model}, + Author = {C. Giry and F. Dufour and J. Mazars}, + Journal = {International Journal of Solids and Structures}, + Year = {2011}, + Pages = {3431 - 3443}, + Volume = {48} +} + +@Article{hughes-83a, + Title = {A precis of developments in computational methods for transient analysis}, + Author = {Hughes, T.J.R. and Belytschko, T.}, + Journal = {Journal. of Applied Mechanics (ASME)}, + Year = {1983}, + Pages = {1033-1041}, + Volume = {50} +} + +@Article{hughes83a, + Title = {{A Pr{\'e}cis of Developments in Computational Methods for Transient Analysis}}, + Author = {{Hughes}, T.~J.~R. and {Belytschko}, T.}, + Journal = {Journal of Applied Mechanics}, + Year = {1983}, + + Month = {December}, + Pages = {1033}, + Volume = {50}, + + Optdoi = {10.1115/1.3167186} +} + +@Article{jirasek07a, + Title = {Mathematical analysis of strain localization}, + Author = {Milan Jir\'asek}, + Journal = {Revue Europeenne de Genie Civil}, + Year = {2007}, + Pages = {977 - 991}, + Volume = {11} +} + +@Article{jirasek07b, + Title = {Nonlocal damage mechanics}, + Author = {Milan Jir\'asek}, + Journal = {Revue Europeenne de Genie Civil}, + Year = {2007}, + Pages = {993 - 1021}, + Volume = {11} +} + +@Article{jirasek98a, + Title = {Nonlocal models for damage and fracture: comparison of approaches}, + Author = {Milan Jir\'asek}, + Journal = {International Journal of Solids and Structures}, + Year = {1998}, + Pages = {4133 - 4145}, + Volume = {35} +} + +@Article{jirasek03a, + Title = {Comparison of integral-type nonlocal plasticity models for strain-softening materials}, + Author = {Milan Jir\'asek and Simon Rolshoven}, + Journal = {International Journal of Engineering Science}, + Year = {2003}, + Pages = {1553-1602}, + Volume = {41} +} + +@Article{jirasek04a, + Title = {Size effect on fracture energy induced by non-locality}, + Author = {Milan Jir\'asek and S Rolshoven and P Grassl}, + Journal = {International Journal for numerical and analytical methods in geomechanics}, + Year = {2004}, + Pages = {653 - 670}, + Volume = {28} +} + +@Article{ladeveze92a, + Title = {A damage computational method for composite structures}, + Author = {P. Ladeveze}, + Journal = {Computational \& Structures}, + Year = {1992}, + Pages = {79-87}, + Volume = {44} +} + +@Book{Laursen:2002, + Title = {Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis}, + Author = {Laursen, T.A.}, + Publisher = {Springer}, + Year = {2002}, + Series = {Engineering online library}, + + Bdsk-url-1 = {http://books.google.ch/books?id=umzsErNuyFgC}, + Date-added = {2014-03-21 13:24:56 +0000}, + Date-modified = {2014-03-21 13:24:56 +0000}, + ISBN = {9783540429067}, + Lccn = {2002511027} +} + +@Book{lemaitre96a, + Title = {A Course on Damage Mechanics}, + Author = {Lemaitre, Jean}, + Publisher = {Springer Berlin Heidelberg}, + Year = {1996}, + + ISBN = {978-3-540-60980-3} +} + +@PhdThesis{levy10a, + Title = {Exploring the {P}hysics behind {D}ynamic {F}ragmentation through {P}arallel {S}imulations}, + Author = {Levy, Sarah}, + School = {ENAC}, + Year = {2010}, + + Address = {Lausanne}, + + Affiliation = {EPFL}, + Doctoral = {EDME}, + Institute = {IIC}, + Original-unit = {LSMS}, + Publisher = {EPFL}, + Unit = {LSMS} +} + +@Article{marigo81a, + Title = {{Formulation d'une loi d'endommagement d'un mat\'eriau \'elastique}}, + Author = {Marigo, Jean-Jacques}, + Journal = {{C. R. Acad. Sci., Paris, S\'er. II}}, + Year = {1981}, + Number = {19}, + Pages = {1309-1312}, + Volume = {292}, + + ISSN = {0249-6305} +} + +@PhdThesis{mazars84a, + Title = {{Application de la m\'ecanique de l'endommagement au comportement non lin\'eaire et \`a la rupture du b\'eton de structure}}, + Author = {Mazars, Jacky}, + School = {Universit\'e Paris 6}, + Year = {1984} +} + +@Article{needleman88a, + Title = {Material rate dependence and mesh sensitivity in localization problems}, + Author = {A. Needleman}, + Journal = {Computer Methods in Applied Mechanics and Engineering}, + Year = {1988}, + Number = {1}, + Pages = {69 - 85}, + Volume = {67}, + + ISSN = {0045-7825}, + Optdoi = {10.1016/0045-7825(88)90069-2} +} + +@Article{newmark59a, + Title = {{A Method of Computation for Structural Dynamics}}, + Author = {Nathan M. Newmark}, + Journal = {Journal of the Engineering Mechanics Division}, + Year = {1959}, + + Month = {July}, + Number = {3}, + Pages = {67-94}, + Volume = {85}, + + Editor = {ASCE} +} + +@Article{nguyen2001, + Title = {A cohesive model of fatigue crack growth}, + Author = {Nguyen, O. and Repetto, E. A. and Ortiz, M. and Radovitzky, R. A.}, + Journal = {International Journal of Fracture}, + Year = {2001}, + + Month = aug, + Number = {4}, + Pages = {351--369}, + Volume = {110}, + + Doi = {10.1023/A:1010839522926}, + ISSN = {0376-9429, 1573-2673}, + Language = {en}, + Urldate = {2015-02-17} +} + +@TechReport{Omohundro:1989, + Title = {Five Balltree Construction Algorithms}, + Author = {Stephen M. Omohundro}, + Institution = {International Computer Science Institute, University of California at Berkeley}, + Year = {1989}, + Number = {TR-89-063}, + Type = {Technical Report}, + + Bdsk-url-1 = {http://ftp.icsi.berkeley.edu/ftp/pub/techreports/1989/tr-89-063.pdf}, + Date-added = {2014-05-22 09:47:58 +0000}, + Date-modified = {2014-05-22 09:47:58 +0000}, + Url = {http://ftp.icsi.berkeley.edu/ftp/pub/techreports/1989/tr-89-063.pdf} +} + +@Article{ortiz1999, + Title = {Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis}, + Author = {Ortiz, M. and Pandolfi, A.}, + Journal = {International Journal for Numerical Methods in Engineering (IJNME)}, + Year = {1999}, + Pages = {1267-1282}, + Volume = {44} +} + +@Article{ortiz99a, + Title = {Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis}, + Author = {M Ortiz and A Pandolfi}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {1999}, + Pages = {1267-1282}, + Volume = {44} +} + +@Article{pandolfi12a, + Title = {An eigenerosion approach to brittle fracture}, + Author = {A Pandolfi and M Ortiz}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {2012}, + Pages = {694-714}, + Volume = {92} +} + +@Article{patzak01a, + Title = {Parallel explicit finite element dynamics with nonlocal constitutive models}, + Author = {B. Patz{\'a}k and D. Rypl and Z. Bittnar}, + Journal = {{Computers \& Structures}}, + Year = {2001}, + Number = {26--28}, + Pages = {2287 - 2297}, + Volume = {79}, + + ISSN = {0045-7949} +} + +@PhdThesis{Pietrzak:1997, + Title = {Continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems}, + Author = {Pietrzak, Grzegorz}, + School = {{\'E}cole {P}olytechnique {F}{\'e}d{\'e}rale de {L}ausanne}, + Year = {1997}, + + Date-added = {2014-09-16 12:23:01 +0000}, + Date-modified = {2014-09-16 12:41:09 +0000}, + Doi = {10.5075/epfl-thesis-1656} +} + +@Article{pijaudier87a, + Title = {Nonlocal damage theory}, + Author = {Gilles Pijaudier-Cabot and Zden\v{e}k P. Ba\v{z}ant}, + Journal = {Journal of Engineering Mechanics}, + Year = {1987}, + Number = {10}, + Pages = {1512 - 1533}, + Volume = {113} +} + +@Article{rice76a, + Title = {The Localisation of Plastic Deformation}, + Author = {James R. Rice}, + Journal = {Theoretical and Applied Mechanics (14th International Congress on Theoretical and Applied Mechanics, Delft, 1976, ed. W.T. Koiter)}, + Year = {1976}, + Pages = {207 - 220}, + Volume = {1} +} + +@Article{sharon96a, + Title = {Microbranching instability and the dynamic fracture of brittle materials}, + Author = {Eran Sharon and Jay Fineberg}, + Journal = {Physical Review B}, + Year = {1996}, + Number = {10}, + Pages = {7128 - 7139}, + Volume = {54} +} + +@Article{silling00a, + Title = {Reformulation of elasticity theory for discontinuities and long-range forces}, + Author = {S A Silling}, + Journal = {Journal of the Mechanics and Physics of Solids}, + Year = {2000}, + Pages = {175-209}, + Volume = {48} +} + +@Book{simo92, + Title = {Computational Inelasticity}, + Author = {Simo, J.C. and Hughes, T.J.R.}, + Publisher = {Springer}, + Year = {1992} +} + +@Article{snozzi_cohesive_2013, + Title = {A cohesive element model for mixed mode loading with frictional contact capability}, + Author = {Snozzi, Leonardo and Molinari, Jean-Francois}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {2013}, + + Month = feb, + Number = {5}, + Pages = {510--526}, + Volume = {93} +} + +@Book{Belytschko:2000, + Title = {Nonlinear Finite Elements for Continua and Structures}, + Author = {Ted Belytschko, Wing Kam Liu, Brian Moran}, + Publisher = {Wiley}, + Year = {2000} +} + +@Article{devree95, + Title = {Comparison of nonlocal approaches in continuum damage mechanics}, + Author = {J.H.P. de Vree and W.A.M. Brekelmans and M.A.J. van Gils}, + Journal = {Computers \&; Structures}, + Year = {1995}, + Number = {4}, + Pages = {581 - 588}, + Volume = {55}, + + ISSN = {0045-7949}, + Optdoi = {10.1016/0045-7949(94)00501-S} +} + +@Unpublished{vocialta15, + Title = {3D dynamic fragmentation with parallel dynamic insertion of cohesive elements}, + Author = {M. Vocialta and N. Richart and J.-F. Molinari}, + Note = {Submitted to IJNME}, + Year = {2015} +} + +@Unpublished{wolff14a, + Title = {A non-local continuum damage approach to model dynamic crack branching}, + Author = {C. Wolff and N. Richart and J.-F. Molinari}, + Note = {Submitted to IJNME}, + Year = {2014} +} + +@Article{Zhou_Molinari_2004, + Title = {Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency}, + Author = {F. Zhou and J. F. Molinari}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {2004}, + + Timestamp = {2015.07.30} +} + +@Misc{abaqus, + Title = {Simulia ABAQUS FEA}, + + Bdsk-url-1 = {http://www.3ds.com/products-services/simulia/portfolio/abaqus/}, + Key = {Unified FEA}, + Url = {\url{http://www.3ds.com/products-services/simulia/portfolio/abaqus/}} +} + +@Misc{cmake, + Title = {CMake - Cross Platform Make}, + + Bdsk-url-1 = {http://www.cmake.org/}, + Url = {\url{http://www.cmake.org/}} +} + +@Misc{diana, + Title = {TNO DIANA}, + + Bdsk-url-1 = {http://tnodiana.com/content/DIANA}, + Key = {FEM}, + Url = {\url{http://tnodiana.com/content/DIANA}} +} + +@Misc{mayavi, + Title = {The MayaVi Data Visualizer}, + + Bdsk-url-1 = {http://mayavi.sourceforge.net/}, + Url = {\url{http://mayavi.sourceforge.net/}} +} + +@Misc{mumps, + Title = {MUMPS : a parallel sparse direct solver}, + + Bdsk-url-1 = {http://graal.ens-lyon.fr/MUMPS/}, + Key = {sparse matrix, direct solver, parallelisme}, + Url = {\url{http://graal.ens-lyon.fr/MUMPS/}} +} + +@Misc{numpy, + Title = {NumPy - Fundamental package for scientific computing with Python}, + + Bdsk-url-1 = {http://www.numpy.org/}, + Url = {\url{http://www.numpy.org/}} +} +@Misc{paraview, + Title = {ParaView - Open Source Scientific Visualization}, + + Bdsk-url-1 = {http://www.paraview.org/}, + Url = {\url{http://www.paraview.org/}} +} + +@Misc{scotch, + Title = {SCOTCH: Static Mapping, Graph, Mesh and Hypergraph Partitioning}, + + Bdsk-url-1 = {http://www.labri.fr/perso/pelegrin/scotch/}, + Url = {\url{http://www.labri.fr/perso/pelegrin/scotch/}} +} + +@Misc{visit, + Title = {VisIt Visualization Tool}, + + Bdsk-url-1 = {http://wci.llnl.gov/codes/visit/}, + Url = {\url{http://wci.llnl.gov/codes/visit/}} +} + diff --git a/doc/dev-doc/manual/new-constitutive-laws.rst b/doc/dev-doc/manual/new-constitutive-laws.rst new file mode 100644 index 000000000..efb4d2392 --- /dev/null +++ b/doc/dev-doc/manual/new-constitutive-laws.rst @@ -0,0 +1,375 @@ +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: + +.. 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(); + + + /* ---------------------------------------------------------------------- */ + #include "material.hh" + /* ---------------------------------------------------------------------- */ + + #ifndef __AKANTU_MATERIAL_XXX_HH__ + #define __AKANTU_MATERIAL_XXX_HH__ + + namespace akantu { + + class MaterialXXX : public Material { + + /// declare here the interface of your material + + }; + +In the header file the user also needs to declare all the members of the new +material. These include the parameters that a read from the +material input file, as well as any other material parameters that will be +computed during the simulation and internal variables. + + +In the following the example of adding a new damage material will be +presented. In this case the parameters in the material will consist of the +Young's modulus, the Poisson coefficient, the resistance to damage and the +damage threshold. The material will then from these values compute its Lamé +coefficients and its bulk modulus. Furthermore, the user has to add a new +internal variable ``damage`` in order to store the amount of damage at each +quadrature point in each step of the simulation. For this specific material the +member declaration inside the class will look as follows:: + + class LocalMaterialDamage : public Material { + + /// declare constructors/destructors here + + /// declare methods and accessors here + + /* -------------------------------------------------------------------- */ + /* Class Members */ + /* -------------------------------------------------------------------- */ + + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); + private: + + /// the young modulus + Real E; + + /// Poisson coefficient + Real nu; + + /// First Lame coefficient + Real lambda; + + /// Second Lame coefficient (shear modulus) + Real mu; + + /// resistance to damage + Real Yd; + + /// damage threshold + Real Sd; + + /// Bulk modulus + Real kpa; + + /// damage internal variable + InternalField damage; + + }; + +In order to enable to print the material parameters at any point in +the user's example file using the standard output stream by typing:: + + for (UInt m = 0; m < model.getNbMaterials(); ++m) + std::cout << model.getMaterial(m) << std::endl; + +the standard output stream operator has to be redefined. This should be done at the end of the header file:: + + class LocalMaterialDamage : public Material { + + /// declare here the interace of your material + + }: + /* ---------------------------------------------------------------------- */ + /* inline functions */ + /* ---------------------------------------------------------------------- */ + /// standard output stream operator + inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this) + { + _this.printself(stream); + return stream; + } + +However, the user still needs to register the material parameters that +should be printed out. The registration is done during the call of the +constructor. Like all definitions the implementation of the +constructor has to be written in the ``material_XXX.cc`` +file. However, the declaration has to be provided in the +``material_XXX.hh`` file:: + + class LocalMaterialDamage : public Material { + /* -------------------------------------------------------------------- */ + /* Constructors/Destructors */ + /* -------------------------------------------------------------------- */ + public: + + LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); + }; + +The user can now define the implementation of the constructor in the +``material_XXX.cc`` file:: + + /* ---------------------------------------------------------------------- */ + #include "local_material_damage.hh" + #include "solid_mechanics_model.hh" + + namespace akantu { + + /* ---------------------------------------------------------------------- */ + LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model, + const ID & id) : + Material(model, id), + damage("damage", *this) { + AKANTU_DEBUG_IN(); + + this->registerParam("E", E, 0., _pat_parsable, "Young's modulus"); + this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio"); + this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient"); + this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient"); + this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient"); + this->registerParam("Yd", Yd, 50., _pat_parsmod); + this->registerParam("Sd", Sd, 5000., _pat_parsmod); + + damage.initialize(1); + + AKANTU_DEBUG_OUT(); + } + +During the intializer list the reference to the model and the material id are +assigned and the constructor of the internal field is called. Inside the scope +of the constructor the internal values have to be initialized and the +parameters, that should be printed out, are registered with the function: +``registerParam``:: + + void registerParam(name of the parameter (key in the material file), + member variable, + default value (optional parameter), + access permissions, + description); + +The available access permissions are as follows: +- ``_pat_internal``: Parameter can only be output when the material is printed. +- ``_pat_writable``: User can write into the parameter. The parameter is output when the material is printed. +- ``_pat_readable``: User can read the parameter. The parameter is output when the material is printed. +- ``_pat_modifiable``: Parameter is writable and readable. +- ``_pat_parsable``: Parameter can be parsed, *i.e.* read from the input file. +- ``_pat_parsmod``: Parameter is modifiable and parsable. + +In order to implement the new constitutive law the user needs to +specify how the additional material parameters, that are not +defined in the input material file, should be calculated. Furthermore, +it has to be defined how stresses and the stable time step should be +computed for the new local material. In the case of implicit +simulations, in addition, the computation of the tangent stiffness needs +to be defined. Therefore, the user needs to redefine the following +functions of the parent material:: + + void initMaterial(); + + // for explicit and implicit simulations void + computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); + + // for implicit simulations + void computeTangentStiffness(const ElementType & el_type, + Array & tangent_matrix, + GhostType ghost_type = _not_ghost); + + // for explicit and implicit simulations + Real getStableTimeStep(Real h, const Element & element); + +In the following a detailed description of these functions is provided: + +- ``initMaterial``: This method is called after the material file is fully read + and the elements corresponding to each material are assigned. Some of the + frequently used constant parameters are calculated in this method. For + example, the Lam\'{e} constants of elastic materials can be considered as such + parameters. + +- ``computeStress``: In this method, the stresses are computed based on the + constitutive law as a function of the strains of the quadrature points. For + example, the stresses for the elastic material are calculated based on the + following formula: + + .. math:: + + \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon} + + Therefore, this method contains a loop on all quadrature points assigned to + the material using the two macros: + ``MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN`` and + ``MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END`` + + .. code:: + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type); + + // sigma <- f(grad_u) + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + + The strain vector in Akantu contains the values of :math:`\nabla \vec{u}`, + i.e. it is really the *displacement gradient*, + +- ``computeTangentStiffness``: This method is called when the tangent to the + stress-strain curve is desired (see Fig \ref {fig:smm:AL:K}). For example, + it is called in the implicit solver when the stiffness matrix for the + regular elements is assembled based on the following formula: + + .. math:: + \label{eqn:smm:constitutive_elasc} \mat{K } + =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}} + + Therefore, in this method, the ``tangent`` matrix (\mat{D}) is + computed for a given strain. + + The ``tangent`` matrix is a :math:`4^{th}` order tensor which is stored as + a matrix in Voigt notation. + + .. _fig:smm:AL:K: + .. figure:: figures/tangent.svg + :align: center + :width: 60% + + Tangent to the stress-strain curve. + +.. + \begin{figure}[!htb] + \begin{center} + \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf} + \caption{Tangent to the stress-strain curve.} + \label{fig:smm:AL:K} + \end{center} + \end{figure} + +- ``getCelerity``: The stability criterion of the explicit integration scheme + depend on the fastest wave celerity~\eqref{eqn:smm:explicit:stabletime}. This + celerity depend on the material, and therefore the value of this velocity + should be defined in this method for each new material. By default, the + fastest wave speed is the compressive wave whose celerity can be defined in ``getPushWaveSpeed``. + +Once the declaration and implementation of the new material has been +completed, this material can be used in the user's example by including the header file:: + + #include "material_XXX.hh" + +For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by +default, the materials are initialized inside the method +``initFull``. If a local material should be used instead, the +initialization of the material has to be postponed until the local +material is registered in the model. Therefore, the model is +initialized with the boolean for skipping the material initialization +equal to true:: + + /// model initialization + model.initFull(_analysis_method = _explicit_lumped_mass); + +Once the model has been initialized, the local material needs +to be registered in the model:: + + model.registerNewCustomMaterials("name_of_local_material"); + +Only at this point the material can be initialized:: + + model.initMaterials(); + +A full example for adding a new damage law can be found in +``examples/new_material``. + +Adding a New Non-Local Constitutive Law +``````````````````````````````````````` + +In order to add a new non-local material we first have to add the local +constitutive law in Akantu (see above). We can then add the non-local version +of the constitutive law by adding the two files (``material_XXX_non_local.hh`` +and ``material_XXX_non_local.cc``) where ``XXX`` is the name of the +corresponding local material. The new law must inherit from the two classes, +non-local parent class, such as the ``MaterialNonLocal`` class, and from the +local version of the constitutive law, *i.e.* ``MaterialXXX``. It is therefore +necessary to include the interface of those classes in the header file of your +custom material and indicate the inheritance in the declaration of the class:: + + /* ---------------------------------------------------------------------- */ + #include "material_non_local.hh" // the non-local parent + #include "material_XXX.hh" + /* ---------------------------------------------------------------------- */ + + #ifndef __AKANTU_MATERIAL_XXX_HH__ + #define __AKANTU_MATERIAL_XXX_HH__ + + namespace akantu { + + class MaterialXXXNonLocal : public MaterialXXX, + public MaterialNonLocal { + + /// declare here the interface of your material + + }; + +As members of the class we only need to add the internal fields to store the +non-local quantities, which are obtained from the averaging process:: + + /* -------------------------------------------------------------------------- */ + /* Class members */ + /* -------------------------------------------------------------------------- */ + protected: + InternalField grad_u_nl; + +The following four functions need to be implemented in the non-local material:: + + /// initialization of the material + void initMaterial(); + /// loop over all element and invoke stress computation + virtual void computeNonLocalStresses(GhostType ghost_type); + /// compute stresses after local quantities have been averaged + virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type) + /// compute all local quantities + void computeStress(ElementType el_type, GhostType ghost_type); + +In the intialization of the non-local material we need to register the local +quantity for the averaging process. In our example the internal field +*grad_u_nl* is the non-local counterpart of the gradient of the displacement +field (*grad_u_nl*):: + + void MaterialXXXNonLocal::initMaterial() { + MaterialXXX::initMaterial(); + MaterialNonLocal::initMaterial(); + /// register the non-local variable in the manager + this->model->getNonLocalManager().registerNonLocalVariable( + this->grad_u.getName(), + this->grad_u_nl.getName(), + spatial_dimension * spatial_dimension); + } + +The function to register the non-local variable takes as parameters the name of +the local internal field, the name of the non-local counterpart and the number +of components of the field we want to average. In the *computeStress* we now +need to compute all the quantities we want to average. We can then write a loop +for the stress computation in the function *computeNonLocalStresses* and then +provide the constitutive law on each integration point in the function +*computeNonLocalStress*. diff --git a/doc/dev-doc/manual/solidmechanicsmodel.rst b/doc/dev-doc/manual/solidmechanicsmodel.rst index 93b03ac09..a3cdd0745 100644 --- a/doc/dev-doc/manual/solidmechanicsmodel.rst +++ b/doc/dev-doc/manual/solidmechanicsmodel.rst @@ -1,1221 +1,848 @@ Solid Mechanics Model ===================== -The solid mechanics model is a specific implementation of the ``Model`` -interface dedicated to handle the equations of motion or equations of -equilibrium. The model is created for a given mesh. It will create its own -``FEEngine`` object to compute the interpolation, gradient, integration and -assembly operations. A :cpp:class:`SolidMechanicsModel -` object can simply be created like this:: +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``: -blocked_dofs +: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. -displacement +: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). -velocity +: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). -acceleration +: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). -external_force +:cpp:func:`external_force `: contains the external forces applied on the nodes (:math:`\vec{f}_{\st{ext}}` in the following). -internal_force +: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. Model Setup ----------- Setting Initial Conditions `````````````````````````` For a unique solution of the equations of motion, initial displacements and velocities for all degrees of freedom must be specified: .. math:: \vec{u}(t=0) & = \vec{u}_0\\ \dot{\vec u}(t=0) & = \vec{v}_0 The solid mechanics model can be initialized as follows:: model.initFull() This function initializes the internal arrays and sets them to zero. Initial displacements and velocities that are not equal to zero can be prescribed by running a loop over the total number of nodes. Here, the initial displacement in :math:`x`-direction and the initial velocity in :math:`y`-direction for all nodes is set to :math:`0.1` and :math:`1`, respectively:: auto & disp = model.getDisplacement(); auto & velo = model.getVelocity(); for (UInt node = 0; node < mesh.getNbNodes(); ++node) { disp(node, 0) = 0.1; velo(node, 1) = 1.; } +.. _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 - :width: 75% Problem domain :math:`\Omega` with boundary in three dimensions. The Dirchelet and the Neumann regions of the boundary are denoted with :math:`\Gamma_u` and :math:`\Gamma_t`, respecitvely. Different ways of imposing these boundary conditions exist. A basic way is to loop over nodes or elements at the boundary and apply local values. A more advanced method consists of using the notion of the boundary of the mesh. In the following both ways are presented. Starting with the basic approach, as mentioned, the Dirichlet boundary conditions can be applied by looping over the nodes and assigning the required values. :numref:`fig:smm:dirichlet_bc` shows a beam with a fixed support on the left side. On the right end of the beam, a load is applied. At the fixed support, the displacement has a given value. For this example, the displacements in both the :math:`x` and the :math:`y`-direction are set to zero. Implementing this displacement boundary condition is similar to the implementation of initial displacement conditions described above. However, in order to impose a displacement boundary condition for all time steps, the corresponding nodes need to be marked as boundary nodes using the function ``blocked``. While, in order to impose a load on the right side, the nodes are not marked. -The detail codes are shown as follows:: +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 (UInt node = 0; node < nb_nodes; ++node) { if(Math::are_float_equal(pos(node, _x), 0)) { blocked(node, _x) = true; // block dof in x-direction blocked(node, _y) = true; // block dof in y-direction disp(node, _x) = 0.; // fixed displacement in x-direction disp(node, _y) = 0.; // fixed displacement in y-direction } else if (Math::are_float_equal(pos(node, _y), 0)) { blocked(node, _x) = false; // unblock dof in x-direction forces(node, _x) = 10.; // force in x-direction } } .. _fig:smm:dirichlet_bc: .. figure:: figures/dirichlet.svg :align: center Beam with fixed support and load. For the more advanced approach, one needs the notion of a boundary in the mesh. Therefore, the boundary should be created before boundary condition functors can be applied. Generally the boundary can be specified from the mesh file or the geometry. For the first case, the function ``createGroupsFromMeshData`` is called. This function can read any types of mesh data which are provided in the mesh file. If the mesh file is created with Gmsh, the function takes one input strings which is either ``tag_0``, ``tag_1`` or ``physical_names``. The first two tags are assigned by Gmsh to each element which shows the physical group that they belong to. In Gmsh, it is also possible to consider strings for different groups of elements. These elements can be separated by giving a string ``physical_names`` to the function -``createGroupsFromMeshData``:: +``createGroupsFromMeshData`` - mesh.createGroupsFromMeshData("physical_names"). +.. 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 ``applyBC`` needs to be -called on a :cpp:class:`SolidMechanicsModel `. This -function gets a Dirichlet or Neumann functor and a string which specifies the -desired boundary on which the boundary conditions is to be applied. The functors -specify the type of conditions to apply. Three built-in functors for Dirichlet -exist: ``FlagOnly, FixedValue,`` and ``IncrementValue``. The functor -``FlagOnly`` is used if a point is fixed in a given direction. Therefore, the -input parameter to this functor is only the fixed direction. The ``FixedValue`` -functor is used when a displacement value is applied in a fixed direction. The -``IncrementValue`` applies an increment to the displacement in a given -direction. The following code shows the utilization of three functors for the -top, bottom and side surface of the mesh which were already defined in the Gmsh -file:: +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 -``FromTraction`` of Neumann boundary conditions is called. Otherwise, the -functor ``FromStress`` should be called which gets the stress tensor as an input -parameter:: +: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.}; + 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:: +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 Fig.~\ref{fig:smm:bc_and_ic}. It consists of a plate that is fixed with movable supports on the left and bottom side. On the right side, a traction, which increases linearly with the number of time steps, is applied. The initial displacement and velocity in :math:`x`-direction at all free nodes is zero and two respectively. .. _fig:smm:bc_and_ic: .. figure:: figures/bc_and_ic_example.svg :align: center :width: 75% Plate on movable supports. .. \begin{figure}[!htb] \centering \includegraphics[scale=0.8]{figures/bc_and_ic_example} \caption{Plate on movable supports.\label{fig:smm:bc_and_ic}} \end{figure} As it is mentioned in Section \ref{sect:common:groups}, node and element groups can be used to assign the boundary conditions. A generic example is given below with a Dirichlet boundary condition:: // create a node group NodeGroup & node_group = mesh.createNodeGroup("nodes_fix"); /* fill the node group with the nodes you want */ // create an element group using the existing node group mesh.createElementGroupFromNodeGroup("el_fix", "nodes_fix", spatial_dimension-1); // boundary condition can be applied using the element group name model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "el_fix"); Material Selector ````````````````` If the user wants to assign different materials to different -finite elements groups in \akantu, a material selector has to be -used. By default, \akantu assigns the first valid material in the +finite elements groups in ``Akantu``, a material selector has to be +used. By default, ``Akantu`` assigns the first valid material in the material file to all elements present in the model (regular continuum materials are assigned to the regular elements and cohesive materials are assigned to cohesive elements or element facets). -To assign different materials to specific elements, mesh data -information such as tag information or specified physical names can be -used. ``MeshDataMaterialSelector`` class uses this information to -assign different materials. With the proper physical name or tag name -and index, different materials can be assigned as demonstrated in the -examples below:: +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. ``MaterialSelector`` and ``DefaultMaterialSelector`` is -used to assign a material to regular elements by default. For the -regular elements, as in the example above, -``MeshDataMaterialSelector`` can be used to assign different -materials to different elements. +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 +Apart from the ``Akantu``'s default material selectors, users can always develop their own classes in the main code to tackle various multi-material assignment situations. Static Analysis --------------- The :cpp:class:`SolidMechanicsModel ` class can handle different analysis methods, the first one being presented is the static case. In this case, the equation to solve is .. math:: \mat{K} \vec{u} = \vec{f}_{\mathrm{ext}} + :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 ``initFull`` method by giving the -``SolidMechanicsModelOptions``. These options specify the type of analysis to be -performed and whether the materials should be initialized with ``initMaterials`` -or not:: +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 -``_static`` to the method. By default, the Boolean for no -initialization of the materials is set to false, so that they are -initialized during the ``initFull``. The method ``initFull`` -also initializes all appropriate vectors to zero. Once the model is -created and initialized, the boundary conditions can be set as -explained in Section~\ref{sect:smm:boundary}. Boundary conditions -will prescribe the external forces for some free degrees of freedom -:math:`\vec{f}_{\st{ext}}` and displacements for some others. At this point -of the analysis, the function -``solveStep``\index{SolidMechanicsModel!solveStep} can be called:: - - model.solveStep<_scm_newton_raphson_tangent_modified, - SolveConvergenceCriteria::_residual>(1e-4, 1); +: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 (\ref{eqn:smm:static}) is modified in order to apply a +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 Equation~(\ref{eqn:smm:static}) is -automatically assembled. A Newton-Raphson iteration is subsequently started, -:math:`\mat{K}` is updated according to the displacement computed at the -previous iteration and one loops until the forces are balanced -(``SolveConvergenceCriteria::_residual``), i.e. :math:`||\vec{r}|| < -\texttt{SolveConvergenceCriteria::_residual}`. One can also iterate until the -increment of displacement is zero (``SolveConvergenceCriteria::_increment``) -which also means that the equilibrium is found. For a linear elastic problem, -the solution is obtained in one iteration and therefore the maximum number of -iterations can be set to one. But for a non-linear case, one needs to iterate as -long as the norm of the residual exceeds the tolerance threshold and therefore -the maximum number of iterations has to be higher, e.g. :math:`100`:: - - model.solveStep<_scm_newton_raphson_tangent_modified, - SolveConvergenceCriteria::_residual>(1e-4, 100) +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`. .. \begin{figure}[!htb] \centering \includegraphics[width=.7\linewidth]{figures/static_analysis} \caption{Solution of the static analysis. Left: the initial condition, right: the solution (deformation magnified 50 times)} \label{fig:smm:implicit:static_solution} \end{figure} .. _fig:smm:implicit:static_solution: .. figure:: figures/static_analysis.png :align: center :width: 75% Solution of the static analysis. Left: the initial condition, right: the solution (deformation magnified 50 times). -Static implicit analysis with dynamic insertion of cohesive elements -```````````````````````````````````````````````````````````````````` - - -In order to solve problems with the extrinsic cohesive method in the -static implicit solution scheme, the function ``solveStepCohesive`` -has to be used:: - - model.solveStepCohesive<_scm_newton_raphson_tangent, - SolveConvergenceCriteria::_increment>(1e-13, error, - 25, false, 1e5, - true); - - -in which the arguments are: tolerance, error, max_iteration, -load_reduction, tol_increase_factor, do_not_factorize. This -function, first applies the Newton-Raphson procedure to solve the -problem. Then, it calls the method ``checkCohesiveStress`` to -check if cohesive elements have to be inserted. Since the approach is -implicit, only one element is added, the most stressed one (see -Section \ref{extrinsic_insertion}). After insertion, the -Newton-Raphson procedure is applied again to solve the same -incremental loading step, with the new inserted cohesive element. The -procedure loops in this way since no new cohesive elements have to be -inserted. At that point, the solution is saved, and the simulation -can advance to the next incremental loading step. In case the -convergence is not reached, the obtained solution is not saved and the -simulation return to the main file with the error given by the -solution saved in the argument of the function *error*. In this -way, the user can intervene in the simulation in order to find anyhow -convergence. A possibility is, for instance, to reduce the last -incremental loading step. The variable *load_reduction* can be -used to identify if the load has been already reduced or not. At the -same time, with the variable *tol_increase_factor* it is -possible to increase the tolerance by a factor defined by the user in -the main file, in order to accept a solution even with an error bigger -than the tolerance set at the beginning. It is possible to increase -the tolerance only in the phase of loading reduction, i.e., when -load_reduction = true. A not converged solution is never saved. In -case the convergence is not reached even after the loading reduction -procedure, the displacement field is not updated and remains the one -of the last converged incremental steps. Also, cohesive elements are -inserted only if convergence is reached. An example of the extrinsic -cohesive method in the static implicit solution scheme is presented in -``examples/cohesive_element/cohesive_extrinsic_implicit``. +.. + Static implicit analysis with dynamic insertion of cohesive elements + ```````````````````````````````````````````````````````````````````` + + In order to solve problems with the extrinsic cohesive method in the + static implicit solution scheme, the function ``solveStepCohesive`` + has to be used:: + + model.solveStepCohesive<_scm_newton_raphson_tangent, + SolveConvergenceCriteria::_increment>(1e-13, error, + 25, false, 1e5, + true); + + + in which the arguments are: tolerance, error, max_iteration, + load_reduction, tol_increase_factor, do_not_factorize. This + function, first applies the Newton-Raphson procedure to solve the + problem. Then, it calls the method ``checkCohesiveStress`` to + check if cohesive elements have to be inserted. Since the approach is + implicit, only one element is added, the most stressed one (see + Section \ref{extrinsic_insertion}). After insertion, the + Newton-Raphson procedure is applied again to solve the same + incremental loading step, with the new inserted cohesive element. The + procedure loops in this way since no new cohesive elements have to be + inserted. At that point, the solution is saved, and the simulation + can advance to the next incremental loading step. In case the + convergence is not reached, the obtained solution is not saved and the + simulation return to the main file with the error given by the + solution saved in the argument of the function *error*. In this + way, the user can intervene in the simulation in order to find anyhow + convergence. A possibility is, for instance, to reduce the last + incremental loading step. The variable *load_reduction* can be + used to identify if the load has been already reduced or not. At the + same time, with the variable *tol_increase_factor* it is + possible to increase the tolerance by a factor defined by the user in + the main file, in order to accept a solution even with an error bigger + than the tolerance set at the beginning. It is possible to increase + the tolerance only in the phase of loading reduction, i.e., when + load_reduction = true. A not converged solution is never saved. In + case the convergence is not reached even after the loading reduction + procedure, the displacement field is not updated and remains the one + of the last converged incremental steps. Also, cohesive elements are + inserted only if convergence is reached. An example of the extrinsic + cohesive method in the static implicit solution scheme is presented in + ``examples/cohesive_element/cohesive_extrinsic_implicit``. Dynamic Methods --------------- Different ways to solve the equations of motion are implemented in the solid mechanics model. The complete equations that should be solved are: -.. _eqn:equation-motion: -.. math:: - \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = - \vec{f}_{\mathrm{ext}} +.. 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. +equation in the static case, where :math:`\ddot{\vec{u}} = \dot{\vec{u}} = 0`. +Here the method to solve this equation in the general case will be presented. +For this purpose, a time discretization has to be specified. The most common +discretization method in solid mechanics is the Newmark-:math:`\beta` method, +which is also the default in ``Akantu``. -For the Newmark-:math:`\beta` method, (:numref:`eqn:equation-motion`) becomes a -system of three equations (see \cite{curnier92a} \cite{hughes-83a} for -more details): +For the Newmark-:math:`\beta` method, (:eq:`eqn-equation-motion`) becomes a +system of three equations (see :cite:`curnier92a,hughes-83a` for more details): -.. _eqn:finite-difference: .. math:: - \mat{M} \ddot{\vec{u}}_{n+1} + \mat{C}\dot{\vec{u}}_{n+1} + \mat{K} \vec{u}_{n+1} &={\vec{f}_{\st{ext}}}_{\, n+1} - \label{eqn:equation-motion-discret} \\ - \vec{u}_{n+1} &=\vec{u}_{n} + \left(1 - \alpha\right) \Delta t \dot{\vec{u}}_{n} + - \alpha \Delta t \dot{\vec{u}}_{n+1} + \left(\frac{1}{2} - - \alpha\right) \Delta t^2 - \ddot{\vec{u}}_{n} \label{eqn:finite-difference-1}\\ - \dot{\vec{u}}_{n+1} &= \dot{\vec{u}}_{n} + \left(1 - \beta\right) - \Delta t \ddot{\vec{u}}_{n} + \beta \Delta t - \ddot{\vec{u}}_{n+1} - + \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~(\ref{eqn:equation-motion-discret}) is the equation of motion discretized in space (finite-element discretization), and the equations above are discretized in both space and time (Newmark discretization). The :math:`\alpha` and :math:`\beta` parameters determine the stability and the accuracy of the algorithm. Classical values for :math:`\alpha` and :math:`\beta` are usually :math:`\beta = 1/2` for no numerical damping and :math:`0 < \alpha < 1/2`. -+-------------------+-------------+----------+ -|:math:`\alpha` |Method |Type | -| |(:math:`\beta| | -| |= 1/2`) | | -+-------------------+-------------+----------+ -|:math:`0` |central |explicit | -| |difference | | -+-------------------+-------------+----------+ -|:math:`\frac{1}{6}`|Fox-Goodwin |implicit | -| |(royal road) | | -+-------------------+-------------+----------+ -|:math:`\frac{1}{3}`|Linear |implicit | -| |acceleration | | -+-------------------+-------------+----------+ -|:math:`\frac{1}{2}`|Average |implicit | -| |acceleration | | -| |(trapeziodal | | -| |rule) | | -+-------------------+-------------+----------+ +.. 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, -(\ref{eqn:equation-motion-discret})-(\ref{eqn:finite-difference-2}) is +(: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:: +- *Predictor*: - \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} + .. math:: \vec{u}_{n+1}^{0} &= \vec{u}_{n} + \Delta t \dot{\vec{u}}_{n} + \frac{\Delta t^2}{2} \ddot{\vec{u}}_{n} \\ + \dot{\vec{u}}_{n+1}^{0} &= \dot{\vec{u}}_{n} + \Delta t \ddot{\vec{u}}_{n} \\ + \ddot{\vec{u}}_{n+1}^{0} &= \ddot{\vec{u}}_{n} + +- *Solve*: + + .. math:: \left(c \mat{M} + d \mat{C} + e \mat{K}_{n+1}^i\right) + \vec{w} &= {\vec{f}_{\st{ext}}}_{\,n+1} - {\vec{f}_{\st{int}}}_{\,n+1}^i - + \mat{C} \dot{\vec{u}}_{n+1}^i - \mat{M} \ddot{\vec{u}}_{n+1}^i\\ + &= \vec{r}_{n+1}^i + +- *Corrector*: + + .. math:: \ddot{\vec{u}}_{n+1}^{i+1} &= \ddot{\vec{u}}_{n+1}^{i} +c \vec{w} \\ + \dot{\vec{u}}_{n+1}^{i+1} &= \dot{\vec{u}}_{n+1}^{i} + d\vec{w} \\ + \vec{u}_{n+1}^{i+1} &= \vec{u}_{n+1}^{i} + e \vec{w} where :math:`i` is the Newton-Raphson iteration counter and :math:`c`, :math:`d` and :math:`e` are parameters depending on the method used to solve the equations -+--------------------+----------------------------+------------------------+----------------------+----------------------+ -| |:math:`\vec{w}` |:math:`e` |:math:`d` |:math:`c` | -| | | | | | -| | | | | | -+--------------------+----------------------------+------------------------+----------------------+----------------------+ -|in acceleration |:math:`\delta\ddot{\vec{u}}`|:math:`\alpha\beta\Delta|:math:`\beta\Delta |:math:`1` | -| | |t^2` |t` | | -| | | | | | -+--------------------+----------------------------+------------------------+----------------------+----------------------+ -|in velocity |:math:`\delta\dot{\vec{u}}` |:math:`\alpha\Delta t` |:math:`1` |:math:`\frac{1}{\beta | -| | | | |\Delta t}` | -| | | | | | -+--------------------+----------------------------+------------------------+----------------------+----------------------+ -|in displacement |:math:`\delta\vec{u}` |:math:`1` |:math:`\frac{1}{\alpha|:math:`\frac{1}{\alpha| -| | | |\Delta t}` |\beta \Delta t^2}` | -| | | | | | -+--------------------+----------------------------+------------------------+----------------------+----------------------+ +.. csv-table:: + :header: "", ":math:`\vec{w}`", ":math:`e`", ":math:`d`", ":math:`c`" -.. - \begin{center} - \begin{tabular}{lcccc} - \toprule - & :math:`\vec{w}` & :math:`e` & :math:`d` & :math:`c`\\ - \midrule - in acceleration &:math:` \delta\ddot{\vec{u}}` & :math:`\alpha \beta\Delta t^2` &:math:`\beta \Delta t` &:math:`1`\\ - in velocity & :math:` \delta\dot{\vec{u}}`& :math:`\alpha\Delta t` & :math:`1` & :math:`\frac{1}{\beta \Delta t}`\\ - in displacement &:math:`\delta\vec{u}` & :math:` 1` & :math:`\frac{1}{\alpha \Delta t}` & :math:`\frac{1}{\alpha \beta \Delta t^2}`\\ - \bottomrule - \end{tabular} - \end{center} + "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}.} + +.. 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 +``Akantu`` implements a trapezoidal rule by default. That is to say :math:`\alpha = 1/2` and :math:`\beta = 1/2` which is unconditionally stable. Therefore the value of the time step can be chosen arbitrarily within reason:: model.setTimeStep(time_step); Since the system has to be solved for a given amount of time steps, the method ``solveStep()``, (which has already been used in the static example in Section~\ref{sect:smm:static}), is called inside a time loop:: /// time loop Real time = 0.; auto & solver = model.getNonLinearSolver(); solver.set("max_iterations", 100); solver.set("threshold", 1e-12); solver.set("convergence_type", SolveConvergenceCriteria::_solution); for (UInt s = 1; time ` class. In the initialization, the explicit scheme is selected using the ``_explicit_lumped_mass`` constant:: SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _explicit_lumped_mass); .. note:: Writing ``model.initFull()`` or ``model.initFull();`` is equivalent to use the ``_explicit_lumped_mass`` keyword, as this is the default case. -The explicit time integration scheme implemented in \akantu uses a +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 +Newmark:math:`-\beta` equations with :math:`\beta=1/2` and :math:`\alpha=0`. In ``Akantu`` these computations at each time step are invoked by calling the function ``solveStep``:: for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) { model.solveStep(); } The method ``solveStep`` wraps the four following functions: - ``model.explicitPred()`` allows to compute the displacement field at :math:`t+1` and a part of the velocity field at :math:`t+1`, denoted by :math:`\vec{\dot{u}^{\st{p}}}_{n+1}`, which will be used later in the method ``model.explicitCorr()``. The equations are: .. math:: \vec{u}_{n+1} &= \vec{u}_{n} + \Delta t \vec{\dot{u}}_{n} + \frac{\Delta t^2}{2} \vec{\ddot{u}}_{n}\\ \vec{\dot{u}^{\st{p}}}_{n+1} &= \vec{\dot{u}}_{n} + \Delta t \vec{\ddot{u}}_{n} - ``model.updateResidual()`` and ``model.updateAcceleration()`` compute the acceleration increment :math:`\delta \vec{\ddot{u}}`: .. math:: \left(\mat{M} + \frac{1}{2} \Delta t \mat{C}\right) \delta \vec{\ddot{u}} = \vec{f_{\st{ext}}} - \vec{f}_{\st{int}\, n+1} - \mat{C} \vec{\dot{u}^{\st{p}}}_{n+1} - \mat{M} \vec{\ddot{u}}_{n} The internal force :math:`\vec{f}_{\st{int}\, n+1}` is computed from the displacement :math:`\vec{u}_{n+1}` based on the constitutive law. - ``model.explicitCorr()`` computes the velocity and acceleration fields at :math:`t+1`: .. math:: \vec{\dot{u}}_{n+1} &= \vec{\dot{u}^{\st{p}}}_{n+1} + \frac{\Delta t}{2} \delta \vec{\ddot{u}} \\ \vec{\ddot{u}}_{n+1} &= \vec{\ddot{u}}_{n} + \delta \vec{\ddot{u}} The use of an explicit time integration scheme is illustrated by the example: ``examples/explicit/explicit_dynamic.cc``. This example models the propagation of a wave in a steel beam. The beam and the applied displacement in the :math:`x` direction are shown in :numref:`fig:smm:explicit`. .. _fig:smm:explicit: .. figure:: figures/explicit.svg :align: center :width: 90% Numerical setup. .. \begin{figure}[!htb] \centering \begin{tikzpicture} \coordinate (c) at (0,2); \draw[shift={(c)},thick, color=blue] plot [id=x, domain=-5:5, samples=50] ({\x, {(40 * sin(0.1*pi*3*\x) * exp(- (0.1*pi*3*\x)*(0.1*pi*3*\x) / 4))}}); \draw[shift={(c)},-latex] (-6,0) -- (6,0) node[right, below] {:math:`x`}; \draw[shift={(c)},-latex] (0,-0.7) -- (0,1) node[right] {:math:`u`}; \draw[shift={(c)}] (-0.1,0.6) node[left] {:math:`A`}-- (1.5,0.6); \coordinate (l) at (0,0.6); \draw[shift={(0,-0.7)}] (-5, 0) -- (5,0) -- (5, 1) -- (-5, 1) -- cycle; \draw[shift={(l)}, latex-latex] (-5,0)-- (5,0) node [midway, above] {:math:`L`}; \draw[shift={(l)}] (5,0.2)-- (5,-0.2); \draw[shift={(l)}] (-5,0.2)-- (-5,-0.2); \coordinate (h) at (5.3,-0.7); \draw[shift={(h)}, latex-latex] (0,0)-- (0,1) node [midway, right] {:math:`h`}; \draw[shift={(h)}] (-0.2,1)-- (0.2,1); \draw[shift={(h)}] (-0.2,0)-- (0.2,0); \end{tikzpicture} \caption{Numerical setup \label{fig:smm:explicit}} \end{figure} The length and height of the beam are :math:`L={10}\textrm{m}` and :math:`h = {1}\textrm{m}`, respectively. The material is linear elastic, homogeneous and isotropic (density: \SI{7800}{\kilo\gram\per\cubic\metre}, Young's modulus: \SI{210}{\giga\pascal} and Poisson's ratio: :math:`0.3`). The imposed displacement follow a Gaussian function with a maximum amplitude of :math:`A = {0.01}\textrm{m}`. The potential, kinetic and total energies are computed. The safety factor is equal to :math:`0.8`. - \input{manual-constitutive-laws} - -Adding a New Constitutive Law ------------------------------ - -There are several constitutive laws in \akantu as described in the -previous Section~\ref{sect:smm:CL}. It is also possible to use a -user-defined material for the simulation. These materials are referred -to as local materials since they are local to the example of the user -and not part of the \akantu library. To define a new local material, -two files (``material_XXX.hh`` and ``material_XXX.cc``) have -to be provided where ``XXX`` is the name of the new material. The -header file ``material_XXX.hh`` defines the interface of your -custom material. Its implementation is provided in the -``material_XXX.cc``. The new law must inherit from the -``Material`` class or any other existing material class. It is -therefore necessary to include the interface of the parent material -in the header file of your local material and indicate the inheritance -in the declaration of the class:: - - /* ---------------------------------------------------------------------- */ - #include "material.hh" - /* ---------------------------------------------------------------------- */ - - #ifndef __AKANTU_MATERIAL_XXX_HH__ - #define __AKANTU_MATERIAL_XXX_HH__ - - namespace akantu { - - class MaterialXXX : public Material { - - /// declare here the interface of your material - - }; - -In the header file the user also needs to declare all the members of the new -material. These include the parameters that a read from the -material input file, as well as any other material parameters that will be -computed during the simulation and internal variables. - - -In the following the example of adding a new damage material will be -presented. In this case the parameters in the material will consist of the -Young's modulus, the Poisson coefficient, the resistance to damage and the -damage threshold. The material will then from these values compute its Lamé -coefficients and its bulk modulus. Furthermore, the user has to add a new -internal variable ``damage`` in order to store the amount of damage at each -quadrature point in each step of the simulation. For this specific material the -member declaration inside the class will look as follows:: - - class LocalMaterialDamage : public Material { - - /// declare constructors/destructors here - - /// declare methods and accessors here - - /* -------------------------------------------------------------------- */ - /* Class Members */ - /* -------------------------------------------------------------------- */ - - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); - private: - - /// the young modulus - Real E; - - /// Poisson coefficient - Real nu; - - /// First Lame coefficient - Real lambda; - - /// Second Lame coefficient (shear modulus) - Real mu; - - /// resistance to damage - Real Yd; - - /// damage threshold - Real Sd; - - /// Bulk modulus - Real kpa; - - /// damage internal variable - InternalField damage; - - }; - -In order to enable to print the material parameters at any point in -the user's example file using the standard output stream by typing:: - - for (UInt m = 0; m < model.getNbMaterials(); ++m) - std::cout << model.getMaterial(m) << std::endl; - -the standard output stream operator has to be redefined. This should be done at the end of the header file:: - - class LocalMaterialDamage : public Material { - - /// declare here the interace of your material - - }: - /* ---------------------------------------------------------------------- */ - /* inline functions */ - /* ---------------------------------------------------------------------- */ - /// standard output stream operator - inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this) - { - _this.printself(stream); - return stream; - } - -However, the user still needs to register the material parameters that -should be printed out. The registration is done during the call of the -constructor. Like all definitions the implementation of the -constructor has to be written in the ``material_XXX.cc`` -file. However, the declaration has to be provided in the -``material_XXX.hh`` file:: - - class LocalMaterialDamage : public Material { - /* -------------------------------------------------------------------- */ - /* Constructors/Destructors */ - /* -------------------------------------------------------------------- */ - public: - - LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); - }; - -The user can now define the implementation of the constructor in the -``material_XXX.cc`` file:: - - /* ---------------------------------------------------------------------- */ - #include "local_material_damage.hh" - #include "solid_mechanics_model.hh" - - namespace akantu { - - /* ---------------------------------------------------------------------- */ - LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model, - const ID & id) : - Material(model, id), - damage("damage", *this) { - AKANTU_DEBUG_IN(); - - this->registerParam("E", E, 0., _pat_parsable, "Young's modulus"); - this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio"); - this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient"); - this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient"); - this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient"); - this->registerParam("Yd", Yd, 50., _pat_parsmod); - this->registerParam("Sd", Sd, 5000., _pat_parsmod); - - damage.initialize(1); - - AKANTU_DEBUG_OUT(); - } - -During the intializer list the reference to the model and the material id are -assigned and the constructor of the internal field is called. Inside the scope -of the constructor the internal values have to be initialized and the -parameters, that should be printed out, are registered with the function: -``registerParam``:: - - void registerParam(name of the parameter (key in the material file), - member variable, - default value (optional parameter), - access permissions, - description); - -The available access permissions are as follows: -- ``_pat_internal``: Parameter can only be output when the material is printed. -- ``_pat_writable``: User can write into the parameter. The parameter is output when the material is printed. -- ``_pat_readable``: User can read the parameter. The parameter is output when the material is printed. -- ``_pat_modifiable``: Parameter is writable and readable. -- ``_pat_parsable``: Parameter can be parsed, *i.e.* read from the input file. -- ``_pat_parsmod``: Parameter is modifiable and parsable. - -In order to implement the new constitutive law the user needs to -specify how the additional material parameters, that are not -defined in the input material file, should be calculated. Furthermore, -it has to be defined how stresses and the stable time step should be -computed for the new local material. In the case of implicit -simulations, in addition, the computation of the tangent stiffness needs -to be defined. Therefore, the user needs to redefine the following -functions of the parent material:: - - void initMaterial(); - - // for explicit and implicit simulations void - computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); - - // for implicit simulations - void computeTangentStiffness(const ElementType & el_type, - Array & tangent_matrix, - GhostType ghost_type = _not_ghost); - - // for explicit and implicit simulations - Real getStableTimeStep(Real h, const Element & element); - -In the following a detailed description of these functions is provided: - -- ``initMaterial``: This method is called after the material file is fully read - and the elements corresponding to each material are assigned. Some of the - frequently used constant parameters are calculated in this method. For - example, the Lam\'{e} constants of elastic materials can be considered as such - parameters. - -- ``computeStress``: In this method, the stresses are computed based on the - constitutive law as a function of the strains of the quadrature points. For - example, the stresses for the elastic material are calculated based on the - following formula: - - .. math:: - - \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon} - - Therefore, this method contains a loop on all quadrature points assigned to - the material using the two macros: - ``MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN`` and - ``MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END`` - - .. code:: - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type); - - // sigma <- f(grad_u) - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - - The strain vector in Akantu contains the values of :math:`\nabla \vec{u}`, - i.e. it is really the *displacement gradient*, - -- ``computeTangentStiffness``: This method is called when the tangent to the - stress-strain curve is desired (see Fig \ref {fig:smm:AL:K}). For example, - it is called in the implicit solver when the stiffness matrix for the - regular elements is assembled based on the following formula: - - .. math:: - \label{eqn:smm:constitutive_elasc} \mat{K } - =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}} - - Therefore, in this method, the ``tangent`` matrix (\mat{D}) is - computed for a given strain. - - The ``tangent`` matrix is a :math:`4^{th}` order tensor which is stored as - a matrix in Voigt notation. - - .. _fig:smm:AL:K: - .. figure:: figures/tangent.svg - :align: center - :width: 60% - - Tangent to the stress-strain curve. - -.. - \begin{figure}[!htb] - \begin{center} - \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf} - \caption{Tangent to the stress-strain curve.} - \label{fig:smm:AL:K} - \end{center} - \end{figure} - -- ``getCelerity``: The stability criterion of the explicit integration scheme - depend on the fastest wave celerity~\eqref{eqn:smm:explicit:stabletime}. This - celerity depend on the material, and therefore the value of this velocity - should be defined in this method for each new material. By default, the - fastest wave speed is the compressive wave whose celerity can be defined in ``getPushWaveSpeed``. - -Once the declaration and implementation of the new material has been -completed, this material can be used in the user's example by including the header file:: - - #include "material_XXX.hh" - -For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by -default, the materials are initialized inside the method -``initFull``. If a local material should be used instead, the -initialization of the material has to be postponed until the local -material is registered in the model. Therefore, the model is -initialized with the boolean for skipping the material initialization -equal to true:: - - /// model initialization - model.initFull(_analysis_method = _explicit_lumped_mass); - -Once the model has been initialized, the local material needs -to be registered in the model:: - - model.registerNewCustomMaterials("name_of_local_material"); - -Only at this point the material can be initialized:: - - model.initMaterials(); - -A full example for adding a new damage law can be found in -``examples/new_material``. - -Adding a New Non-Local Constitutive Law -``````````````````````````````````````` - -In order to add a new non-local material we first have to add the local -constitutive law in Akantu (see above). We can then add the non-local version -of the constitutive law by adding the two files (``material_XXX_non_local.hh`` -and ``material_XXX_non_local.cc``) where ``XXX`` is the name of the -corresponding local material. The new law must inherit from the two classes, -non-local parent class, such as the ``MaterialNonLocal`` class, and from the -local version of the constitutive law, *i.e.* ``MaterialXXX``. It is therefore -necessary to include the interface of those classes in the header file of your -custom material and indicate the inheritance in the declaration of the class:: - - /* ---------------------------------------------------------------------- */ - #include "material_non_local.hh" // the non-local parent - #include "material_XXX.hh" - /* ---------------------------------------------------------------------- */ - - #ifndef __AKANTU_MATERIAL_XXX_HH__ - #define __AKANTU_MATERIAL_XXX_HH__ - - namespace akantu { - - class MaterialXXXNonLocal : public MaterialXXX, - public MaterialNonLocal { - - /// declare here the interface of your material - - }; - -As members of the class we only need to add the internal fields to store the -non-local quantities, which are obtained from the averaging process:: - - /* -------------------------------------------------------------------------- */ - /* Class members */ - /* -------------------------------------------------------------------------- */ - protected: - InternalField grad_u_nl; - -The following four functions need to be implemented in the non-local material:: - - /// initialization of the material - void initMaterial(); - /// loop over all element and invoke stress computation - virtual void computeNonLocalStresses(GhostType ghost_type); - /// compute stresses after local quantities have been averaged - virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type) - /// compute all local quantities - void computeStress(ElementType el_type, GhostType ghost_type); - -In the intialization of the non-local material we need to register the local -quantity for the averaging process. In our example the internal field -*grad_u_nl* is the non-local counterpart of the gradient of the displacement -field (*grad_u_nl*):: - - void MaterialXXXNonLocal::initMaterial() { - MaterialXXX::initMaterial(); - MaterialNonLocal::initMaterial(); - /// register the non-local variable in the manager - this->model->getNonLocalManager().registerNonLocalVariable( - this->grad_u.getName(), - this->grad_u_nl.getName(), - spatial_dimension * spatial_dimension); - } +.. include:: ./constitutive-laws.rst -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*. +.. include:: ./new-constitutive-laws.rst diff --git a/doc/dev-doc/reference.rst b/doc/dev-doc/reference.rst index afd68de73..bb7dc0415 100644 --- a/doc/dev-doc/reference.rst +++ b/doc/dev-doc/reference.rst @@ -1,10 +1,71 @@ .. _reference: Reference --------- -.. doxygennamespace:: akantu - :project: Akantu - :members: - :protected-members: - :private-members: +Common +`````` + +.. doxygenfunction:: akantu::initialize(const std::string &input_file, int &argc, char **&argv) +.. doxygenfunction:: akantu::initialize(int &argc, char **&argv) + +.. doxygentypedef:: akantu::UInt +.. doxygentypedef:: akantu::Int +.. doxygentypedef:: akantu::Real + +.. doxygenenum:: akantu::ElementType +.. doxygenenum:: akantu::ModelType +.. doxygenenum:: akantu::AnalysisMethod +.. doxygenenum:: akantu::SolveConvergenceCriteria + +.. doxygenclass:: akantu::ArrayBase +.. doxygenclass:: akantu::ArrayDataLayer +.. doxygenclass:: akantu::Array + +.. doxygenclass:: akantu::ElementTypeMapArray + +.. doxygenclass:: akantu::Vector +.. doxygenclass:: akantu::Matrix + +Mesh +```` +.. doxygenclass:: akantu::Mesh +.. doxygenclass:: akantu::FEEngine + +Models +`````` + +Common +...... + +.. doxygenclass:: akantu::BC::Dirichlet::FixedValue +.. doxygenclass:: akantu::BC::Dirichlet::FlagOnly +.. doxygenclass:: akantu::BC::Dirichlet::IncrementValue +.. doxygenclass:: akantu::BC::Neumann::FromStress +.. doxygenclass:: akantu::BC::Neumann::FromTraction +.. doxygenclass:: akantu::BoundaryCondition +.. doxygenclass:: akantu::BoundaryConditionFunctor +.. doxygenclass:: akantu::EventHandlerManager +.. doxygenclass:: akantu::Model +.. doxygenclass:: akantu::NonLocalManagerCallback + +Solvers +....... + +.. doxygenclass:: akantu::ModelSolver +.. doxygenclass:: akantu::DOFManager +.. doxygenclass:: akantu::NonLinearSolver +.. doxygenclass:: akantu::NonLinearSolverNewtonRaphson + +Solid Mechanics Model +..................... + +.. doxygenclass:: akantu::SolidMechanicsModel +.. doxygenclass:: akantu::SolidMechanicsModelOptions +.. doxygenclass:: akantu::MaterialSelector +.. doxygenclass:: akantu::MeshDataMaterialSelector +.. doxygenclass:: akantu::Material + +Synchronizers +````````````` +.. doxygenclass:: akantu::DataAccessor diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh index d85d525a6..6390db524 100644 --- a/src/common/aka_array.hh +++ b/src/common/aka_array.hh @@ -1,439 +1,442 @@ /** * @file aka_array.hh * * @author Till Junge * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jan 16 2018 * * @brief Array container for Akantu * This container differs from the std::vector from the fact it as 2 dimensions * a main dimension and the size stored per entries * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "aka_common.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ARRAY_HH_ #define AKANTU_ARRAY_HH_ namespace akantu { /// class that afford to store vectors in static memory // NOLINTNEXTLINE(cppcoreguidelines-special-member-functions) class ArrayBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit ArrayBase(const ID &id = "") : id(id) {} ArrayBase(const ArrayBase & other, const ID & id = "") { this->id = (id.empty()) ? other.id : id; } ArrayBase(ArrayBase && other) = default; ArrayBase & operator=(const ArrayBase & other) = default; ArrayBase & operator=(ArrayBase && other) noexcept = default; virtual ~ArrayBase() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the amount of space allocated in bytes virtual UInt getMemorySize() const = 0; // changed empty to match std::vector empty inline bool empty() const __attribute__((warn_unused_result)) { return size_ == 0; } /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the Size of the Array UInt size() const { return size_; } /// Get the number of components AKANTU_GET_MACRO(NbComponent, nb_component, UInt); /// Get the name of th array AKANTU_GET_MACRO(ID, id, const ID &); /// Set the name of th array AKANTU_SET_MACRO(ID, id, const ID &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the vector ID id; /// the size used UInt size_{0}; /// number of components UInt nb_component{1}; }; /* -------------------------------------------------------------------------- */ namespace { template struct IteratorHelper {}; template struct IteratorHelper<0, T> { using type = T; }; template struct IteratorHelper<1, T> { using type = Vector; }; template struct IteratorHelper<2, T> { using type = Matrix; }; template struct IteratorHelper<3, T> { using type = Tensor3; }; template using IteratorHelper_t = typename IteratorHelper::type; } // namespace /* -------------------------------------------------------------------------- */ /* Memory handling layer */ /* -------------------------------------------------------------------------- */ enum class ArrayAllocationType { _default, _pod, }; template struct ArrayAllocationTrait : public std::conditional_t< std::is_scalar::value, std::integral_constant, std::integral_constant> {}; /* -------------------------------------------------------------------------- */ template ::value> class ArrayDataLayer : public ArrayBase { public: using value_type = T; using reference = value_type &; using pointer_type = value_type *; using const_reference = const value_type &; public: ~ArrayDataLayer() override = default; /// Allocation of a new vector explicit ArrayDataLayer(UInt size = 0, UInt nb_component = 1, const ID & id = ""); /// Allocation of a new vector with a default value ArrayDataLayer(UInt size, UInt nb_component, const_reference value, const ID & id = ""); /// Copy constructor (deep copy) ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = ""); /// Copy constructor (deep copy) explicit ArrayDataLayer(const std::vector & vect); // copy operator ArrayDataLayer & operator=(const ArrayDataLayer & other); // move constructor ArrayDataLayer(ArrayDataLayer && other) noexcept = default; // move assign ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default; protected: // deallocate the memory virtual void deallocate() {} // allocate the memory virtual void allocate(UInt size, UInt nb_component); // allocate and initialize the memory virtual void allocate(UInt size, UInt nb_component, const T & value); public: /// append a tuple of size nb_component containing value inline void push_back(const_reference value); /// append a vector // inline void push_back(const value_type new_elem[]); /// append a Vector or a Matrix template