diff --git a/.gitignore b/.gitignore index 9948f2a58..3a53eaa50 100644 --- a/.gitignore +++ b/.gitignore @@ -1,20 +1,25 @@ build* .dir-locals.el TAGS third-party/*/ !third-party/cmake/* !third-party/akantu-iterators !third-party/iohelper *~ release .*.swp *.tar.gz *.tgz *.tbz *.tar.bz2 .idea __pycache__ .mailmap +paraview/* +*.vtu +*.pvd +*.pvtu +*.vtk compile_commands.json .clangd .iwyu.imp diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3ca043942..b0c336627 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,349 +1,350 @@ stages: - configure - build - check-warnings - test - deploy .docker_build: image: 'docker:19.03.11' stage: .pre services: - docker:19.03.11-dind variables: # Use TLS https://docs.gitlab.com/ee/ci/docker/using_docker_build.html#tls-enabled DOCKER_HOST: tcp://docker:2376 DOCKER_TLS_CERTDIR: "/certs" before_script: - docker login -u $CI_REGISTRY_USER -p $CI_REGISTRY_PASSWORD $CI_REGISTRY script: - cd test/ci/${IMAGE_NAME}/ - docker build -t registry.gitlab.com/akantu/akantu/${IMAGE_NAME} . - docker push registry.gitlab.com/akantu/akantu/${IMAGE_NAME} docker build:debian-testing: variables: IMAGE_NAME: debian:testing extends: .docker_build rules: - changes: - test/ci/debian:testing/Dockerfile docker build:ubuntu-lts: variables: IMAGE_NAME: ubuntu:lts extends: .docker_build rules: - changes: - test/ci/ubuntu:lts/Dockerfile .configure: stage: configure except: - tags variables: BLA_VENDOR: 'Generic' script: - cmake -E make_directory build - cd build - cmake -DAKANTU_COHESIVE_ELEMENT:BOOL=TRUE -DAKANTU_IMPLICIT:BOOL=TRUE -DAKANTU_PARALLEL:BOOL=TRUE -DAKANTU_STRUCTURAL_MECHANICS:BOOL=TRUE -DAKANTU_HEAT_TRANSFER:BOOL=TRUE -DAKANTU_DAMAGE_NON_LOCAL:BOOL=TRUE -DAKANTU_PHASE_FIELD:BOOL=TRUE -DAKANTU_PYTHON_INTERFACE:BOOL=TRUE + -DAKANTU_CONTACT_MECHANICS:BOOL=TRUE -DAKANTU_EXAMPLES:BOOL=TRUE -DAKANTU_BUILD_ALL_EXAMPLES:BOOL=TRUE -DAKANTU_TEST_EXAMPLES:BOOL=TRUE -DAKANTU_TESTS:BOOL=TRUE -DAKANTU_RUN_IN_DOCKER:BOOL=TRUE -DCMAKE_BUILD_TYPE:STRING=Coverage .. - cp compile_commmands.json .. artifacts: when: on_success paths: - build - compile_commands.json expire_in: 10h .build: stage: build script: - cmake --build build/src > >(tee -a ${output}-out.log) 2> >(tee -a ${output}-err.log >&2) - cmake --build build/python > >(tee -a ${output}-out.log) 2> >(tee -a ${output}-err.log >&2) - cmake --build build/test/ > >(tee -a ${output}-out.log) 2> >(tee -a ${output}-err.log >&2) - cmake --build build/examples > >(tee -a ${output}-out.log) 2> >(tee -a ${output}-err.log >&2) artifacts: when: on_success paths: - build/ #- ${output}-out.log - ${output}-err.log - compile_commands.json expire_in: 10h .tests: stage: test script: - cd build - ctest -T test --no-compress-output --timeout 1800 after_script: - cd build - tag=$(head -n 1 < Testing/TAG) - if [ -e Testing/${tag}/Test.xml ]; then - xsltproc -o ./juint.xml ${CI_PROJECT_DIR}/test/ci/ctest2junit.xsl Testing/${tag}/Test.xml; - fi - gcovr --xml --gcov-executable "${GCOV_EXECUTABLE}" --output coverage.xml --object-directory ${CI_PROJECT_DIR}/build --root ${CI_PROJECT_DIR} -s || true artifacts: when: always paths: - build/juint.xml - build/coverage.xml reports: junit: - build/juint.xml cobertura: - build/coverage.xml .analyse_build: stage: check-warnings script: - if [[ $(cat ${output}-err.log | grep warning -i) ]]; then - cat ${output}-err.log; - exit 1; - fi allow_failure: true artifacts: when: on_failure paths: - "$output-err.log" # ------------------------------------------------------------------------------ .cache_build: variables: CCACHE_BASEDIR: ${CI_PROJECT_DIR}/build CCACHE_DIR: ${CI_PROJECT_DIR}/.ccache CCACHE_NOHASHDIR: 1 CCACHE_COMPILERCHECK: content cache: key: ${output} policy: pull-push paths: - .ccache/ - third-party/google-test - third-party/pybind11 before_script: - ccache --zero-stats || true after_script: - ccache --show-stats || true # ------------------------------------------------------------------------------ .image_debian_testing: image: registry.gitlab.com/akantu/akantu/debian:testing .image_ubuntu_lts: image: registry.gitlab.com/akantu/akantu/ubuntu:lts # ------------------------------------------------------------------------------ .compiler_gcc: variables: CC: /usr/lib/ccache/gcc CXX: /usr/lib/ccache/g++ FC: gfortran GCOV_EXECUTABLE: gcov .compiler_clang: variables: CC: /usr/lib/ccache/clang CXX: /usr/lib/ccache/clang++ FC: gfortran GCOV_EXECUTABLE: llvm-cov gcov # ------------------------------------------------------------------------------ .debian_testing_gcc: variables: output: debian_testing_gcc extends: - .compiler_gcc - .image_debian_testing - .cache_build .debian_testing_clang: variables: output: debian_testing_clang extends: - .compiler_clang - .image_debian_testing - .cache_build .ubuntu_lts_gcc: variables: output: ubuntu_lts_gcc extends: - .compiler_gcc - .image_ubuntu_lts - .cache_build # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ configure:debian_testing_gcc: extends: - .debian_testing_gcc - .configure cache: policy: pull-push build:debian_testing_gcc: extends: - .debian_testing_gcc - .build dependencies: - configure:debian_testing_gcc test:debian_testing_gcc: extends: - .debian_testing_gcc - .tests dependencies: - build:debian_testing_gcc analyse_build:debian_testing_gcc: extends: - .debian_testing_gcc - .analyse_build dependencies: - build:debian_testing_gcc # ------------------------------------------------------------------------------ configure:debian_testing_clang: extends: - .debian_testing_clang - .configure cache: policy: pull-push build:debian_testing_clang: extends: - .debian_testing_clang - .build dependencies: - configure:debian_testing_clang test:debian_testing_clang: extends: - .debian_testing_clang - .tests dependencies: - build:debian_testing_clang analyse_build:debian_testing_clang: extends: - .debian_testing_clang - .analyse_build dependencies: - build:debian_testing_clang # ------------------------------------------------------------------------------ configure:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .configure cache: policy: pull-push build:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build dependencies: - configure:ubuntu_lts_gcc analyse_build:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .analyse_build dependencies: - build:ubuntu_lts_gcc test:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .tests dependencies: - build:ubuntu_lts_gcc # ------------------------------------------------------------------------------ code_quality: stage: test image: docker:19.03.12 allow_failure: true services: - docker:19.03.12-dind variables: DOCKER_DRIVER: overlay2 DOCKER_HOST: tcp://docker:2376 DOCKER_TLS_CERTDIR: "/certs" CODECLIMATE_DEV: 1 CODE_QUALITY_IMAGE: "registry.gitlab.com/gitlab-org/ci-cd/codequality:0.85.22" needs: [] script: - export SOURCE_CODE=$PWD - | # this is required to avoid undesirable reset of Docker image ENV variables being set on build stage function propagate_env_vars() { CURRENT_ENV=$(printenv) for VAR_NAME; do echo $CURRENT_ENV | grep "${VAR_NAME}=" > /dev/null && echo "--env $VAR_NAME " done } - docker pull --quiet "$CODE_QUALITY_IMAGE" - | - docker build -t codeclimate/codeclimate-clang-tidy test/ci/codeclimate/codeclimate-clang-tidy - | docker run \ $(propagate_env_vars \ SOURCE_CODE \ TIMEOUT_SECONDS \ CODECLIMATE_DEBUG \ CODECLIMATE_DEV \ REPORT_STDOUT \ REPORT_FORMAT \ ENGINE_MEMORY_LIMIT_BYTES \ ) \ --volume "$PWD":/code \ --volume /var/run/docker.sock:/var/run/docker.sock \ "$CODE_QUALITY_IMAGE" /code artifacts: paths: - gl-code-quality-report.json reports: codequality: gl-code-quality-report.json expire_in: 1 week dependencies: [] rules: - if: '$CODE_QUALITY_DISABLED' when: never - if: '$CI_COMMIT_TAG || $CI_COMMIT_BRANCH' # ------------------------------------------------------------------------------ pages: stage: deploy extends: - .debian_testing_gcc script: - cd build - cmake -DAKANTU_DOCUMENTATION=ON .. - cmake --build . -t sphinx-doc - mv doc/dev-doc/html ../public dependencies: - build:debian_testing_gcc artifacts: paths: - public only: - master diff --git a/doc/dev-doc/akantu.dox.j2 b/doc/dev-doc/akantu.dox.j2 index 6a0439dbe..7a944900d 100644 --- a/doc/dev-doc/akantu.dox.j2 +++ b/doc/dev-doc/akantu.dox.j2 @@ -1,47 +1,52 @@ PROJECT_NAME = Akantu PROJECT_NUMBER = {{ akantu_version }} STRIP_FROM_PATH = {{ akantu_source_path }} STRIP_FROM_INC_PATH = {{ akantu_source_path }} TAB_SIZE = 4 ALIASES = "rst=\verbatim embed:rst" \ "endrst=\endverbatim" QUIET = NO WARN_IF_UNDOCUMENTED = NO WARN_IF_DOC_ERROR = NO WARN_AS_ERROR = NO -INPUT = {{ akantu_source_path }}/src +INPUT = {{ akantu_source_path }}/src \ + {{ akantu_source_path }}/test/ci/includes_for_ci FILE_PATTERNS = *.c *.cc *.hh *.py -EXCLUDE = {{ akantu_source_path }}/src/common/aka_fwd.hh +EXCLUDE = {{ akantu_source_path }}/src/common/aka_fwd.hh \ + {{ akantu_source_path }}/src/io/dumper/dumpable_dummy.hh + {{ akantu_source_path }}/src/common/aka_config.hh.in RECURSIVE = YES EXAMPLE_PATH = {{ akantu_source_path }}/examples EXAMPLE_RECURSIVE = YES SOURCE_BROWSER = NO GENERATE_HTML = NO GENERATE_HTMLHELP = NO USE_MATHJAX = YES GENERATE_LATEX = NO GENERATE_XML = YES XML_OUTPUT = xml -ENABLE_PREPROCESSING = YES -MACRO_EXPANSION = YES -PREDEFINED = DOXYGEN DOXYGEN_INCLUDE_PATH = {{ akantu_source_path }}/src/common \ {{ akantu_source_path }}/src/fe_engine \ {{ akantu_source_path }}/src/mesh \ - {{ akantu_source_path }}/src/model -PREDEFINED = "AKANTU_GET_MACRO(name, value, type) = type get##name() const;" \ + {{ akantu_source_path }}/src/model \ + {{ akantu_source_path }}/test/ci/includes_for_ci +ENABLE_PREPROCESSING = YES +MACRO_EXPANSION = YES +PREDEFINED = DOXYGEN \ + __cplusplus=201402L \ + "AKANTU_GET_MACRO(name, value, type) = type get##name() const;" \ "AKANTU_GET_MACRO_NOT_CONST(name, value, type) = type get##name();" EXPAND_AS_DEFINED = __BEGIN_AKANTU__ \ __END_AKANTU__ \ __BEGIN_AKANTU_DUMPER__ \ __END_AKANTU_DUMPER__ \ AKANTU_SET_MACRO \ AKANTU_GET_MACRO_DEREF_PTR \ AKANTU_GET_MACRO_BY_ELEMENT_TYPE \ AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST COLLABORATION_GRAPH = NO UML_LOOK = YES TEMPLATE_RELATIONS = YES CALL_GRAPH = YES CALLER_GRAPH = YES LOOKUP_CACHE_SIZE = 0 diff --git a/doc/dev-doc/conf.py b/doc/dev-doc/conf.py index d7b77b805..12886ef12 100644 --- a/doc/dev-doc/conf.py +++ b/doc/dev-doc/conf.py @@ -1,334 +1,335 @@ # -*- 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 shutil import jinja2 import git import re import subprocess # -- General configuration --------------------------------------------------- # If your documentation needs a minimal Sphinx version, state it here. # # needs_sphinx = '1.0' # Number figures numfig = True # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ 'sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.ifconfig', 'sphinx.ext.viewcode', 'sphinxcontrib.bibtex', 'breathe', ] read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True' if read_the_docs_build: akantu_path = "." akantu_source_path = "../../" else: akantu_path = "@CMAKE_CURRENT_BINARY_DIR@" akantu_source_path = "@CMAKE_SOURCE_DIR@" # 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', 'manual/appendix/elements.rst', 'manual/appendix/material-parameters.rst', 'manual/constitutive-laws.rst', 'manual/new-constitutive-laws.rst'] # The name of the Pygments (syntax highlighting) style to use. pygments_style = 'sphinx' primary_domain = 'cpp' highlight_language = 'cpp' bibtex_bibfiles = ['manual/manual-bibliography.bib'] # -- Project information ----------------------------------------------------- project = 'Akantu' copyright = '2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)' + \ ' Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)' author = 'Nicolas Richart' with open(os.path.join(akantu_source_path, 'VERSION'), 'r') as fh: version_file = fh.readlines() file_release = version_file[0].strip() try: tag_prefix = 'v' git_repo = git.Repo(akantu_source_path) git_describe = git_repo.git.describe('--tags', '--dirty', '--always', '--long', '--match', '{}*'.format(tag_prefix)) print("GIT Describe: {}".format(git_describe)) # git describe to PEP404 version describe_matches = re.search( (r'^{}(?P.+?)' + r'(?:-(?P\d+)-g(?P[0-9a-f]+)' + r'(?:-(?Pdirty))?)?$').format(tag_prefix), git_describe) if describe_matches: describe_matches = describe_matches.groupdict() release = describe_matches['version'] if describe_matches['distance']: release += '.' if '+' in release else '+' release += '{distance}.{sha}'.format(**describe_matches) if describe_matches['dirty']: release += '.dirty' else: count = git_repo.git.rev_list('HEAD', '--count') describe_matches = re.search( (r'^(?P[0-9a-f]+)' + r'(?:-(?Pdirty))?$').format(tag_prefix), git_describe).groupdict() release = '{}.{}+{}'.format(file_release, count, describe_matches['sha']) except git.InvalidGitRepositoryError: release = file_release version = re.sub(r'^([0-9]+)\.([0-9+]).*', r'\1.\2', release) print('Release: {} - Version: {}'.format(release, version)) # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # if read_the_docs_build: html_theme = 'default' else: html_theme = 'sphinx_rtd_theme' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. # # html_theme_options = {} # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] html_logo = '_static/logo_only_akantu.svg' # Custom sidebar templates, must be a dictionary that maps document names # to template names. # # The default sidebars (for documents that don't match any pattern) are # defined by theme itself. Builtin themes are using these templates by # default: ``['localtoc.html', 'relations.html', 'sourcelink.html', # 'searchbox.html']``. # # html_sidebars = {} html_sidebars = { '**': [ 'relations.html', # needs 'show_related': True theme option to display 'searchbox.html', ]} math_eqref_format = "Eq. {number}" # MathJax configuration if not read_the_docs_build: mathjax_config = { 'extensions': [ "tex2jax.js", "siunitx.js" ], 'TeX': { 'Macros': { 'st': [r'\mathrm{#1}', 1], 'mat': [r'\mathbf{#1}', 1], 'half': [r'\frac{1}{2}', 0], }, 'extensions': ["AMSmath.js", "AMSsymbols.js", "sinuitx.js"], }, } else: mathjax3_config = { 'tex': { 'macros': { 'st': [r'\mathrm{#1}', 1], 'mat': [r'\mathbf{#1}', 1], 'half': [r'\frac{1}{2}', 0], }, 'packages': ['base', 'ams'], }, 'loader': { 'load': ['[tex]/ams'] }, } # -- Options for HTMLHelp output --------------------------------------------- # Output file base name for HTML help builder. htmlhelp_basename = 'Akantudoc' # -- Options for LaTeX output ------------------------------------------------ latex_elements = { # The paper size ('letterpaper' or 'a4paper'). # # 'papersize': 'letterpaper', # The font size ('10pt', '11pt' or '12pt'). # # 'pointsize': '10pt', # Additional stuff for the LaTeX preamble. # 'preamble': r'''\usepackage{amsmath}''', # Latex figure (float) alignment # # 'figure_align': 'htbp', } # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, # author, documentclass [howto, manual, or own class]). latex_documents = [ (master_doc, 'Akantu.tex', 'Akantu Documentation', 'Nicolas Richart', 'manual'), ] # -- Options for manual page output ------------------------------------------ # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). man_pages = [ (master_doc, 'akantu', 'Akantu Documentation', [author], 1) ] # -- Options for Texinfo output ---------------------------------------------- # Grouping the document tree into Texinfo files. List of tuples # (source start file, target name, title, author, # dir menu entry, description, category) texinfo_documents = [ (master_doc, 'Akantu', 'Akantu Documentation', author, 'Akantu', 'One line description of project.', 'Miscellaneous'), ] # -- Options for Epub output ------------------------------------------------- # Bibliographic Dublin Core info. epub_title = project epub_author = author epub_publisher = author epub_copyright = copyright # The unique identifier of the text. This can be a ISBN number # or the project homepage. # epub_identifier = '' # A unique identification for the text. # # epub_uid = '' # A list of files that should not be packed into the epub file. epub_exclude_files = ['search.html'] # -- Extension configuration ------------------------------------------------- j2_args = {} if read_the_docs_build: j2_template_path = '.' else: j2_template_path = '@CMAKE_CURRENT_SOURCE_DIR@' os.makedirs(os.path.join(akantu_path, '_static'), exist_ok=True) shutil.copyfile( os.path.join('@CMAKE_CURRENT_SOURCE_DIR@', html_logo), os.path.join(akantu_path, html_logo)) j2_args = { 'akantu_source_path': akantu_source_path, 'akantu_version': version.replace('v', ''), } print(akantu_path) j2_env = jinja2.Environment( loader=jinja2.FileSystemLoader(j2_template_path), undefined=jinja2.DebugUndefined) j2_template = j2_env.get_template('akantu.dox.j2') with open(os.path.join(akantu_path, 'akantu.dox'), 'w') as fh: fh.write(j2_template.render(j2_args)) subprocess.run(['doxygen', 'akantu.dox'], cwd=akantu_path) # print("akantu_path = '{}'".format(akantu_path)) breathe_projects = {"Akantu": os.path.join(akantu_path, "xml")} breathe_default_project = "Akantu" breathe_default_members = ('members', 'undoc-members') breathe_implementation_filename_extensions = ['.c', '.cc', '.cpp'] breathe_show_enumvalue_initializer = True +breathe_debug_trace_directives = 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 34de864af..35c8ab203 100644 --- a/doc/dev-doc/index.rst +++ b/doc/dev-doc/index.rst @@ -1,50 +1,51 @@ .. 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 ./manual/heattransfermodel.rst + ./manual/contactmechanicsmodel.rst ./manual/phasefieldmodel.rst ./manual/structuralmechanicsmodel.rst ./manual/io.rst .. toctree:: :maxdepth: 2 :caption: Changelog ./changelog.rst .. toctree:: :maxdepth: 2 :caption: API Reference ./reference.rst .. toctree:: :maxdepth: 2 :caption: Appendix ./manual/appendix.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/constitutive-laws.rst b/doc/dev-doc/manual/constitutive-laws.rst index e7b3ff616..6b76d5957 100644 --- a/doc/dev-doc/manual/constitutive-laws.rst +++ b/doc/dev-doc/manual/constitutive-laws.rst @@ -1,1409 +1,1407 @@ .. _sect-smm-cl: - - Constitutive Laws ----------------- In order to compute an element’s response to deformation, one needs to use an appropriate constitutive relationship. The constitutive law is used to compute the element’s stresses from the element’s strains. In the finite-element discretization, the constitutive formulation is applied to every quadrature point of each element. When the implicit formulation is used, the tangent matrix has to be computed. | The chosen materials for the simulation have to be specified in the mesh file or, as an alternative, they can be assigned using the at ``element_material`` vector. For every material assigned to the problem one has to specify the material characteristics (constitutive behavior and material properties) using the text input file (see :ref:`sect-io-material`). | In order to conveniently store values at each quadrature in a material point ``Akantu`` provides a special data structure, the at :cpp:class:`InternalField `. The internal fields are inheriting from the at :cpp:class:`ElementTypeMapArray `. Furthermore, it provides several functions for initialization, auto-resizing and auto removal of quadrature points. The constitutive law is precised within an input file. The dedicated material section is then read by :cpp:func:`initFull ` method of :cpp:class:`SolidMechanicsModel ` which initializes the different materials specified with the following convention .. code-block:: python material constitutive_law [ name = value rho = value ... ] where *constitutive_law* is the adopted constitutive law, followed by the material properties listed one by line in the bracket (*e.g.*, ``name`` and density :math:``rho``. Some constitutive laws can also have an *optional flavor*. For certain materials, it is possible to activate the large deformation strain and stress evaluations. Internally the strain measure becomes the right Cauchy–Green deformation tensor and the evaluated stress measure becomes the Piola-Kirchhoff stress tensor. This is activated by using: .. code-block:: python material constitutive_law [ finite_deformation = true # Activates the large deformation routines (bool) ... ] Sometimes it is also desired to generate random distributions of internal parameters. An example might be the critical stress at which the material fails. To generate such a field, in the text input file, a random quantity needs be added to the base value: All parameters are real numbers. For the uniform distribution, minimum and maximum values have to be specified. Random parameters are defined as a :math:`base` value to which we add a random number that follows the chosen distribution. The `Uniform `__ distribution is gives a random values between in :math:`[min, max)`. The `Weibull `__ distribution is characterized by the following cumulative distribution function: .. math:: F(x) = 1- e^{-\left({x/\lambda}\right)^m} which depends on :math:`m` and :math:`\lambda`, which are the shape parameter and the scale parameter. These random distributions are different each time the code is executed. In order to obtain always the same one, it possible to manually set the *seed* that is the number from which these pseudo-random distributions are created. This can be done by adding the following line to the input file *outside* the material parameters environments: .. code-block:: seed = 1.0 where the value 1.0 can be substituted with any number. Currently ``Akantu`` can reproduce always the same distribution when the seed is specified *only* in serial. The value of the *seed* can be also specified directly in the code (for instance in the main file) with the command: .. code-block:: RandGenerator::seed(1.0) The same command, with empty brackets, can be used to check the value of the *seed* used in the simulation. The following sections describe the constitutive models implemented in ``Akantu``. ----- Elastic ``````` The elastic law is a commonly used constitutive relationship that can be used for a wide range of engineering materials (*e.g.*, metals, concrete, rock, wood, glass, rubber, etc.) provided that the strains remain small (*i.e.*, small deformation and stress lower than yield strength). The elastic laws are often expressed as :math:`\boldsymbol{\sigma} = \boldsymbol{C}:\boldsymbol{\varepsilon}` with where :math:`\boldsymbol{\sigma}` is the Cauchy stress tensor, :math:`\boldsymbol{\varepsilon}` represents the infinitesimal strain tensor and :math:`\boldsymbol{C}` is the elastic modulus tensor. .. _sect-smm-linear-elastic-isotropic: Linear isotropic '''''''''''''''' Keyword: **elastic** Material description with input file: .. code-block:: python #input.dat material elastic [ name = steel rho = 7800 # density (Real) E = 2.1e11 # young's modulus (Real) nu = 0.3 # poisson's ratio (Real) Plane_stress = false # Plane stress simplification (only 2D problems) (bool) finite_deformation = false # activates the evaluation of strains with green's tensor (bool) ] Available energies Energies: - ``potential``: elastic potential energy ----- The linear isotropic elastic behavior is described by Hooke’s law, which states that the stress is linearly proportional to the applied strain (material behaves like an ideal spring), as illustrated in  :numref:`fig:smm:cl:el`. .. figure:: figures/cl/stress_strain_el.svg :alt: Elastic strain-stress curve :name: fig:smm:cl:el :width: 60.0% Stress-strain curve of elastic material and schematic representation of Hooke's law, denoted as a spring. The equation that relates the strains to the displacements is: point) from the displacements as follows: .. math:: \label{eqn:smm:strain_inf} \boldsymbol{\varepsilon} = \frac{1}{2} \left[ \nabla_0 \boldsymbol{u}+\nabla_0 \boldsymbol{u}^T \right] where :math:`\boldsymbol{\varepsilon}` represents the infinitesimal strain tensor, :math:`\nabla_{0}\boldsymbol{u}` the displacement gradient tensor according to the initial configuration. The constitutive equation for isotropic homogeneous media can be expressed as: .. math:: \label{eqn:smm:material:constitutive_elastic} \boldsymbol{\sigma } =\lambda\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{I}+2 \mu\boldsymbol{\varepsilon} where :math:`\boldsymbol{\sigma}` is the Cauchy stress tensor (:math:`\lambda` and :math:`\mu` are the the first and second Lame’s coefficients). In Voigt notation this correspond to .. math:: \begin{aligned} \left[\begin{array}{c} \sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{23}\\ \sigma_{13}\\ \sigma_{12}\\ \end{array}\right] &= \frac{E}{(1+\nu)(1-2\nu)}\left[ \begin{array}{cccccc} 1-\nu & \nu & \nu & 0 & 0 & 0\\ \nu & 1-\nu & \nu & 0 & 0 & 0\\ \nu & \nu & 1-\nu & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \\ \end{array}\right] \left[\begin{array}{c} \varepsilon_{11}\\ \varepsilon_{22}\\ \varepsilon_{33}\\ 2\varepsilon_{23}\\ 2\varepsilon_{13}\\ 2\varepsilon_{12}\\ \end{array}\right]\end{aligned} This formulation is not sufficient to represent all elastic material behavior. Some materials have characteristic orientation that have to be taken into account. To represent this anisotropy a more general stress-strain law has to be used, as shown below. ----- .. _sect-smm-linear-elastic-anisotropic: Linear anisotropic '''''''''''''''''' Keyword: **elastic_anisotropic** Material description with input file: .. code-block:: python #input.dat material elastic_anisotropic [ name = aluminum rho = 1.6465129043044597 # density (Real) C11 = 105.092023 # Coefficient ij of material tensor C (Real) C12 = 59.4637759 # all the 36 values C13 = 59.4637759 # in Voigt notation can be entered C14 = 0 # zero coefficients can be omited C15 = 0 C16 = 0 C22 = 105.092023 C23 = 59.4637759 C24 = 0 C25 = 0 C26 = 0 C33 = 105.092023 C34 = 0 C35 = 0 C36 = 0 C44 = 30.6596356 C45 = 0 C46 = 0 C55 = 30.6596356 C56 = 0 C66 = 30.6596356 n1 = [-1, 1, 0] # Direction of first material axis (Vector) n2 = [ 1, 1, 1] # Direction of second material axis (Vector) n3 = [ 1, 1, -2] # Direction of thrid material axis (Vector) ] ---- We define the elastic modulus tensor as follows: .. math:: \begin{aligned} \left[\begin{array}{c} \sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{23}\\ \sigma_{13}\\ \sigma_{12}\\ \end{array}\right] &= \left[ \begin{array}{cccccc} c_{11} & c_{12} & c_{13} & c_{14} & c_{15} & c_{16}\\ c_{21} & c_{22} & c_{23} & c_{24} & c_{25} & c_{26}\\ c_{31} & c_{32} & c_{33} & c_{34} & c_{35} & c_{36}\\ c_{41} & c_{42} & c_{43} & c_{44} & c_{45} & c_{46}\\ c_{51} & c_{52} & c_{53} & c_{54} & c_{55} & c_{56}\\ c_{61} & c_{62} & c_{63} & c_{64} & c_{65} & c_{66}\\ \end{array}\right] \left[\begin{array}{c} \varepsilon_{11}\\ \varepsilon_{22}\\ \varepsilon_{33}\\ 2\varepsilon_{23}\\ 2\varepsilon_{13}\\ 2\varepsilon_{12}\\ \end{array}\right]\end{aligned} To simplify the writing of input files the :math:`\boldsymbol{C}` tensor is expressed in the material basis. And this basis as to be given too. This basis :math:`\Omega_{{\mathrm{mat}}} = \{\boldsymbol{n_1}, \boldsymbol{n_2}, \boldsymbol{n_3}\}` is used to define the rotation :math:`R_{ij} = \boldsymbol{n_j} . \boldsymbol{e_i}`. And :math:`\boldsymbol{C}` can be rotated in the global basis :math:`\Omega = \{\boldsymbol{e_1}, \boldsymbol{e_2}, \boldsymbol{e_3}\}` as follow: .. math:: \begin{aligned} \boldsymbol{C}_{\Omega} &= \boldsymbol{R}_1 \boldsymbol{C}_{\Omega_{{\mathrm{mat}}}} \boldsymbol{R}_2\\ \boldsymbol{R}_1 &= \left[ \begin{array}{cccccc} R_{11} R_{11} & R_{12} R_{12} & R_{13} R_{13} & R_{12} R_{13} & R_{11} R_{13} & R_{11} R_{12}\\ R_{21} R_{21} & R_{22} R_{22} & R_{23} R_{23} & R_{22} R_{23} & R_{21} R_{23} & R_{21} R_{22}\\ R_{31} R_{31} & R_{32} R_{32} & R_{33} R_{33} & R_{32} R_{33} & R_{31} R_{33} & R_{31} R_{32}\\ R_{21} R_{31} & R_{22} R_{32} & R_{23} R_{33} & R_{22} R_{33} & R_{21} R_{33} & R_{21} R_{32}\\ R_{11} R_{31} & R_{12} R_{32} & R_{13} R_{33} & R_{12} R_{33} & R_{11} R_{33} & R_{11} R_{32}\\ R_{11} R_{21} & R_{12} R_{22} & R_{13} R_{23} & R_{12} R_{23} & R_{11} R_{23} & R_{11} R_{22}\\ \end{array}\right]\\ \boldsymbol{R}_2 &= \left[ \begin{array}{cccccc} R_{11} R_{11} & R_{21} R_{21} & R_{31} R_{31} & R_{21} R_{31} & R_{11} R_{31} & R_{11} R_{21}\\ R_{12} R_{12} & R_{22} R_{22} & R_{32} R_{32} & R_{22} R_{32} & R_{12} R_{32} & R_{12} R_{22}\\ R_{13} R_{13} & R_{23} R_{23} & R_{33} R_{33} & R_{23} R_{33} & R_{13} R_{33} & R_{13} R_{23}\\ R_{12} R_{13} & R_{22} R_{23} & R_{32} R_{33} & R_{22} R_{33} & R_{12} R_{33} & R_{12} R_{23}\\ R_{11} R_{13} & R_{21} R_{23} & R_{31} R_{33} & R_{21} R_{33} & R_{11} R_{33} & R_{11} R_{23}\\ R_{11} R_{12} & R_{21} R_{22} & R_{31} R_{32} & R_{21} R_{32} & R_{11} R_{32} & R_{11} R_{22}\\ \end{array}\right]\\\end{aligned} ----- .. _sect-smm-linear-elastic-orthotropic: Linear orthotropic '''''''''''''''''' Keyword: **elastic_orthotropic** Inherits from **elastic_anisotropic** Material description with input file: .. code-block:: python #input.dat material elastic_orthotropic [ name = test_mat_1 rho = 1 # density n1 = [-1, 1, 0] # Direction of first material axis (Vector) n2 = [ 1, 1, 1] # Direction of second material axis (Vector) n3 = [ 1, 1, -2] # Direction of thrid material axis (Vector) E1 = 1 # Young's modulus in direction n1 (Real) E2 = 2 # Young's modulus in direction n2 (Real) E3 = 3 # Young's modulus in direction n3 (Real) nu12 = 0.1 # Poisson's ratio 12 (Real) nu13 = 0.2 # Poisson's ratio 13 (Real) nu23 = 0.3 # Poisson's ratio 23 (Real) G12 = 0.5 # Shear modulus 12 (Real) G13 = 1 # Shear modulus 13 (Real) G23 = 2 # Shear modulus 23 (Real) ] ----- A particular case of anisotropy is when the material basis is orthogonal in which case the elastic modulus tensor can be simplified and rewritten in terms of 9 independents material parameters. .. math:: \begin{aligned} \left[\begin{array}{c} \sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{23}\\ \sigma_{13}\\ \sigma_{12}\\ \end{array}\right] &= \left[ \begin{array}{cccccc} c_{11} & c_{12} & c_{13} & 0 & 0 & 0 \\ & c_{22} & c_{23} & 0 & 0 & 0 \\ & & c_{33} & 0 & 0 & 0 \\ & & & c_{44} & 0 & 0 \\ & \text{sym.} & & & c_{55} & 0 \\ & & & & & c_{66}\\ \end{array}\right] \left[\begin{array}{c} \varepsilon_{11}\\ \varepsilon_{22}\\ \varepsilon_{33}\\ 2\varepsilon_{23}\\ 2\varepsilon_{13}\\ 2\varepsilon_{12}\\ \end{array}\right]\end{aligned} .. math:: \begin{aligned} c_{11} &= E_1 (1 - \nu_{23}\nu_{32})\Gamma \qquad c_{22} = E_2 (1 - \nu_{13}\nu_{31})\Gamma \qquad c_{33} = E_3 (1 - \nu_{12}\nu_{21})\Gamma\\ c_{12} &= E_1 (\nu_{21} - \nu_{31}\nu_{23})\Gamma = E_2 (\nu_{12} - \nu_{32}\nu_{13})\Gamma\\ c_{13} &= E_1 (\nu_{31} - \nu_{21}\nu_{32})\Gamma = E_2 (\nu_{13} - \nu_{21}\nu_{23})\Gamma\\ c_{23} &= E_2 (\nu_{32} - \nu_{12}\nu_{31})\Gamma = E_3 (\nu_{23} - \nu_{21}\nu_{13})\Gamma\\ c_{44} &= \mu_{23} \qquad c_{55} = \mu_{13} \qquad c_{66} = \mu_{12} \\ \Gamma &= \frac{1}{1 - \nu_{12} \nu_{21} - \nu_{13} \nu_{31} - \nu_{32} \nu_{23} - 2 \nu_{21} \nu_{32} \nu_{13}}\end{aligned} The Poisson ratios follow the rule :math:`\nu_{ij} = \nu_{ji} E_i / E_j`. ----- .. _sect-smm-cl-neohookean: Neo-Hookean ''''''''''' Keyword: **neohookean** Inherits from **elastic** Material description with input file: .. code-block:: python #input.dat material neohookean [ name = material_name rho = 7800 # density (Real) E = 2.1e11 # young's modulus (Real) nu = 0.3 # poisson's ratio (Real) ] ----- The hyperelastic Neo-Hookean constitutive law results from an extension of the linear elastic relationship (Hooke’s Law) for large deformation. Thus, the model predicts nonlinear stress-strain behavior for bodies undergoing large deformations. .. figure:: figures/cl/stress_strain_neo.svg :alt: Neo-hookean Stress-strain curve. :name: fig:smm:cl:neo_hookean :width: 40.0% Neo-hookean Stress-strain curve. As illustrated in :numref:`fig:smm:cl:neo_hookean`, the behavior is initially linear and the mechanical behavior is very close to the corresponding linear elastic material. This constitutive relationship, which accounts for compressibility, is a modified version of the one proposed by Ronald Rivlin :cite:`Belytschko:2000`. The strain energy stored in the material is given by: .. math:: \label{eqn:smm:constitutive:neohookean_potential} \Psi(\boldsymbol{C}) = \frac{1}{2}\lambda_0\left(\ln J\right)^2-\mu_0\ln J+\frac{1}{2} \mu_0\left(\mathrm{tr}(\boldsymbol{C})-3\right) where :math:`\lambda_0` and :math:`\mu_0` are, respectively, Lamé’s first parameter and the shear modulus at the initial configuration. :math:`J` is the jacobian of the deformation gradient (:math:`\boldsymbol{F}=\nabla_{\!\!\boldsymbol{X}}\boldsymbol{x}`): :math:`J=\text{det}(\boldsymbol{F})`. Finally :math:`\boldsymbol{C}` is the right Cauchy-Green deformation tensor. Since this kind of material is used for large deformation problems, a finite deformation framework should be used. Therefore, the Cauchy stress (:math:`\boldsymbol{\sigma}`) should be computed through the second Piola-Kirchhoff stress tensor :math:`\boldsymbol{S}`: .. math:: \boldsymbol{\sigma } = \frac{1}{J}\boldsymbol{F}\boldsymbol{S}\boldsymbol{F}^T Finally the second Piola-Kirchhoff stress tensor is given by: .. math:: \boldsymbol{S} = 2\frac{\partial\Psi}{\partial\boldsymbol{C}} = \lambda_0\ln J \boldsymbol{C}^{-1}+\mu_0\left(\boldsymbol{I}-\boldsymbol{C}^{-1}\right) The parameters to indicate in the material file are the same as those for the elastic case: ``E`` (Young’s modulus), ``nu`` (Poisson’s ratio). ----- .. _sect-smm-cl-ve: Visco-Elastic ````````````` .. _sect-smm-cl-sls: Standard-Linear Solid ''''''''''''''''''''' Keyword: **sls_deviatoric** Inherits from **elastic** Material description with input file: .. code-block:: python #input.dat material sls_deviatoric [ name = material_name rho = 1000 # density (Real) E = 2.1e9 # young's modulus (Real) nu = 0.4 # poisson's ratio (Real) Eta = 1. # Viscosity (Real) Ev = 0.5 # Stiffness of viscous element (Real) Plane_stress = false # Plane stress simplification (bool, only 2D problems) ] Since this material inherits from :cpp:class:`MaterialElastic ` the parameter :math:`E` preceeds :math:`E_{\mathrm{inf}}`. So only :math:`E` and :math:`E_v` can be set, :math:`E_{\mathrm{inf}}` is deduced by :math:`E_{\mathrm{inf}} = E - E_{v}`. ----- Visco-elasticity is characterized by strain rate dependent behavior. Moreover, when such a material undergoes a deformation it dissipates energy. This dissipation results in a hysteresis loop in the stress-strain curve at every loading cycle (see :numref:`fig:smm:cl:visco-elastic:hyst`). In principle, it can be applied to many materials, since all materials exhibit a visco-elastic behavior if subjected to particular conditions (such as high temperatures). .. figure:: figures/cl/stress_strain_visco.svg :name: fig:smm:cl:visco-elastic:hyst :align: center :width: 40.0% Characteristic stress-strain behavior of a visco-elastic material with hysteresis loop .. figure:: figures/cl/visco_elastic_law.svg :name: fig:smm:cl:visco-elastic:model :align: center :width: 40.0% Schematic representation of the standard rheological linear solid visco-elastic model The standard rheological linear solid model (see Sections 10.2 and 10.3 of :cite:`simo92`) has been implemented in ``Akantu``. This model results from the combination of a spring mounted in parallel with a spring and a dashpot connected in series, as illustrated in :numref:`fig:smm:cl:visco-elastic:model`. The advantage of this model is that it allows to account for creep or stress relaxation. The equation that relates the stress to the strain is (in 1D): .. math:: \frac{d\varepsilon(t)}{dt} = \left ( E_{\mathrm{inf}} + E_{v} \right ) ^ {-1} \cdot \left [ \frac{d\sigma(t)}{dt} + \frac{E_{v}}{\eta}\sigma(t) - \frac{E_{\mathrm{inf}}E_V}{\eta}\varepsilon(t) \right ] where :math:`\eta` is the viscosity. The equilibrium condition is unique and is attained in the limit, as :math:`t \to \infty`. At this stage, the response is elastic and depends on the Young’s modulus :math:`E`. The mandatory parameters for the material file are the following: ``rho`` (density), ``E`` (Young’s modulus), ``nu`` (Poisson’s ratio), ``Plane_Stress`` (if set to zero plane strain, otherwise plane stress), ``eta`` (dashpot viscosity) and ``Ev`` (stiffness of the viscous element). Note that the current standard linear solid model is applied only on the deviatoric part of the strain tensor. The spheric part of the strain tensor affects the stress tensor like an linear elastic material. ----- .. _sect-smm-cl-maxwell: Maxwell Chain Visco-Elasticity '''''''''''''''''''''''''''''' Keyword: **viscoelastic_maxwell** Inherits from **elastic** Material description with input file: .. code-block:: python #input.dat material viscoelastic_maxwell [ name = material_name rho = 1000 # density (Real) Einf = 5.e9 # Infinite time Young's modulus (Real) nu = 0.4 # poisson's ratio (Real) Ev = [1.e9, 2.e9, 3.e9] # Maxwell elements' stiffness values (Vector) Eta = [1.e14, 2.e16, 3.e16] # Dashpot elements' viscosity values (Vector) Plane_stress = false # Plane stress simplification (bool, only 2D problems) ] ---- .. figure:: figures/cl/maxwell_chain.png :name: fig-smm-cl-visco-elastic-maxwell :align: center :width: 40.0% Schematic representation of the Maxwell chain A different visco-elastic rheological model available to users is the generalized Maxwell chain (see :cite:`de_borst_finiteelement_1994` and Section 46.7.4 of :cite:`diana_manual`). It consists of a series of sequential spring-dashpots (Maxwell elements) placed in parallel with one single spring (see :numref:`fig-smm-cl-visco-elastic-maxwell`). The relation between stresses and strain comes from .. math:: \sigma \left ( t \right ) = \int_{-\infty}^{t} E \left ( t, \tau \right ) \mathbf{D} \dot{\varepsilon} d\tau where :math:`E(t,\tau)` is the time-dependent relaxation function, :math:`\tau` is the loading age, and :math:`\mathbf{D}` is the dimensionless matrix relating a 3D deformation state to a 1D relaxation function. The relaxation function is expanded in the exponential series .. math:: :label: eqn-relaxation-function E \left ( t, \tau \right ) = E_{0} + \sum_{\alpha=1}^{n} E_{\alpha} e^{- \frac{t- \tau}{\lambda_{\alpha}}} where the relaxation time of each Maxwell element is defined as :math:`\lambda_{\alpha}=\eta_{\alpha} / E_{\alpha}` with :math:`\eta_{\alpha}` being the viscosity of a dash-pot. Assuming a constant strain rate within each time step, the analytical integration of the right-hand side of :eq:`eqn-relaxation-function` leads to the following form .. math:: \sigma \left ( t + \Delta t \right ) = E_{0} \mathbf{D} \varepsilon + \sum_{\alpha=1}^n \left ( \left ( 1 - e^{\frac{- \Delta t}{\lambda_{\alpha}}} \right ) \frac{E_{\alpha} \lambda_{\alpha}}{\Delta t} \mathbf{D} \delta \varepsilon + e^{\frac{-\Delta t}{\lambda_{\alpha}}} \sigma_{\alpha} \left ( t \right ) \right ) with :math:`\sigma_{\alpha}(t)` being the internal stress within each Maxwell element, defined as .. math:: \sigma_{\alpha} \left ( t \right ) = \mathbf{D} \int_0^t E_{\alpha} e^{\frac{-t- \tau}{\lambda_{\alpha}}} \dot{\varepsilon} d \tau The first term under the sum sign in above equation could be seen as the effective stiffness of a single Maxwell element multiplied by the matrix :math:`\mathbf{D}` and the strain increment :math:`\Delta \varepsilon`: .. math:: E_{\alpha}^{ef} = \left ( 1- e^{\frac{-Δt}{λ_α}} \right ) \frac{E_α λ_α}{Δt} Time increment :math:`Δt` controls the rate dependency of the effective stiffness. By limit analysis, we find the limiting values of the effective stiffness which are equal to :math:`E_0` for infinitely slow loading (:math:`Δt` tending to infinity) and :math:`E_0+ΣE_α` for infinitely fast (:math:`Δt` tending to 0). At the end of each converged time step, the internal stress :math:`σ_α(t)` is updated according to .. math:: σ_α \left ( t \right ) = σ_α \left ( t - Δt \right ) e^{\frac{-Δt}{λ_α}} + E_α^{ef} \mathbf{D} Δε The mandatory parameters for the material file are the following: ``rho`` (density), ``nu`` (Poisson’s ratio), ``Plane_Stress`` (if set to zero plane strain, otherwise plane stress), ``Einf`` (infinite time Young’s modulus), ``Ev`` (Maxwell elements' stiffness values stored in a vector), ``Eta`` (dashpots' viscosity values stored in a vector). The Maxwell model is applied on the entire strain tensor and does not distinguish between its deviatoric and hydrostatic components. Note that the time step has to be specified for the model using current material both for static and dynamic simulations: .. code-block:: c++ model.setTimeStep(time_step_value); .. _sect-smm-cl-plastic: Plastic ``````` Small-Deformation Plasticity '''''''''''''''''''''''''''' Keyword: **plastic_linear_isotropic_hardening** Inherits from **elastic** Material description with input file: .. code-block:: python #input.dat material plastic_linear_isotropic_hardening [ name = material_name rho = 1000 # density (Real) E = 2.1e9 # young's modulus (Real) nu = 0.4 # poisson's ratio (Real) h = 0.1 # Hardening modulus (Real) sigma_y = 1e6 # Yield stress (Real) ] Energies: - ``potential``: elastic part of the potential energy - ``plastic``: dissipated plastic energy (integrated over time) ----- The small-deformation plasticity is a simple plasticity material formulation which accounts for the additive decomposition of strain into elastic and plastic strain components. This formulation is applicable to infinitesimal deformation where the additive decomposition of the strain is a valid approximation. In this formulation, plastic strain is a shearing process where hydrostatic stress has no contribution to plasticity and consequently plasticity does not lead to volume change. :numref:`fig:smm:cl:Lin-strain-hard` shows the linear strain hardening elasto-plastic behavior according to the additive decomposition of strain into the elastic and plastic parts in infinitesimal deformation as .. math:: \boldsymbol{\varepsilon} &= \boldsymbol{\varepsilon}^e +\boldsymbol{\varepsilon}^p\\ \boldsymbol{\sigma} &= 2G(\boldsymbol{\varepsilon}^e) + \lambda \mathrm{tr}(\boldsymbol{\varepsilon}^e)\boldsymbol{I} .. figure:: figures/cl/isotropic_hardening_plasticity.svg :name: fig:smm:cl:Lin-strain-hard :align: center Stress-strain curve for the small-deformation plasticity with linear isotropic hardening. In this class, the von Mises yield criterion is used. In the von Mises yield criterion, the yield is independent of the hydrostatic stress. Other yielding criteria such as Tresca and Gurson can be easily implemented in this class as well. In the von Mises yield criterion, the hydrostatic stresses have no effect on the plasticity and consequently the yielding occurs when a critical elastic shear energy is achieved. .. math:: \label{eqn:smm:constitutive:von Mises} f = \sigma_{{\mathrm{eff}}} - \sigma_y = \left(\frac{3}{2} {\boldsymbol{\sigma}}^{{\mathrm{tr}}} : {\boldsymbol{\sigma}}^{{\mathrm{tr}}}\right)^\frac{1}{2}-\sigma_y (\boldsymbol{\varepsilon}^p) .. math:: \label{eqn:smm:constitutive:yielding} f < 0 \quad \textrm{Elastic deformation,} \qquad f = 0 \quad \textrm{Plastic deformation} where :math:`\sigma_y` is the yield strength of the material which can be function of plastic strain in case of hardening type of materials and :math:`{\boldsymbol{\sigma}}^{{\mathrm{tr}}}` is the deviatoric part of stress given by .. math:: \label{eqn:smm:constitutive:deviatoric stress} {\boldsymbol{\sigma}}^{{\mathrm{tr}}}=\boldsymbol{\sigma} - \frac{1}{3} \mathrm{tr}(\boldsymbol{\sigma}) \boldsymbol{I} After yielding :math:`(f = 0)`, the normality hypothesis of plasticity determines the direction of plastic flow which is normal to the tangent to the yielding surface at the load point. Then, the tensorial form of the plastic constitutive equation using the von Mises yielding criterion (see equation 4.34) may be written as .. math:: \label{eqn:smm:constitutive:plastic contitutive equation} \Delta {\boldsymbol{\varepsilon}}^p = \Delta p \frac {\partial{f}}{\partial{\boldsymbol{\sigma}}}=\frac{3}{2} \Delta p \frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}} In these expressions, the direction of the plastic strain increment (or equivalently, plastic strain rate) is given by :math:`\frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}}` while the magnitude is defined by the plastic multiplier :math:`\Delta p`. This can be obtained using the *consistency condition* which impose the requirement for the load point to remain on the yielding surface in the plastic regime. Here, we summarize the implementation procedures for the small-deformation plasticity with linear isotropic hardening: #. Compute the trial stress: .. math:: {\boldsymbol{\sigma}}^{{\mathrm{tr}}} = {\boldsymbol{\sigma}}_t + 2G\Delta \boldsymbol{\varepsilon} + \lambda \mathrm{tr}(\Delta \boldsymbol{\varepsilon})\boldsymbol{I} #. Check the Yielding criteria: .. math:: f = (\frac{3}{2} {\boldsymbol{\sigma}}^{{\mathrm{tr}}} : {\boldsymbol{\sigma}}^{{\mathrm{tr}}})^{1/2}-\sigma_y (\boldsymbol{\varepsilon}^p) #. Compute the Plastic multiplier: .. math:: \begin{aligned} d \Delta p &= \frac{\sigma^{tr}_{eff} - 3G \Delta P^{(k)}- \sigma_y^{(k)}}{3G + h}\\ \Delta p^{(k+1)} &= \Delta p^{(k)}+ d\Delta p\\ \sigma_y^{(k+1)} &= (\sigma_y)_t+ h\Delta p \end{aligned} #. Compute the plastic strain increment: .. math:: \Delta {\boldsymbol{\varepsilon}}^p = \frac{3}{2} \Delta p \frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}} #. Compute the stress increment: .. math:: {\Delta \boldsymbol{\sigma}} = 2G(\Delta \boldsymbol{\varepsilon}-\Delta \boldsymbol{\varepsilon}^p) + \lambda \mathrm{tr}(\Delta \boldsymbol{\varepsilon}-\Delta \boldsymbol{\varepsilon}^p)\boldsymbol{I} #. Update the variables: .. math:: \begin{aligned} {\boldsymbol{\varepsilon^p}} &= {\boldsymbol{\varepsilon}}^p_t+{\Delta {\boldsymbol{\varepsilon}}^p}\\ {\boldsymbol{\sigma}} &= {\boldsymbol{\sigma}}_t+{\Delta \boldsymbol{\sigma}} \end{aligned} We use an implicit integration technique called *the radial return method* to obtain the plastic multiplier. This method has the advantage of being unconditionally stable, however, the accuracy remains dependent on the step size. The plastic parameters to indicate in the material file are: :math:`\sigma_y` (Yield stress) and ``h`` (Hardening modulus). In addition, the elastic parameters need to be defined as previously mentioned: ``E`` (Young’s modulus), ``nu`` (Poisson’s ratio). ----- Damage `````` In the simplified case of a linear elastic and brittle material, isotropic damage can be represented by a scalar variable :math:`d`, which varies from :math:`0` to :math:`1` for no damage to fully broken material respectively. The stress-strain relationship then becomes: .. math:: \boldsymbol{\sigma} = (1-d)\, \boldsymbol{C}:\boldsymbol{\varepsilon} where :math:`\boldsymbol{\sigma}`, :math:`\boldsymbol{\varepsilon}` are the Cauchy stress and strain tensors, and :math:`\boldsymbol{C}` is the elastic stiffness tensor. This formulation relies on the definition of an evolution law for the damage variable. In ``Akantu``, many possibilities exist and they are listed below. ---- .. _sect-smm-cl-damage-marigo: Marigo '''''' Keyword: **marigo** Inherits from **elastic** Material description with input file: .. code-block:: python #input.dat material marigo [ name = material_name rho = 1000 # density (Real) E = 2.1e9 # young's modulus (Real) nu = 0.4 # poisson's ratio (Real) Plane_stress = false # Plane stress simplification (bool, only 2D problems) Yd = 0.1 # Hardening modulus (Random) Sd = 1. # Damage energy (Real) ] Energies: - ``dissipated``: energy dissipated in damage ----- This damage evolution law is energy based as defined by Marigo :cite:`marigo81a`, :cite:`lemaitre96a`. It is an isotropic damage law. .. math:: \begin{aligned} Y &= \frac{1}{2}\boldsymbol{\varepsilon}:\boldsymbol{C}:\boldsymbol{\varepsilon}\\ F &= Y - Y_d - S d\\ d &= \left\{ \begin{array}{l l} \mathrm{min}\left(\frac{Y-Y_d}{S},\;1\right) & \mathrm{if}\; F > 0\\ \mathrm{unchanged} & \mathrm{otherwise} \end{array} \right.\end{aligned} In this formulation, :math:`Y` is the strain energy release rate, :math:`Y_d` the rupture criterion and :math:`S` the damage energy. The non-local version of this damage evolution law is constructed by averaging the energy :math:`Y`. .. _sect-smm-cl-damage-mazars: Mazars '''''' Keyword: **mazars** Inherits from **elastic** Material description with input file: .. code-block:: python #input.dat material mazars [ name = concrete rho = 3000 # density (Real) E = 32e9 # young's modulus (Real) nu = 0.2 # poisson's ratio (Real) K0 = 9.375e-5 # Damage threshold (Real) At = 1.15 # Parameter damage traction 1 (Real) Bt = 10000 # Parameter damage traction 2 (Real) Ac = 0.8 # Parameter damage compression 1 (Real) Bc = 1391.3 # Parameter damage compression 2 (Real) beta = 1.00 # Parameter for shear (Real) ] Energies: - ``dissipated``: energy dissipated in damage ----- This law introduced by Mazars :cite:`mazars84a` is a behavioral model to represent damage evolution in concrete. This model does not rely on the computation of the tangent stiffness, the damage is directly evaluated from the strain. The governing variable in this damage law is the equivalent strain :math:`\varepsilon_{{\mathrm{eq}}} = \sqrt{<\boldsymbol{\varepsilon}>_+:<\boldsymbol{\varepsilon}>_+}`, with :math:`<.>_+` the positive part of the tensor. This part is defined in the principal coordinates (I, II, III) as :math:`\varepsilon_{{\mathrm{eq}}} = \sqrt{<\boldsymbol{\varepsilon_I}>_+^2 + <\boldsymbol{\varepsilon_{II}}>_+^2 + <\boldsymbol{\varepsilon_{III}}>_+^2}`. The damage is defined as: .. math:: \begin{aligned} D &= \alpha_t^\beta D_t + (1-\alpha_t)^\beta D_c\\ D_t &= 1 - \frac{\kappa_0 (1- A_t)}{\varepsilon_{{\mathrm{eq}}}} - A_t \exp^{-B_t(\varepsilon_{{\mathrm{eq}}}-\kappa_0)}\\ D_c &= 1 - \frac{\kappa_0 (1- A_c)}{\varepsilon_{{\mathrm{eq}}}} - A_c \exp^{-B_c(\varepsilon_{{\mathrm{eq}}}-\kappa_0)}\\ \alpha_t &= \frac{\sum_{i=1}^3<\varepsilon_i>_+\varepsilon_{{\mathrm{nd}}\;i}}{\varepsilon_{{\mathrm{eq}}}^2}\end{aligned} With :math:`\kappa_0` the damage threshold, :math:`A_t` and :math:`B_t` the damage parameter in traction, :math:`A_c` and :math:`B_c` the damage parameter in compression, :math:`\beta` is the shear parameter. :math:`\alpha_t` is the coupling parameter between traction and compression, the :math:`\varepsilon_i` are the eigenstrain and the :math:`\varepsilon_{{\mathrm{nd}}\;i}` are the eigenvalues of the strain if the material were undamaged. The coefficients :math:`A` and :math:`B` are the post-peak asymptotic value and the decay shape parameters. .. _sect:smm:CLNL: Non-Local Constitutive Laws ``````````````````````````` Continuum damage modeling of quasi-brittle materials undergo significant softening after the onset of damage. This fast growth of damage causes a loss of ellipticity of partial differential equations of equilibrium. Therefore, the numerical simulation results won't be objective anymore, because the dissipated energy will depend on mesh size used in the simulation. One way to avoid this effect is the use of non-local damage formulations. In this approach a local quantity such as the strain is replaced by its non-local average, where the size of the domain, over which the quantitiy is averaged, depends on the underlying material microstructure. ``Akantu`` provides non-local versions of many constitutive laws for damage. Examples are for instance the material :ref:`sect-smm-cl-damage-mazars` and the material :ref:`sect-smm-cl-damage-marigo`, that can be used in a non-local context. In order to use the corresponding non-local formulation the user has to define the non-local material he wishes to use in the text input file: .. code-block:: none material constitutive_law_non_local [ name = material_name rho = $value$ ... ] where ``constitutive_law_non_local`` is the name of the non-local constitutive law, *e.g.* `marigo_non_local`. In addition to the material the non-local neighborhood, that should be used for the averaging process needs to be defined in the material file as well: .. code-block:: none non_local neighborhood_name weight_function_type [ radius = $value$ ... weight_function weight_parameter [ damage_limit = $value$ ... ] ] for the non-local averaging, *e.g.* ``base_wf``, followed by the properties of the non-local neighborhood, such as the radius, and the weight function parameters. It is important to notice that the non-local neighborhood must have the same name as the material to which the neighborhood belongs! The following two sections list the non-local constitutive laws and different type of weight functions available in ``Akantu``. \subsection{Non-local constitutive laws} Let us consider a body having a volume :math:`V` and a boundary :math:`\Gamma`. The stress-strain relation for a non-local damage model can be described as follows: .. _eq:non-local-const: .. math:: \vec{\sigma} = (1-\bar{d}) \vec{D}:\epsilon with :math:`\vec{D}` the elastic moduli tensor, :math:`\sigma` the stress tensor, :math:`\epsilon` the strain tensor and :math:`\bar{d}` the non-local damage variable. Note that this stres-strain relationship is similar to the relationship defined in Damage model except :math:`\bar{d}`. The non-local damage model can be extended to the damage constitutive laws: :ref:`sect-smm-cl-damage-marigo` and :ref:`sect-smm-cl-damage-mazars`. The non-local damage variable :math:`\bar{d}` is defined as follows: .. _eq:non-local-damage: .. math:: \bar{d}(\vec{x}) = \int_{V}W(\vec{x}, \vec{y}) d(\vec{y}) dV(\vec{y}) with :math:`W(\vec{x},\vec{y})` the weight function which averages local damage variables to describe the non-local interactions. A list of available weight functions and its functionalities in \akantu are explained in the next section. Non-local weight functions '''''''''''''''''''''''''' The available weight functions in ``Akantu`` are follows: - ``base_weight_function``: This weight function averages local damage variables by using a bell-shape function on spatial dimensions. - ``damaged_weight_function``: A linear-shape weight function is applied to average local damage variables. Its slope is determined by damage variables. For example, the damage variables for an element which is highly damaged are averaged over large spatial dimension (linear function including a small slope). - ``remove_damaged_weight_function``: This weight function averages damage values by using a bell-shape function as ``base_weight_function``, but excludes elements which are fully damaged. - ``remove_damaged_with_damage_rate_weight_function``: A bell-shape function is applied to average local damage variables for elements having small damage rates. - ``stress_based_weight_function``: Non local integral takes stress states, and use the states to construct weight function: an ellipsoid shape. Detailed explanations of this weight function are given in Giry et al. :cite:`giry13a`. ----- .. _sec-cohesive-laws: Cohesive Constitutive laws `````````````````````````` .. _ssect-smm-cl-coh-snozzi: Linear Irreversible Law ''''''''''''''''''''''' Keyword: **cohesive_linear** Material description with input file: .. code-block:: python #input.dat material cohesive_linear [ name = cohesive sigma_c = 0.1 # critical stress sigma_c (default: 0) G_c = 1e-2 # Mode I fracture energy beta = 0. # weighting parameter for sliding and normal opening (default: 0) penalty = 0. # stiffness in compression to prevent penetration (α in the text) kappa = 1. # ration between mode-I and mode-II fracture energy (Gc_II/Gc_I) contact_after_breaking = true # Activation of contact when the elements are fully damaged max_quad_stress_insertion = false # Insertion of cohesive element when stress is high # enough just on one quadrature point # if false the average stress on facet's quadrature points is used ] ----- .. figure:: figures/cl/linear_cohesive_law.svg :alt: Irreversible cohesive laws for explicit simulations. :name: fig:smm:coh:linear_cohesive_law :align: center :width: 60.0% Irreversible cohesive laws for explicit simulations. `Akantu` includes the Snozzi-Molinari :cite:`snozzi_cohesive_2013` linear irreversible cohesive law (see :numref:`fig:smm:coh:linear_cohesive_law`). It is an extension to the Camacho-Ortiz :cite:`camacho_computational_1996` cohesive law in order to make dissipated fracture energy path-dependent. The concept of free potential energy is dropped and a new independent parameter :math:`\kappa` is introduced: .. math:: \kappa = \frac{G_\mathrm{c, II}}{G_\mathrm{c, I}} where :math:`G_\mathrm{c, I}` and :math:`G_\mathrm{c, II}` are the necessary works of separation per unit area to open completely a cohesive zone under mode I and mode II, respectively. Their model yields to the following equation for cohesive tractions :math:`\vec{T}` in case of crack opening :math:`{\delta}`: .. math:: \vec{T} = \left( \frac{\beta^2}{\kappa} \Delta_\mathrm{t} \vec{t} + \Delta_\mathrm{n} \vec{n} \right) \frac{\sigma_\mathrm{c}}{\delta} \left( 1- \frac{\delta}{\delta_\mathrm{c}} \right) = \hat{\vec T}\, \frac{\sigma_\mathrm{c}}{\delta} \left( 1- \frac{\delta}{\delta_\mathrm{c}} \right) :label: eq-smm-coh-tractions where :math:`\sigma_\mathrm{c}` is the material strength along the fracture, :math:`\delta_\mathrm{c}` the critical effective displacement after which cohesive tractions are zero (complete decohesion), :math:`\Delta_\mathrm{t}` and :math:`\Delta_\mathrm{n}` are the tangential and normal components of the opening displacement vector :math:`\vec{\Delta}`, respectively. The parameter :math:`\beta` is a weight that indicates how big the tangential opening contribution is. The effective opening displacement is: .. math:: \delta = \sqrt{\frac{\beta^2}{\kappa^2} \Delta_\mathrm{t}^2 + \Delta_\mathrm{n}^2} In case of unloading or reloading :math:`\delta < \delta_\mathrm{max}`, tractions are calculated as: .. math:: \begin{eqnarray} T_\mathrm{n} &= \Delta_\mathrm{n}\, \frac{\sigma_\mathrm{c}}{\delta_\mathrm{max}} \left( 1- \frac{\delta_\mathrm{max}}{\delta_\mathrm{c}} \right) \\ T_\mathrm{t} &= \frac{\beta^2}{\kappa}\, \Delta_\mathrm{t}\, \frac{\sigma_\mathrm{c}}{\delta_\mathrm{max}} \left( 1- \frac{\delta_\mathrm{max}}{\delta_\mathrm{c}} \right) \end{eqnarray} so that they vary linearly between the origin and the maximum attained tractions. As shown in :numref:`fig:smm:coh:linear_cohesive_law`, in this law, the dissipated and reversible energies are: .. math:: \begin{eqnarray} E_\mathrm{diss} &= \frac{1}{2} \sigma_\mathrm{c}\, \delta_\mathrm{max}\\[1ex] E_\mathrm{rev} &= \frac{1}{2} T\, \delta \end{eqnarray} Moreover, a damage parameter :math:`D` can be defined as: .. math:: D = \min \left( \frac{\delta_\mathrm{max}}{\delta_\mathrm{c}},1 \right) which varies from 0 (undamaged condition) and 1 (fully damaged condition). This variable can only increase because damage is an irreversible process. A simple penalty contact model has been incorporated in the cohesive law so that normal tractions can be returned in case of compression: .. math:: T_\mathrm{n} = \alpha \Delta_\mathrm{n} \quad\text{if}\quad \Delta_\mathrm{n}\quad <\quad 0 where :math:`\alpha` is a stiffness parameter that defaults to zero. The relative contact energy is equivalent to reversible energy but in compression. The material name of the linear decreasing cohesive law is ``material_cohesive_linear`` and its parameters with their respective default values are: - ``sigma_c = 0`` - ``delta_c = 0`` - ``beta = 0`` - ``G_c = 0`` - ``kappa = 1`` - ``penalty = 0`` where ``G_c`` corresponds to :math:`G_\mathrm{c, I}`. A random number generator can be used to assign a random :math:`\sigma_\mathrm{c}` to each facet following a given distribution (see Section :ref:`sect-smm-cl`). Only one parameter between ``delta_c`` and ``G_c`` has to be specified. For random :math:`\sigma_\mathrm{c}` distributions, the chosen parameter of these two is kept fixed and the other one is varied. The bi-linear constitutive law works exactly the same way as the linear one, except for the additional parameter ``delta_0`` that by default is zero. Two examples for the extrinsic and intrinsic cohesive elements and also an example to assign different properties to inter-granular and trans-granular cohesive elements can be found in the folder ``examples/cohesive_element/``. ---- .. _ssect:smm:cl:coh-friction: Linear Cohesive Law with Friction ''''''''''''''''''''''''''''''''' Keyword: **cohesive_linear_friction** Material description with input file: .. code-block:: python #input.dat material cohesive_linear_friction [ name = interface beta = 1 # weighting parameter for sliding and normal opening (default: 0) G_c = 30e-3 # Mode I fracture energy penalty = 1.0e6 # stiffness in compression to prevent penetration (α in the text) sigma_c = 2.0 # critical stress sigma_c (default: 0) contact_after_breaking = true # Activation of contact when the elements are fully damaged mu = 0.5 # Maximum value of the friction coefficient penalty_for_friction = 5.0e3 # Penalty parameter for the friction behavior ] ----- This law represents a variation of the linear irreversible cohesive of the previous section, which adds friction. The friction behavior is approximated with an elasto-plastic law, which relates the friction force to the relative sliding between the two faces of the cohesive element. The slope of the elastic branch is called ``penalty_for_friction``, and is defined by the user, together with the friction coefficient, as a material property. The friction contribution evolves with the damage of the cohesive law: it is null when the damage is zero, and it becomes maximum when the damage is equal to one. This is done by defining a current value of the friction coefficient (mu) that increases linearly with the damage, up to the value of the friction coefficient defined by the user. The yielding plateau of the friction law is given by the product of the current friction coefficient and the local compression stress acting in the cohesive element. Such an approach is equivalent to a node-to-node contact friction. Its accuracy is acceptable only for small displacements. The material name of the linear cohesive law with friction is ``material_cohesive_linear_friction``. Its additional parameters with respect to those of the linear cohesive law without friction, with the respective default values, are: - ``mu = 0`` - ``penalty_for_friction = 0`` .. _ssect:smm:cl:coh-fatigue: Linear Cohesive Law with Fatigue '''''''''''''''''''''''''''''''' Keyword: **cohesive_linear_fatigue** Material description with input file: .. code-block:: python #input.dat material cohesive_linear_fatigue [ name = cohesive sigma_c = 1 # critical stress sigma_c (default: 0) beta = 1 # weighting parameter for sliding and normal opening (default: 0) delta_c = 1 # Critical displacement delta_f = 1 # delta_f (normalization of opening rate to alter reloading stiffness after fatigue) count_switches = true # Count the opening/closing switches per element ] ----- This law represents a variation of the linear irreversible cohesive law of the previous section, that removes the hypothesis of elastic unloading-reloading cycles. With this law, some energy is dissipated also during unloading and reloading with hysteresis. The implementation follows the work of :cite:`nguyen2001`. During the unloading-reloading cycle, the traction increment is computed as .. math:: \dot{T} = \begin{cases} K^- \, \dot{\delta} & \text{if $\dot{\delta} < 0$} \\ K^+ \, \dot{\delta} & \text{if $\dot{\delta} > 0$} \\ \end{cases} where :math:`\dot{\delta}` and :math:`\dot{T}` are respectively the effective opening displacement and the cohesive traction increments with respect to time, while :math:`K^-` and :math:`K^+` are respectively the unloading and reloading incremental stiffness. The unloading path is linear and results in an unloading stiffness .. math:: K^- = \frac{T_\mathrm{max}}{\delta_\mathrm{max}} where :math:`T_\mathrm{max}` and :math:`\delta_\mathrm{max}` are the maximum cohesive traction and the effective opening displacement reached during the precedent loading phase. The unloading stiffness remains constant during the unloading phase. On the other hand the reloading stiffness increment :math:`\dot{K}^+` is calculated as .. math:: \dot{K}^+ = \begin{cases} - K^+ \, \dot{\delta} / \delta_\mathrm{f} & \text{if $\dot{\delta} > 0$} \\ \left( K^+ - K^- \right) \, \dot{\delta} / \delta_\mathrm{f} & \text{if $\dot{\delta}$ < $0$} \end{cases} where :math:`\delta_\mathrm{f}` is a material parameter (refer to :cite:`vocialta15` for more details). During unloading the stiffness :math:`K^+` tends to :math:`K^-`, while during reloading :math:`K^+` gets decreased at every time step. If the cohesive traction during reloading exceeds the upper limit given by equation :eq:`eq-smm-coh-tractions`, it is recomputed following the behavior of the linear decreasing cohesive law for crack opening. .. _ssect:smm:cl:coh-exponential: Exponential Cohesive Law ''''''''''''''''''''''''' Keyword: **cohesive_exponential** Material description with input file: .. code-block:: python #input.dat material cohesive_exponential [ name = coh1 sigma_c = 1.5e6 # critical stress sigma_c (default: 0) beta = 1 # weighting parameter for sliding and normal opening (default: 0) delta_c = 1e-4 # Critical displacement exponential_penalty = true # Is contact penalty following the exponential law? contact_tangent = 1.0 # Ratio of contact tangent over the initial exponential tangent ] ----- Ortiz and Pandolfi proposed this cohesive law in 1999 :cite:`ortiz1999`. The traction-opening equation for this law is as follows: .. math:: T = e \sigma_c \frac{\delta}{\delta_c}e^{-\delta/ \delta_c} :label: eq:exponential_law This equation is plotted in :numref:`fig:smm:cl:ecl`. The term :math:`\partial{\vec{T}}/ \partial{\delta}` after the necessary derivation can expressed as .. math:: \frac{\partial{\vec{T}}} {\partial{\delta}} = \hat{\vec{T}} \otimes \frac {\partial{(T/\delta)}}{\partial{\delta}} \frac{\hat{\vec{T}}}{\delta}+ \frac{T}{\delta} \left[ \beta^2 \mat{I} + \left(1-\beta^2\right) \left(\vec{n} \otimes \vec{n}\right)\right] :label: eq:tangent_cohesive where .. math:: \frac{\partial{(T/ \delta)}}{\partial{\delta}} = \left\{\begin{array} {l l} -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad \text{if} \delta \geq \delta_{max}\\ 0 & \quad \text{if} \delta < \delta_{max}, \delta_n > 0 \end{array} \right. As regards the behavior in compression, two options are available: a contact penalty approach with stiffness following the formulation of the exponential law and a contact penalty approach with constant stiffness. In the second case, the stiffness is defined as a function of the tangent of the exponential law at the origin. .. figure:: figures/cl/cohesive_exponential.png :alt: Exponential cohesive law :name: fig:smm:cl:ecl :align: center Exponential cohesive law diff --git a/doc/dev-doc/manual/contactmechanicsmodel.rst b/doc/dev-doc/manual/contactmechanicsmodel.rst new file mode 100644 index 000000000..617895736 --- /dev/null +++ b/doc/dev-doc/manual/contactmechanicsmodel.rst @@ -0,0 +1,472 @@ +Contact Mechanics Model +======================= + +The contact mechanics model is a specific implementation of +:cpp:class:`Model ` interface to handle contact between +bodies. + +Theory +------ + +.. _fig-contactmechanicsmodel-schematic: + +.. figure:: figures/contact_mechanics_schematic.png + :align: center + :width: 90% + + Basic notation for the contact between two bodies. + +Let us consider two deformable bodies, represented as :math:`\Omega_\alpha`, +:math:`\alpha=1, 2`. The boundary :math:`\Gamma_\alpha` of a body is divided +into three non-intersecting regions : :math:`\Gamma^D_\alpha` with +prescribed displacements, :math:`\Gamma^N_\alpha` with prescribed tractions +and :math:`\Gamma^C_\alpha` where the two bodies :math:`\Omega_1` can potentially +come into contact such that: + +.. math:: + \Gamma^D_\alpha \cup \Gamma^N_\alpha \cup \Gamma^C_\alpha = + \Gamma_\alpha, \quad \Gamma^D_\alpha \cap \Gamma^N_\alpha \cap + \Gamma^C_\alpha = \emptyset + + +The motion of the two bodies is described in a fixed spatial frame +defined by orthonormal basis :math:`[\boldsymbol{e}_x, +\boldsymbol{e}_y, \boldsymbol{e}_z]` by mapping +:math:`\mathcal{M}_\alpha^t` in time interval :math:`t \in [0, T]`. In +*reference* configuration :math:`i.e.~ t=0`, the position vector for +an arbitrary point on body :math:`\Omega_\alpha` is represented as +:math:`\boldsymbol{X}_\alpha` and in actual configuration the points +are denoted by small letters for example :math:`\boldsymbol{x}`. +:numref:`fig-contactmechanicsmodel-schematic` shows the motion of two +bodies from the *reference* to actual configuration. During the +motion, the bodies can potentially come in contact along +:math:`\Gamma^C_\alpha` as shown in +:numref:`fig-contactmechanicsmodel-schematic`, termed as *potential +contact zone*. Upon contact, two physical conditions need to be +satisfied :math:`(i)` the two surfaces :math:`\Gamma^C_\alpha` cannot +interpenetrate at any time during motion :math:`(ii)` forces must be +exerted by bodies along the :math:`\Gamma^C_\alpha` to resist the +interpenetration as well as relative motion along the contacting +surfaces. Of the two conditions, the non-penetration of bodies can be +defined by a *gap function* :math:`g` which represents the separation +between the two bodies: + +.. math:: + g = (\boldsymbol{r} -\boldsymbol{\rho}).\boldsymbol{n} + + +where :math:`\boldsymbol{r} \in \mathcal{M}[\Gamma^C_1]` is the +position of point :math:`\boldsymbol{X}_1` at time :math:`t` , +:math:`\boldsymbol{\rho} \in \mathcal{M}[\Gamma^C_2]` is the closest +point projection of :math:`\boldsymbol{r}` and :math:`\boldsymbol{n}` +is the outward normal at :math:`\boldsymbol{\rho}` (see +:numref:`fig-contactmechanicsmodel-schematic` b). To preclude the +interpenetration, the *gap function* is constrained with :math:`g \geq 0`, as it is shown in +\Cref{fig:body-contact}b. When the two bodies eventually come in +contact, the gap vanishes :math:`i.e.~ g=0` which leads to the +development of tractions :math:`\boldsymbol{T}_\alpha` along the +contact interface. Thus, the boundary value problem can be formulated +for the two bodies :math:`\Omega_\alpha` with contact constraints as an +extra boundary condition: + +.. math:: + \boldsymbol{T}\Big |_{{\Gamma^C}_{\alpha}} = \boldsymbol{T}_\alpha + +Solution spaces :math:`\mathcal{U}_\alpha` and weighting spaces +:math:`\mathcal{W}_\alpha` are defined for each body: + +.. math:: + \mathcal{U}_\alpha = \{ \boldsymbol{u}_\alpha: \Omega_\alpha \to \mathbb{R}^d ~|~ \boldsymbol{u}=\boldsymbol{u}[\boldsymbol{x}_\alpha]~ \forall~\boldsymbol{x}_\alpha \in \Gamma^D_\alpha \} + + +.. math:: + \mathcal{W}_\alpha = \{ \boldsymbol{w}_\alpha: \Omega_\alpha \to + \mathbb{R}^d ~|~ \boldsymbol{w}=0~ \forall~ \boldsymbol{x}_\alpha \in + \Gamma^D_\alpha \} + + +The variational form for each body :math:`\alpha` is given as + +.. math:: + \int_{\Omega_\alpha}\boldsymbol{\sigma}[\boldsymbol{u}_\alpha]:\boldsymbol{\epsilon}[\boldsymbol{w}_\alpha]~d\Omega_\alpha + = + \int_{\Omega_\alpha}\boldsymbol{b}_\alpha.\boldsymbol{w}_\alpha~d\Omega_\alpha + + \int_{\Gamma^C_\alpha}\boldsymbol{T}_\alpha.\boldsymbol{w}_\alpha~d\Gamma_\alpha + + \int_{\Gamma^N_\alpha}\boldsymbol{T}_\alpha^D.\boldsymbol{w}_\alpha~d\Gamma^N_\alpha + + +In equilibrium state, following Newton's :math:`3^{rd}` law, it can be +stated the forces exerted by two bodies along :math:`\Gamma^C_\alpha` +are equal and opposite: + +.. math:: + \boldsymbol{T}_1~d\Gamma^C_1 = - \boldsymbol{T}_2~d\Gamma^C_2 + + +This allows to replace the two integrals over contact surfaces +$\Gamma^C_\alpha$ with a single integral over any one of the surfaces: + +.. math:: + \int_\alpha + \boldsymbol{T}_\alpha.\delta\boldsymbol{u}_\alpha~d\Gamma^C_\alpha &= + \int_{\Gamma^C_1}\boldsymbol{T}_1\delta\boldsymbol{u}_1~d\Gamma^C_1 + + \int_{\Gamma^C_2} \boldsymbol{T}_2\delta\boldsymbol{u}_2~d\Gamma^C_2 \\ &= + \int_{\Gamma^C_1} \boldsymbol{T}_1(\delta\boldsymbol{u}_1 - + \delta\boldsymbol{u}_2)~d\Gamma^C_1 + + +The contact traction :math:`\boldsymbol{T}_1` can be decomposed into +its normal and tangential components as: + +.. math:: + \boldsymbol{T}_\alpha = \boldsymbol{T}^n + \boldsymbol{T}^t + = \sigma_n\boldsymbol{n}+ \boldsymbol{T}^t + +where :math:`\boldsymbol{T}^n` is the component along the normal +vector :math:`\boldsymbol{n}` and :math:`\boldsymbol{T}^t` is the +component tangential to :math:`\boldsymbol{n}` developed due to the +friction along the surfaces. Upon contact , the normal pressure +:math:`\sigma_n` associated to :math:`\boldsymbol{T}^n` must be +compressive to preclude the interpenetration of bodies. In general, if +a point is not in contact :math:`g > 0`, then :math:`\sigma_n=0` and +in contact :math:`\sigma_n < 0`. This leads to the non-penetration +condition: + +.. math:: + \sigma_n g = 0 + + +The above set of conditions care called *Hertz-Signorini-Moreau* or +*Karush-Kuhn-Tucker* condition given as: + +.. math:: + g \geq 0, \quad \sigma_n \leq 0 \quad g.\sigma_n = 0 + + +The tangential component of the contact traction is defined as : + +.. math:: + \boldsymbol{T}_t = (\boldsymbol{I} - \boldsymbol{n} \otimes \boldsymbol{n})\boldsymbol{\sigma} + + +where :math:`\boldsymbol{\sigma}` is the Cauchy stress tensor and +$\boldsymbol{n}$ is the outward normal at +:math:`\boldsymbol{\rho}`. The direction of tangential traction, +:math:`\boldsymbol{s}` is determined by the relative sliding velocity +:math:`\boldsymbol{v}_t` of the point :math:`\boldsymbol{r}` and its +projection point :math:`\boldsymbol{\rho}` in contact and is given as: + +.. math:: + \boldsymbol{s} = \begin{cases} \dfrac{\boldsymbol{v}_t}{\| + \boldsymbol{v}_t \|}, & \text{ if } \| \boldsymbol{v}_t\| > 0 \\ & + \\ 0, & \text{ if } \| \boldsymbol{v}_t\| = 0 \end{cases} + + +According to the experimental observations of Amontons +and~\cite{coulomb}, in presence of friction, the interface develops a +frictional strength which governs the sliding between them and thus +constrains the tangential contact. In Akantu implementation, we +restrict to the classical non-associated Coulomb's friction +law~\citep{coulomb} which is widely used in many physical and +engineering applications. The Coulomb's friction law, as a first-order +approximation, states that the frictional strength is proportional to +the normal pressure: + +.. math:: + \sigma_{fric} = \mu |\sigma_n| + + +where :math:`\mu` is coefficient of friction between interfaces. If +the tangential traction :math:`||\boldsymbol{T}^t||` developed is +below the frictional strength, the relative tangential sliding is zero +:math:`i.e.~\boldsymbol{v}_t = 0` : + +.. math:: + ||\boldsymbol{T}^t|| < \mu| \sigma_n |, \quad \boldsymbol{v}_t=0 + + +The above equation denotes a *stick state*. As soon as the tangential +stress reaches the frictional strength, the two surfaces start +slipping relative to each other :math:`i.e.~\boldsymbol{v}_t > 0`. The +slipping of the surfaces ensures that the tangential stress does not +exceeds the frictional strength, :math:`||\boldsymbol{\sigma}_t|| +-\mu|\sigma_n| = 0`. This definition of *slip state* is defined as: + +.. math:: + ||\boldsymbol{\sigma}_t|| -\mu|\sigma_n| = 0, \quad + ||\boldsymbol{v}_t|| > 0 + + +Similar to *Karush-Kuhn-Tucker* condition for normal contact, the +above conditions formulate the necessary conditions for tangential +contact: + +.. math:: + ||\boldsymbol{v}_t|| \geq 0, \quad ||\boldsymbol{T}^t|| + -\mu|\sigma_n| = 0, \quad ||\boldsymbol{v}_t|| + \Big(||\boldsymbol{T}^t|| -\mu|\sigma_n| \Big) = 0 + + +The above contact and frictional constraints are unilateral in nature +*i.e.* they do not behave symmetrically with respect to *gap function, +g*. This renders the balance of work as a variational inequality: + +.. math:: + \sum_{\alpha=1}^2 + \int_{\Omega_\alpha}\boldsymbol{\sigma}[\boldsymbol{u}_\alpha]:\boldsymbol{\epsilon}[\boldsymbol{w}_\alpha]~d\Omega_\alpha + \geq \sum_{\alpha=1}^2 \Big \{ + \int_{\Omega_\alpha}\boldsymbol{b}_\alpha.\boldsymbol{w}_\alpha~d\Omega_\alpha + + \int_{\Gamma^Ct_\alpha}\boldsymbol{T}_\alpha.\boldsymbol{w}_\alpha~d\Gamma_\alpha + + \int_{\Gamma^N_\alpha}\boldsymbol{T}_\alpha^D.\boldsymbol{w}_\alpha~d\Gamma_\alpha\Big + \} + + +which makes it a non-linear optimization problem. The strategy +employed (optimization techniques) to find the solution of variational +inequality depends on the choice of numerical framework employed to +solve the B.V.P. + +To solve the minimization problem, FEM introduces the concept of +active set strategy to overcome the problem. In active set strategy, +it is assumed that at current solution step, the part of potential +contact zone :math:`\Gamma^C_1` that are in contact are known, +:math:`\Gamma^{C\star}_1 \subseteq \Gamma^C_1`. This +is achieved by first allowing the two bodies to interpenetrate and +finding the interpenetrated part of potential contact zone which will +denote the active set. Knowing the active set, transforms the +optimization problem to a variational equality where constraints are +imposed along the active part of contacting interface: + +.. math:: + + \sum_{\alpha=1}^2 + \int_{\Omega_\alpha}\boldsymbol{\sigma}[\boldsymbol{u}_\alpha]:\boldsymbol{\epsilon}[\delta\boldsymbol{u}_\alpha]~d\Omega_\alpha + = \sum_{\alpha=1}^2 \Big \{ + \int_{\Omega_\alpha}\boldsymbol{b}_\alpha.\delta\boldsymbol{u}_\alpha~d\Omega_\alpha + + \int_{\Gamma^N_\alpha}\boldsymbol{T}_\alpha^D.\delta\boldsymbol{u}_\alpha~d\Gamma_\alpha\Big + \} \\ + + \int_{\Gamma^{C\star}_1 + }\boldsymbol{T}_1(\delta\boldsymbol{u}_1 - + \delta\boldsymbol{u}_2)~d\Gamma^{C\star}_1 + + +Thus, the resolution of contact problem in FEM requires two steps: +finding the active set along the contacting interface and then +imposing the contact constraints along the active set only. In the +following section, we describe how to employ the contact detection strategies +implemented within Akantu in order to find the active set and to compute contact forces. + + +Using the Contact Mechanics Model +--------------------------------- + +The :cpp:class:`ContactMechanicsModel ` +object solves the contact problem. An instance of +the class can be created like this:: + + ContactMechanicsModel contact(mesh, spatial_dimension); + +while an existing mesh has been used (see \ref{sect:common:mesh}). To +intialize the model object:: + + contact.initFull(_analysis_method = _explicit_lumped_mass); + +The contact mechanics model contains :cpp:class:`Arrays `: + +:cpp:func:`gaps ` + contains the nodal interpenetrating value :math:`g` (positive for + interpenetration, zero by default after initialization) + +:cpp:func:`normals ` + contains the normal vector at the slave nodes (zero by default + after initialization). + + +In Akantu, the possible contact between surfaces is divided into 3 +categories. + +- Physical Surfaces - The contact occurs between two pyhsically + defined boundaries/surfaces of a body. +- Cohesive Surfaces - The contact occurs between fracturing surfaces + created using :cpp:class:`SolidMechanicsModelCohesive + `. +- All Surfaces - The contact can occur between physical as well as + cohesive surfaces. + + +To select the contacting surfaces, one must define a +:cpp:class:`SurfaceSelector ` of one of the +above defined types. + +To define contact between Physical surfaces, an instance of +:cpp:class:`PhysicalSurfaceSelector ` +is created where the mesh object (see \ref{sect:common:mesh}) is +passed as an argument:: + + auto && surface_selector = std::make_shared(mesh); + + +To define contact between cohesive surfaces, an instance of +:cpp:class:`CohesiveSurfaceSelector ` +must be created. As the contact occurs between the cohesive facets, +therefore the mesh facet object is passed as an argument:: + + auto && surface_selector = std::make_shared(mesh.getMeshFacets()); + + +To defind contact between physical and cohesive surfaces, an instance of +:cpp:class:`AllSurfaceSelector ` +must be created. As the contact occurs between the cohesive facets, +therefore the mesh facet object is passed as an argument:: + + auto && surface_selector = std::make_shared(mesh.getMeshFacets()); + + +Once a surface selector is created it must be assigned to the +:cpp:class:`ContactDetector ` class:: + + contact.getContactDetector().setSurfaceSelector(surface_selector); + + +Contact detection +''''''''''''''''' + +The contact detection algorithm can receive the a few parameters. It is possible to +specify the master/slave surfaces with their string identifier. The geometrical projections +are performed with iterations which can be controlled as a classical optimization problem. +A typical detection configuration is given below: + +.. code-block:: + + contact_detector [ + type = explicit + master = contact_bottom + slave = contact_top + projection_tolerance = 1e-10 + max_iterations = 100 + extension_tolerance = 1e-5 + ] + + + + +Contact resolution +'''''''''''''''''' + +The contact resolution defines the surface tractions due to contact +both for normal and tangential contact. A typical configuration for a +penalization formulation is as follows: + +.. code-block:: + + contact_resolution penalty_linear [ + name = contact_top + mu = 0.0 # friction coefficient + epsilon_n = 4e5 # normal penalization + epsilon_t = 1e5 # friction penalization + is_master_deformable = false # when master is rigid => computational savings + ] + + +Coupling with :cpp:class:`SolidMechanicsModel ` +'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +To couple the +:cpp:class:`ContactMechancisModel` +contact mechanics model with +:cpp:class:`SolidMechanicsModel` a +dedicated coupler class :cpp:class:`CouplerSolidContact` is provided. + +When an instance of a coupler class is created, it automatically +creates the instances of solid mechanics model and contact mechanics +model. The two objects can be retrived from the coupler class. + +.. code-block:: c++ + + CouplerSolidContact coupler(mesh); + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + +Simply initializing the coupler initializes the two models. + +.. code-block:: c++ + + coupler.initFull( _analysis_method = _explicit_lumped_mass); + +However two set the material selector and the contact detector for the +two models, one must set them using directly the instance of the two +model classes. + +.. code-block:: c++ + + auto && selector = std::make_shared>( + "physical_names",solid); + solid.setMaterialSelector(selector); + + +.. code-block:: c++ + + auto && surface_selector = std::make_shared(mesh); + contact.getContactDetector().setSurfaceSelector(surface_selector); + +The dumping fields/vectors belonging to the solid mechanics model and +contact mechanics model can directly be set through the coupler +class. + +.. code-block:: c++ + + coupler.setBaseName("contact-explicit-dynamic"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normal_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("areas"); + coupler.addDumpField("stress"); + + +Finally to solve the two models :cpp:func:`solveStep +` function of coupler class must be +invoked. + +.. code-block:: c++ + + coupler.solveStep(); + + +Coupling with :cpp:class:`SolidMechanicsModelCohesive ` +''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +To use the contact mechanics model with cohesive elements, one must use the +:cpp:class:`CouplerSolidCohesiveContact` to +coupler the :cpp:class:`ContactMechancisModel` +with +:cpp:class:`SolidMechanicsModelCohesive`. +The initialization and invocation of the functions are similar to +:cpp:class:`CouplerSolidContact` except a few +changes. + +.. code-block:: c++ + + solid = coupler.getSolidMechanicsModelCohesive(); + + +While initializing the coupler, the nature of cohesive elements +(extrinsic/intrinsic) should need to be passed. + +.. code-block:: c++ + + coupler.initFull( _analysis_method = _explicit_lumped_mass, _is_extrinsic=true); + + +To ensure that cohesive elements break during an explicit insertion, one must +call the function :cpp:func:`checkCohesiveStress() +` after +:cpp:func:`solveStep() `. + +.. code-block:: c++ + + coupler.solveStep(); + solid.checkCohesiveStress(); + diff --git a/doc/dev-doc/manual/figures/contact_mechanics_schematic.png b/doc/dev-doc/manual/figures/contact_mechanics_schematic.png new file mode 100644 index 000000000..f69ab9b71 Binary files /dev/null and b/doc/dev-doc/manual/figures/contact_mechanics_schematic.png differ diff --git a/doc/dev-doc/manual/figures/contact_mechanics_schematic.svg b/doc/dev-doc/manual/figures/contact_mechanics_schematic.svg new file mode 100644 index 000000000..7d3d0e27e --- /dev/null +++ b/doc/dev-doc/manual/figures/contact_mechanics_schematic.svg @@ -0,0 +1,1195 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + { + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/dev-doc/manual/phasefieldmodel.rst b/doc/dev-doc/manual/phasefieldmodel.rst index 2cbf7b15a..1935beb16 100644 --- a/doc/dev-doc/manual/phasefieldmodel.rst +++ b/doc/dev-doc/manual/phasefieldmodel.rst @@ -1,265 +1,265 @@ PhaseField Model ================ The phasefield model is a specific implementation of :cpp:class:`Model ` interface to handle brittle fracture for infinitesimal strains. Theory ------ The variational formulation of brittle fracture was first proposed by Francfort and Marigo:cite:`franc` so as to overcome the shortcomings of Griffith's criteria. The approach used by Francfort and Margio for variational formulation is based on the principle of global minimality of total energy. Similar to Griffith's work, they define a surface energy corresponding to the discontinuity :math:`\Gamma` as .. math:: \Phi_{s}(\Gamma) = \int\limits_{\Gamma}\mathcal{G}_cd\mathcal{H}^{d-1}(\Gamma) where :math:`\mathcal{H}^{d-1}` is d-1 dimensional Hausdroff measure which is surface measure for smooth hypersurfaces. the material behavior is considered to be linearly elastic with small strains throughout the body. Based on these assumptions the elastic energy of the body is defined as .. math:: \Phi_e(\Gamma,u) = \int\limits_{\Omega_{p}}\psi(\boldsymbol{\epsilon}(u))d\Omega where :math:`\psi` is elastic energy density and it is given as a function of strain as .. math:: \psi_e(\boldsymbol{\epsilon}) = \dfrac{1}{2}\lambda(\text{tr}(\boldsymbol{\epsilon}))^{2}+\mu\boldsymbol{\epsilon}:\boldsymbol{\epsilon} Exponential phasefield law '''''''''''''''''''''''''' Bourdin et. al.:cite:`bourdin` in their work proposed a regularized version of variational formulation where a scalar field variable :math:`d(\vec{x},t)` is used to represent crack. This scalar field variable approximates the sharp crack topology by taking value 1 at crack location and smoothly diffusing into value 0 away from the crack. \figref{fig:diffusive_topology} shows both sharp crack topology and an approximated regularized crack topology for one-dimensional case. The strong form of the phasefield can be expressed as .. math:: 2(1-d)\mathcal{H} - \dfrac{G_c}{l_0}(d-l_0^2d) = 0 \nabla d . \boldsymbol n = 0 Using the PhaseField Model -------------------------- The :cpp:class:`PhaseFieldModel ` object solves the strong form using an inplicit solver. An instance of the class can be created like this:: PhaseFieldModel phase(mesh, spatial_dimension); while ans existing mesh has been used (see \ref{sect:common:mesh}). To intialize the model object:: phase.initFull(); Currently, implicit solver is defined for the phasefield model and no explicit solver is implemented in Akantu. Furthermore, By default, the implicit solver defined is of linear type. One can change the linear solver to a non-linear by initiating a new solver type. The phasefield model contains :cpp:class:`Arrays `: :cpp:func:`damage `: contains the nodal damage :math:`d` (zero by default after the initialization) :cpp:func:`blocked_dofs ` contains a Boolean value specifying whether the damage at a node is to be blocked or not. A Dirichlet boundary condition can be prescribed by setting the **blocked_dofs** value of damage to ``true``. The **damage** ais computed for all nodes where the **blocked_dofs** value is set to ``false``. For the remaining nodes, the imposed values (zero by default after initialization) are kept. :cpp:func:`internal_force ` contains the driving force responsible for the crack to nucelate or propagate. Currently, the phasefield model uses a exponential shaped scalar field variable :math:`d(\vec{x}, t)` to approximate the sharp crack topology. A Phasefield variable thus requires a length scale parameter :math:`l_0` to control the width of the diffusive crack topology along with the material parameters such as elastic modulus, poisson's ratio, critical energy release rate. The data input file provides the parameters for the exponential phasefield law as follows: .. code-block:: phasefield exponential [ name = plate E = 210.0 nu = 0.3 gc = 5e-3 l0 = 0.1 ] :cpp:class:`PhaseFieldModel` can handle phasefield laws for multiple materials. To define so: .. code-block:: phasefield exponential [ name = hard E = 210.0 nu = 0.3 gc = 5e-3 l0 = 0.1 ] phasefield exponential [ name = soft E = 21.0 nu = 0.3 gc = 5e-5 l0 = 0.1 ] In order to assign correct phasefield variable properties based on the names of region as defined in mesh file, :cpp:class:`MeshDataPhaseFieldSelector ` must be set as phasefield selector:: auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); For the crack to nucleate or propagate, phasefield model requires strain measure at each quadrature point. The strain measure computed from solid mechanics model is provided to the phasefield model. Similarly, damage value computed at each quadrature point by the phasefield model is provided to the solid mechanics model. This damage is thus used to degrade the the elastic strain energy. To do so, a new material which is a specific implementation of :cpp:class:`MaterialDamage ` is defined for solid mechanics model. The governing equation is given as .. math:: \boldsymbol \sigma = ((1-d)^2 + \eta)\dfrac{\psi}{\boldsymbol \epsilon} where :math:`\psi` is the elastic energy and :math:`\eta` is a numerical parmaeter avoid numerical difficulties due to full degradation of elastic energy for fully broken state. The material properties are thus provided in the input data filled as: .. code-block:: material phasefield [ name = hard rho = 1. E = 210.0 nu = 0.3 eta = 0.0 Plane_Stress = false ] material phasefield [ name = soft rho = 1. E = 21.0 nu = 0.3 eta = 0.0 Plane_Stress = false ] To simplify the execution of phasefield model coupled with solidmechanis model, a special class :cpp:class:`CouplerSolidPhaseField` is provided. Coupling Phase Field Model and Solid Mechanics Model '''''''''''''''''''''''''''''''''''''''''''''''''''' A dedicated coupler class :cpp:class:`CouplerSolidPhaseField ` is defined in Akantu to ease the coupling of the two models. When an instance of a Coupler class is created, it automatically creates the instances of solid mechanics model and phasefield model. The two objects can be retrived from the coupler class. .. code-block:: c++ CouplerSolidPhaseField coupler(mesh); auto & phase = coupler.getPhaseFieldModel(); auto & solid = coupler.getSolidMechanicsModel(); The two objects must be used to define the solver type and apply boundary conditions. .. code-block:: c++ solid.initFull(_analysis_method = _explicit_lumped_mass); phase.initFull(_analysis_method = _static); -The whole proces of coupling the two models at a given time step is +The whole process of coupling the two models at a given time step is made easy by the :cpp:func:`solve ` function of coupler class. .. code-block:: c++ coupler.solve(, ); Staggered scheme ---------------- To solve the solid mechanics model and phasefield model, staggered scheme is implemented. In staggered scheme, at current time step first solid mechanics model is solved assuming the damage values from the previous time step. The strain thus computed are passed to the phasefield model and now, phasefield model is solved for the damage variable. At each iteration step, convergence in displacement and damage is checked for. If convergence is not reached, the combined newton-raphson iteration continues. To illustrate the staggered scheme solution of phasefield model, one can refere to the examples provided for the phasefield model. Two examples provided are for static `phasefield-static.py` and dynamic crack propgation `phasefield-dynamic.py`. For the static problem, both the solid mechanics model and the phsefield model are solved using alinear implicit solver. Convergence in value of both displacement as well as damage is checked at each loading step. In case of the dynamic problem, solid mechanics model is solved dynamically using an explicit solver and the phasefield model is solved using a linear implcit solver. In this sceanrio, the staggered scheme doesnot check for any convergence. Below is the crack propagation observed for the dyanmic problem. .. _fig-phasefieldmodel-dynamic: .. figure:: figures/phasefield-dynamic.png :align: center Dynamic crack propagation using phasefield model .. LocalWords: phasefield SolidMechanics PhaseField akantu cpp diff --git a/doc/manual/manual-contactmechanicsmodel.tex b/doc/manual/manual-contactmechanicsmodel.tex new file mode 100644 index 000000000..fec7da8b7 --- /dev/null +++ b/doc/manual/manual-contactmechanicsmodel.tex @@ -0,0 +1,9 @@ +\chapter{Contact Mechanics Model\index{ContactMechanicsModel}\label{sect:cmm}} + +The contact mechanics model is a specific implementation of the +\code{Model} interface dedicated to handle the contact between two deformable bodies. + +\section{The contact virtual work} +\begin{align} +G_c(\phi, \phi^{*}) = \int_{\Gamma_c^{(1)}}\Big[t_{N}\delta g + t_{T_{\alpha}}\delta\xi^{\alpha}\Big]d\Gamma +\end{align} diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index bcdace931..5a44de139 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,55 +1,55 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Fri Oct 22 2010 # @date last modification: Fri Jan 22 2016 # # @brief List of examples # # @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 . # # @section DESCRIPTION # #=============================================================================== add_example(new_material "Example on how to add a new material in Akantu" PACKAGE core) add_example(boundary_conditions "Example on hoy to set boundary conditions" PACKAGE core) add_example(explicit "Example on how to run an explicit simulation" PACKAGE core) add_example(io "Example on how to perform Input/Output operations" PACKAGE core) add_example(implicit "Example on how to run an implicit simulation" PACKAGE implicit) add_example(static "Example on how to run a static simulation" PACKAGE implicit) add_example(parallel "Example of how to write a parallel code with Akantu" PACKAGE parallel) add_example(cohesive_element "Cohesive element examples" PACKAGE cohesive_element) add_example(structural_mechanics "Structural mechanics model examples" PACKAGE structural_mechanics) add_example(heat_transfer "Example on how to run heat transfer simulation" PACKAGE heat_transfer) add_example(python "Example on how to use the python interface" PACKAGE python_interface) add_example(embedded "Example on how to run embedded model simulation" PACKAGE embedded) +add_example(contact_mechanics "Example on how to run contact mechanics model simulation" PACKAGE contact_mechanics) add_example(phase_field "Example on how to run phase field model simulation" PACKAGE phase_field) - package_add_files_to_package( examples/README.rst cmake/AkantuExampleMacros.cmake ) diff --git a/packages/model_couplers.cmake b/examples/contact_mechanics/CMakeLists.txt similarity index 60% copy from packages/model_couplers.cmake copy to examples/contact_mechanics/CMakeLists.txt index 47c375b8c..ddcad7d0f 100644 --- a/packages/model_couplers.cmake +++ b/examples/contact_mechanics/CMakeLists.txt @@ -1,47 +1,53 @@ #=============================================================================== -# @file model_couplers.cmake +# @file CMakeLists.txt # -# @author Mohit Pundir +# @author Mohit Pundir # -# @date creation: Sun Sep 30 2018 -# @date last modification: Sun Sep 28 2018 +# @date creation: Mon Jan 18 2016 # -# @brief package description for model couplers +# @brief configuration for heat transfer example # # @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 . # +# @section DESCRIPTION +# #=============================================================================== +add_mesh(hertz hertz.geo 2 1) +add_mesh(cohesive-contact cohesive-contact.geo 2 1) -package_declare(model_couplers ADVANCED - DESCRIPTION "Use Model Couplers package of Akantu") -package_declare_sources(model_couplers - model/model_couplers/coupler_solid_phasefield.hh - model/model_couplers/coupler_solid_phasefield.cc +register_example(contact_explicit_dynamic + SOURCES contact_explicit_dynamic.cc + DEPENDS hertz + FILES_TO_COPY material.dat ) -package_declare_documentation_files(model_couplers - # +register_example(contact_explicit_static + SOURCES contact_explicit_static.cc + DEPENDS hertz + FILES_TO_COPY material.dat ) -package_declare_documentation(model_couplers - "This package activates the modle couplers within Akantu. " - "It has no additional dependencies." + +register_example(cohesive_contact_explicit_dynamic + SOURCES cohesive_contact_explicit_dynamic.cc + DEPENDS cohesive-contact + FILES_TO_COPY material-cohesive.dat ) diff --git a/examples/contact_mechanics/cohesive-contact.geo b/examples/contact_mechanics/cohesive-contact.geo new file mode 100755 index 000000000..101577c68 --- /dev/null +++ b/examples/contact_mechanics/cohesive-contact.geo @@ -0,0 +1,36 @@ +h = 0.5; + +y = 2; +x = 1; + +Point(1) = { x/2., y/2, 0, h}; +Point(2) = {-x/2., y/2, 0, h}; +Point(3) = {-x/2.,-y/2, 0, h}; +Point(4) = { x/2.,-y/2, 0, h}; +Point(5) = {-x/2., 0, 0, h}; +Point(6) = { x/2., 0, 0, h}; + +Line(1) = {1, 2}; +Line(2) = {2, 5}; +Line(3) = {5, 6}; +Line(4) = {6, 1}; + +Line(5) = {5, 3}; +Line(6) = {3, 4}; +Line(7) = {4, 6}; + +Line Loop(1) = {1, 2, 3, 4}; +Line Loop(2) = {-3, 5, 6, 7}; +Plane Surface(1) = {1}; +Plane Surface(2) = {2}; + +Physical Line("fixed") = {6}; +Physical Line("loading") = {1}; +Physical Line("interface") = {3}; +Physical Line("sides") = {2, 5, 7, 4}; + +Physical Surface("body") = {1, 2}; + +Recombine Surface "*"; +Transfinite Surface "*"; +Mesh.SecondOrderIncomplete = 1; diff --git a/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc b/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc new file mode 100644 index 000000000..9a6a82843 --- /dev/null +++ b/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc @@ -0,0 +1,128 @@ +/* -------------------------------------------------------------------------- */ +#include "contact_mechanics_model.hh" +#include "coupler_solid_cohesive_contact.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "surface_selector.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char * argv[]) { + + const UInt spatial_dimension = 2; + initialize("material-cohesive.dat", argc, argv); + + Real time_step{0.}; + Real time_factor = 0.1; + UInt max_steps = 25000; + Real max_displacement = 1e-3; + + Mesh mesh(spatial_dimension); + mesh.read("cohesive-contact.msh"); + + CouplerSolidCohesiveContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModelCohesive(); + auto & contact = coupler.getContactMechanicsModel(); + + auto && material_selector = + std::make_shared(solid); + material_selector->setFallback(solid.getMaterialSelector()); + solid.setMaterialSelector(material_selector); + + auto && surface_selector = std::make_shared(mesh); + contact.getContactDetector().setSurfaceSelector(surface_selector); + + coupler.initFull(_analysis_method = _explicit_lumped_mass, + _is_extrinsic = true); + + coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "sides"); + + time_step = solid.getStableTimeStep(); + time_step *= time_factor; + std::cout << "Time Step = " << time_step << "s (" << time_step << "s)" + << std::endl; + coupler.setTimeStep(time_step); + + coupler.setBaseName("cohesive-contact-explicit-dynamic"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("velocity"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + coupler.addDumpField("gaps"); + coupler.addDumpField("areas"); + + auto & velocity = solid.getVelocity(); + auto & gaps = contact.getGaps(); + + Real damping_ratio = 0.99; + auto increment = max_displacement / max_steps; + + for (auto i : arange(max_steps)) { + + coupler.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "loading"); + coupler.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "fixed"); + + coupler.solveStep(); + + solid.checkCohesiveStress(); + + // damping velocities only along the contacting zone + for (auto && tuple : zip(gaps, make_view(velocity, spatial_dimension))) { + auto & gap = std::get<0>(tuple); + auto & vel = std::get<1>(tuple); + if (gap > 0) { + vel *= damping_ratio; + } + } + + // dumping energies + if (i % 1000 == 0) { + + Real epot = solid.getEnergy("potential"); + Real ekin = solid.getEnergy("kinetic"); + + std::cerr << i << "," << i * increment << "," << epot << "," << ekin + << "," << epot + ekin << "," << std::endl; + } + + if (i % 1000 == 0) { + coupler.dump(); + } + } + + for (auto i : arange(max_steps)) { + + solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "loading"); + solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "fixed"); + + coupler.solveStep(); + + coupler.checkCohesiveStress(); + + // damping velocities only along the contacting zone + for (auto && tuple : zip(gaps, make_view(velocity, spatial_dimension))) { + auto & gap = std::get<0>(tuple); + auto & vel = std::get<1>(tuple); + if (gap > 0) { + vel *= damping_ratio; + } + } + + // dumping energies + if (i % 1000 == 0) { + + Real epot = solid.getEnergy("potential"); + Real ekin = solid.getEnergy("kinetic"); + + std::cerr << i << "," << i * increment << "," << epot << "," << ekin + << "," << epot + ekin << "," << std::endl; + } + + if (i % 1000 == 0) { + coupler.dump(); + } + } +} diff --git a/examples/contact_mechanics/contact_explicit_dynamic.cc b/examples/contact_mechanics/contact_explicit_dynamic.cc new file mode 100644 index 000000000..fd7956f4c --- /dev/null +++ b/examples/contact_mechanics/contact_explicit_dynamic.cc @@ -0,0 +1,104 @@ +/* -------------------------------------------------------------------------- */ +#include "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +#include "surface_selector.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material.dat", argc, argv); + + Real time_step; + Real time_factor = 0.1; + UInt max_steps = 20000; + Real max_displacement = 5e-3; + + Mesh mesh(spatial_dimension); + mesh.read("hertz.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + auto && selector = std::make_shared>( + "physical_names",solid); + solid.setMaterialSelector(selector); + + coupler.initFull( _analysis_method = _explicit_lumped_mass); + + auto && surface_selector = std::make_shared(mesh); + contact.getContactDetector().setSurfaceSelector(surface_selector); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "fixed"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "fixed"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "loading"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "symmetry"); + + time_step = solid.getStableTimeStep(); + time_step *= time_factor; + std::cout << "Time Step = " << time_step << "s (" << time_step + << "s)" << std::endl; + coupler.setTimeStep(time_step); + + coupler.setBaseName("contact-explicit-dynamic"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("velocity"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("areas"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + auto & velocity = solid.getVelocity(); + auto & gaps = contact.getGaps(); + + Real damping_ratio = 0.99; + auto increment = max_displacement/max_steps; + + for (auto i : arange(max_steps)) { + + solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "loading"); + + coupler.solveStep(); + + // damping velocities only along the contacting zone + for(auto && tuple : zip(gaps, + make_view(velocity, spatial_dimension))){ + auto & gap = std::get<0>(tuple); + auto & vel = std::get<1>(tuple); + if(gap > 0) { + vel *= damping_ratio; + } + } + + // dumping energies + if (i % 1000 == 0) { + + Real epot = solid.getEnergy("potential"); + Real ekin = solid.getEnergy("kinetic"); + + std::cerr << i << "," << i * increment << "," << epot << "," << ekin << "," + << epot + ekin << "," << std::endl; + } + + if (i % 1000 == 0) { + coupler.dump(); + } + } + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/examples/contact_mechanics/contact_explicit_static.cc b/examples/contact_mechanics/contact_explicit_static.cc new file mode 100644 index 000000000..2df0bff15 --- /dev/null +++ b/examples/contact_mechanics/contact_explicit_static.cc @@ -0,0 +1,96 @@ +/** + * @file test_contact_mechanics_model.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model class + * + * @section LICENSE + * + * 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +#include "surface_selector.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material.dat", argc, argv); + + Mesh mesh(spatial_dimension); + mesh.read("hertz.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + auto && selector = std::make_shared>( + "physical_names",solid); + solid.setMaterialSelector(selector); + + coupler.initFull(_analysis_method = _static); + + auto && surface_selector = std::make_shared(mesh); + contact.getContactDetector().setSurfaceSelector(surface_selector); + + coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "fixed"); + coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "fixed"); + coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "loading"); + coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "symmetry"); + + coupler.setBaseName("contact-explicit-static"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("areas"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + auto max_steps = 100u; + + for (auto _ [[gnu::unused]] : arange(max_steps)) { + + auto increment = 1e-4; + coupler.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "loading"); + + coupler.solveStep(); + coupler.dump(); + } + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/examples/contact_mechanics/hertz.geo b/examples/contact_mechanics/hertz.geo new file mode 100644 index 000000000..c788ceed8 --- /dev/null +++ b/examples/contact_mechanics/hertz.geo @@ -0,0 +1,47 @@ +cl1 = 0.0075; +cl2 = 0.005; +cl3 = 0.005; +Dy = 0.0; +radius = 0.1; +y = 0.1; +epsilon = 1e-3; + +nb_points = 51; + +Point(1) = {0, y, 0, 0.25*cl1}; +Point(2) = {radius, radius + y, 0, cl1}; + +Point(4) = {0.2, y, 0, 0.5*cl1}; +Point(5) = {0, y, 0, 0.5*cl1}; +Point(6) = {0, 0.025, 0, 10*cl3}; +Point(7) = {0.2, 0.025, 0, 10*cl3}; +Point(8) = {0, radius + y, 0, cl1}; + +Circle(2) = {1, 8, 2}; + +Line(3) = {2, 8}; +Line(13) = {8, 1}; +Line(4) = {6, 7}; +Line(5) = {7, 4}; +Line(6) = {4, 5}; +Line(7) = {5, 6}; + +Line Loop(9) = {2, 3, 13}; +Plane Surface(9) = {9}; + +Line Loop(11) = {6, 7, 4, 5}; +Plane Surface(11) = {11}; + +Physical Line("contact_bottom") = {6}; +Physical Line("contact_top") = {2}; + +Physical Line("loading") = {3}; +Physical Line("fixed") = {4}; + +Physical Surface("upper") = {9}; +Physical Surface("lower") = {11}; + +Physical Line("sides") = {5}; + +Physical Line("symmetry") = {13, 7}; +Recombine Surface "*"; diff --git a/examples/contact_mechanics/material-cohesive.dat b/examples/contact_mechanics/material-cohesive.dat new file mode 100644 index 000000000..073cfb9dd --- /dev/null +++ b/examples/contact_mechanics/material-cohesive.dat @@ -0,0 +1,40 @@ +model solid_mechanics_model_cohesive[ +cohesive_inserter [ + cohesive_surfaces = [interface] + ] + +material elastic [ + name = body + rho = 1.0 # density + E = 1e3 # young's modulu + nu = 0.3 # poisson's ratio +] + +material cohesive_linear [ + name = interface + beta = 1 + G_c = 6.25e-4 + sigma_c = 1 + kappa = 1 + penalty = 0 + contact_after_breaking = false +] + +] + +contact_detector [ + type = explicit + master = interface + slave = interface + projection_tolerance = 1e-10 + max_iterations = 100 + extension_tolerance = 1e-5 +] + +contact_resolution penalty_linear [ + name = interface + mu = 0.0 + epsilon_n = 0 + epsilon_t = 0 + is_master_deformable = true +] diff --git a/examples/contact_mechanics/material.dat b/examples/contact_mechanics/material.dat new file mode 100644 index 000000000..9a299facc --- /dev/null +++ b/examples/contact_mechanics/material.dat @@ -0,0 +1,30 @@ +material elastic [ + name = upper + rho = 1.0 # density + E = 1e3 # young's modulu + nu = 0.3 # poisson's ratio +] + +material elastic [ + name = lower + rho = 1.0 # density + E = 1e3 # young's modulu + nu = 0.3 # poisson's ratio +] + +contact_detector [ + type = explicit + master = contact_bottom + slave = contact_top + projection_tolerance = 1e-10 + max_iterations = 100 + extension_tolerance = 1e-5 +] + +contact_resolution penalty_linear [ + name = contact_top + mu = 0.0 + epsilon_n = 4e5 + epsilon_t = 1e5 + is_master_deformable = false +] diff --git a/examples/python/CMakeLists.txt b/examples/python/CMakeLists.txt index 260733f51..161c19dae 100644 --- a/examples/python/CMakeLists.txt +++ b/examples/python/CMakeLists.txt @@ -1,11 +1,12 @@ add_subdirectory(custom-material) add_subdirectory(dynamics) add_subdirectory(eigen_modes) add_subdirectory(plate-hole) add_subdirectory(stiffness_matrix) add_subdirectory(structural_mechanics) add_subdirectory(fragmentation) +add_subdirectory(contact-mechanics) package_add_files_to_package( examples/python/README.rst ) diff --git a/examples/python/contact-mechanics/CMakeLists.txt b/examples/python/contact-mechanics/CMakeLists.txt new file mode 100644 index 000000000..40d1c1390 --- /dev/null +++ b/examples/python/contact-mechanics/CMakeLists.txt @@ -0,0 +1,7 @@ +add_mesh(compression_mesh ./mesh/compression.geo DIM 2 ORDER 1) + +register_example(compression + SCRIPT compression.py + PYTHON + FILES_TO_COPY ./materials/compression.dat + DEPENDS compression_mesh) diff --git a/examples/python/contact-mechanics/compression.py b/examples/python/contact-mechanics/compression.py new file mode 100755 index 000000000..85b995411 --- /dev/null +++ b/examples/python/contact-mechanics/compression.py @@ -0,0 +1,69 @@ +#!/usr/bin/env python3 + +import akantu as aka + + +max_steps = 20000 +max_displacement = 1e-2 + +damping_interval = 10 +damping_ratio = 0.9 + +spatial_dimension = 2 +aka.parseInput('compression.dat') + +mesh = aka.Mesh(spatial_dimension) +mesh.read('compression.msh') + +coupler = aka.CouplerSolidContact(mesh) + +solid = coupler.getSolidMechanicsModel() +contact = coupler.getContactMechanicsModel() + +material_selector = aka.MeshDataMaterialSelectorString("physical_names", solid) +solid.setMaterialSelector(material_selector) + +coupler.initFull(_analysis_method=aka._explicit_lumped_mass) + +surface_selector = aka.PhysicalSurfaceSelector(mesh) +detector = contact.getContactDetector() +detector.setSurfaceSelector(surface_selector) + +solid.applyBC(aka.FixedValue(0.0, aka._x), "sides") + +time_step = solid.getStableTimeStep() +time_step *= 0.1 +coupler.setTimeStep(time_step) + +coupler.setBaseName("compression") +coupler.addDumpFieldVector("displacement") +coupler.addDumpFieldVector("contact_force") +coupler.addDumpFieldVector("external_force") +coupler.addDumpFieldVector("internal_force") +coupler.addDumpField("gaps") +coupler.addDumpField("areas") +coupler.addDumpField("blocked_dofs") +coupler.addDumpField("grad_u") +coupler.addDumpField("stress") + +coupler.dump() + +velocity = solid.getVelocity() + +increment = max_displacement / max_steps + + +for s in range(0, max_steps): + + print("Step : ", s) + + solid.applyBC(aka.IncrementValue(-increment, aka._y), "loading") + solid.applyBC(aka.IncrementValue(increment, aka._y), "fixed") + + coupler.solveStep() + + if s % damping_interval == 0: + velocity *= damping_ratio + + if s % 100 == 0: + coupler.dump() diff --git a/examples/python/contact-mechanics/detection-explicit.py b/examples/python/contact-mechanics/detection-explicit.py new file mode 100755 index 000000000..bd72d87af --- /dev/null +++ b/examples/python/contact-mechanics/detection-explicit.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python3 + +import akantu as aka +import time + + +spatial_dimension = 2 +aka.parseInput('./materials/detection-explicit.dat') + +mesh = aka.Mesh(spatial_dimension) +mesh.read('./mesh/detection-explicit.msh') + +model = aka.ContactMechanicsModel(mesh) + +model.initFull(_analysis_method=aka._explicit_lumped_mass) +surface_selector = aka.PhysicalSurfaceSelector(mesh) +model.getContactDetector().setSurfaceSelector(surface_selector) + + +model.setBaseName("detection-explicit") +model.addDumpFieldVector("normals") +model.addDumpField("gaps") +model.addDumpField("areas") + + +start_time = time.time() +model.search() +finish_time = time.time() +print('Search time = %s seconds', finish_time - start_time) + +model.dump() + +# by default the contact model creates a group named contact_surface +contact_surface = mesh.getElementGroup("contact_surface") diff --git a/examples/python/contact-mechanics/materials/compression.dat b/examples/python/contact-mechanics/materials/compression.dat new file mode 100644 index 000000000..8b81c3e76 --- /dev/null +++ b/examples/python/contact-mechanics/materials/compression.dat @@ -0,0 +1,30 @@ +material elastic [ + name = upper + rho = 1230 # density + E = 7e7 # young's modulu + nu = 0.0 # poisson's ratio +] + +material elastic [ + name = lower + rho = 1230 # density + E = 7e7 # young's modulu + nu = 0.0 # poisson's ratio +] + +contact_detector [ + type = explicit + max_iterations = 1000 + projection_tolerance = 1e-10 + extension_tolerance = 1e-3 + master = contact + slave = contact +] + +contact_resolution penalty_linear [ + name = contact + mu = 0.0 + epsilon_n = 1050e7 + epsilon_t = 0 + is_master_deformable = true +] \ No newline at end of file diff --git a/examples/python/contact-mechanics/materials/detection-explicit.dat b/examples/python/contact-mechanics/materials/detection-explicit.dat new file mode 100644 index 000000000..c2a46682c --- /dev/null +++ b/examples/python/contact-mechanics/materials/detection-explicit.dat @@ -0,0 +1,33 @@ +material elastic [ + name = upper + rho = 1.0 # density + E = 1.0 # young's modulu + nu = 0.0 # poisson's ratio +] + +material elastic [ + name = lower + rho = 1.0 # density + E = 1.0 # young's modulu + nu = 0.0 # poisson's ratio +] + + +contact_detector [ + type = explicit + slave = contact_top + master = contact_bottom + max_iterations = 1000 + projection_tolerance = 1e-2 + extension_tolerance = 1e-1 +] + + +contact_resolution penalty_linear [ + name = contact_bottom + mu = 0.0 + epsilon_n = 1050e7 + epsilon_t = 0 + is_master_deformable = true +] + diff --git a/examples/python/contact-mechanics/mesh/compression.geo b/examples/python/contact-mechanics/mesh/compression.geo new file mode 100644 index 000000000..b930bff0f --- /dev/null +++ b/examples/python/contact-mechanics/mesh/compression.geo @@ -0,0 +1,47 @@ +nb_points = 11; +nb_points_height = 11; +height = 2.0; +length = 1.0; + +Point(1) = {-length/2, -height/2, 0, 1}; +Point(2) = {length/2, -height/2, 0, 1}; +Point(3) = {length/2, 0, 0, 1}; +Point(4) = {length/2, height/2, 0, 1}; +Point(5) = {-length/2, height/2, 0, 1}; +Point(6) = {-length/2, 0, 0, 1}; +Point(7) = {length/2, 0, 0, 1}; +Point(8) = {-length/2, 0, 0, 1}; + +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {7, 4}; +Line(4) = {4, 5}; +Line(5) = {5, 6}; +Line(6) = {6, 7}; +Line(7) = {8, 1}; +Line(8) = {3, 8}; + +Line Loop(1) = {1, 2, 8, 7}; +Line Loop(2) = {3, 4, 5, 6}; + +Plane Surface(1) = {1}; +Plane Surface(2) = {2}; + +Physical Line ("loading") = {4}; +Physical Line ("fixed") = {1}; + +Physical Line ("sides") = {2, 7, 3, 5}; + +Physical Line ("contact") = {6, 8}; + +Physical Surface("lower") = {1}; +Physical Surface("upper") = {2}; + +Transfinite Line {6,4} = nb_points Using Progression 1; +Transfinite Line {3,5} = nb_points_height Using Progression 1; + +Transfinite Line {1, 8} = nb_points Using Progression 1; +Transfinite Line {7, 2} = nb_points_height Using Progression 1; + +Transfinite Surface "*"; +Recombine Surface "*"; diff --git a/examples/python/contact-mechanics/mesh/detection-explicit.msh b/examples/python/contact-mechanics/mesh/detection-explicit.msh new file mode 100644 index 000000000..c009e610d --- /dev/null +++ b/examples/python/contact-mechanics/mesh/detection-explicit.msh @@ -0,0 +1,1384 @@ +$MeshFormat +4 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 1 "contact_bottom" +1 2 "contact_top" +1 3 "top" +2 4 "top_body" +2 5 "bot_body" +$EndPhysicalNames +$Entities +8 8 2 0 +1 0 0.05 0 0 0.05 0 0 +2 0.5 0.5499999999999999 0 0.5 0.5499999999999999 0 0 +3 -0.5 0.5499999999999999 0 -0.5 0.5499999999999999 0 0 +4 0.5 0.1 0 0.5 0.1 0 0 +5 -0.5 0.1 0 -0.5 0.1 0 0 +6 -0.5 -0.25 0 -0.5 -0.25 0 0 +7 0.5 -0.25 0 0.5 -0.25 0 0 +8 0 0.5499999999999999 0 0 0.5499999999999999 0 0 +1 -0.5 0.04999999999999993 0 -2.775557561562891e-17 0.5499999999999999 0 1 2 2 3 -1 +2 0 0.04999999999999999 0 0.4999999999999999 0.5499999999999999 0 1 2 2 1 -2 +3 0 0.5499999999999999 0 0.5 0.5499999999999999 0 1 3 2 2 -8 +4 -0.5 -0.25 0 0.5 -0.25 0 0 2 6 -7 +5 0.5 -0.25 0 0.5 0.09999999999999998 0 0 2 7 -4 +6 -0.5 0.1 0 0.5 0.1 0 1 1 2 4 -5 +7 -0.5 -0.25 0 -0.5 0.1 0 0 2 5 -6 +13 -0.5 0.5499999999999999 0 0 0.5499999999999999 0 1 3 2 8 -3 +9 -0.5 0.04999999999999993 0 0.5 0.5499999999999999 0 1 4 4 2 3 13 1 +11 -0.5 -0.25 0 0.5 0.09999999999999999 0 1 5 4 6 7 4 5 +$EndEntities +$Nodes +18 453 +1 0 0 1 +1 0 0.05 0 +2 0 0 1 +2 0.5 0.5499999999999999 0 +3 0 0 1 +3 -0.5 0.5499999999999999 0 +4 0 0 1 +4 0.5 0.1 0 +5 0 0 1 +5 -0.5 0.1 0 +6 0 0 1 +6 -0.5 -0.25 0 +7 0 0 1 +7 0.5 -0.25 0 +8 0 0 1 +8 0 0.5499999999999999 0 +1 1 0 15 +9 -0.4975923633227955 0.5009914297001528 0 +10 -0.4903926401443037 0.4524548387038111 0 +11 -0.4784701677367776 0.4048576609464355 0 +12 -0.4619397660209172 0.3586582832507761 0 +13 -0.4409606318379142 0.3143016309578966 0 +14 -0.4157348056713792 0.2722148827719875 0 +15 -0.3865052260291442 0.2328033571234396 0 +16 -0.3535533897738377 0.1964466085872901 0 +17 -0.3171966413112021 0.1634947726861994 0 +18 -0.2777851158009608 0.1342651933750953 0 +19 -0.2356983677510476 0.1090393674720021 0 +20 -0.1913417156075182 0.0880602335061727 0 +21 -0.1451423381776323 0.07152983199751117 0 +22 -0.0975451607028337 0.05960735973767062 0 +23 -0.04900857001285892 0.05240763664893855 0 +2 1 0 15 +24 0.04900857028081373 0.05240763667532983 0 +25 0.09754516126712158 0.0596073598499145 0 +26 0.1451423390194117 0.07152983225286225 0 +27 0.1913417167183027 0.08806023396627477 0 +28 0.2356983690491223 0.1090393681658376 0 +29 0.2777851172325131 0.1342651943316281 0 +30 0.3171966428861189 0.1634947739787003 0 +31 0.3535533914263695 0.1964466102398221 0 +32 0.3865052273304663 0.2328033587091052 0 +33 0.4157348066433621 0.272214884226663 0 +34 0.4409606325469363 0.3143016322843838 0 +35 0.4619397665248591 0.3586582844673996 0 +36 0.4784701680176695 0.4048576618724118 0 +37 0.4903926402734724 0.4524548393531866 0 +38 0.4975923633581155 0.5009914300587632 0 +3 1 0 9 +39 0.4499999999997918 0.5499999999999999 0 +40 0.3999999999999999 0.5499999999999999 0 +41 0.3500000000003467 0.5499999999999999 0 +42 0.3000000000006934 0.5499999999999999 0 +43 0.2500000000010293 0.5499999999999999 0 +44 0.2000000000008322 0.5499999999999999 0 +45 0.1500000000006241 0.5499999999999999 0 +46 0.1000000000004161 0.5499999999999999 0 +47 0.05000000000020799 0.5499999999999999 0 +4 1 0 19 +48 -0.4500000000001387 -0.25 0 +49 -0.4000000000002774 -0.25 0 +50 -0.3500000000004161 -0.25 0 +51 -0.3000000000005548 -0.25 0 +52 -0.2500000000006935 -0.25 0 +53 -0.2000000000008322 -0.25 0 +54 -0.1500000000009709 -0.25 0 +55 -0.1000000000011096 -0.25 0 +56 -0.05000000000124821 -0.25 0 +57 -1.376121439022882e-12 -0.25 0 +58 0.04999999999875171 -0.25 0 +59 0.09999999999889042 -0.25 0 +60 0.1499999999990291 -0.25 0 +61 0.1999999999991677 -0.25 0 +62 0.2499999999993064 -0.25 0 +63 0.2999999999994452 -0.25 0 +64 0.349999999999584 -0.25 0 +65 0.3999999999997227 -0.25 0 +66 0.4499999999998613 -0.25 0 +5 1 0 6 +67 0.5 -0.2000000000001232 0 +68 0.5 -0.1500000000002334 0 +69 0.5 -0.1000000000003698 0 +70 0.5 -0.05000000000037691 0 +71 0.5 -2.481903571549537e-13 0 +72 0.5 0.04999999999986504 0 +6 1 0 19 +73 0.4500000000001387 0.1 0 +74 0.4000000000002774 0.1 0 +75 0.3500000000004161 0.1 0 +76 0.3000000000005548 0.1 0 +77 0.2500000000006935 0.1 0 +78 0.2000000000008322 0.1 0 +79 0.1500000000009709 0.1 0 +80 0.1000000000011096 0.1 0 +81 0.05000000000124821 0.1 0 +82 1.376121439022882e-12 0.1 0 +83 -0.04999999999875171 0.1 0 +84 -0.09999999999889042 0.1 0 +85 -0.1499999999990291 0.1 0 +86 -0.1999999999991677 0.1 0 +87 -0.2499999999993064 0.1 0 +88 -0.2999999999994452 0.1 0 +89 -0.349999999999584 0.1 0 +90 -0.3999999999997227 0.1 0 +91 -0.4499999999998613 0.1 0 +7 1 0 6 +92 -0.5 0.05000000000009655 0 +93 -0.5 2.364219930939271e-13 0 +94 -0.5 -0.04999999999962099 0 +95 -0.5 -0.09999999999961617 0 +96 -0.5 -0.1499999999997472 0 +97 -0.5 -0.1999999999998628 0 +13 1 0 9 +98 -0.04999999999990733 0.5499999999999999 0 +99 -0.09999999999977893 0.5499999999999999 0 +100 -0.1499999999996332 0.5499999999999999 0 +101 -0.1999999999994874 0.5499999999999999 0 +102 -0.2499999999993471 0.5499999999999999 0 +103 -0.2999999999994734 0.5499999999999999 0 +104 -0.349999999999605 0.5499999999999999 0 +105 -0.3999999999997367 0.5499999999999999 0 +106 -0.4499999999998684 0.5499999999999999 0 +9 2 0 186 +107 -0.01228171240411089 0.3 0 +108 0.1945144116810613 0.3406251720943104 0 +109 -0.2100211298736208 0.3479833165134979 0 +110 0.1107452138798006 0.2162994548598184 0 +111 -0.1350490637990908 0.2221111138161371 0 +112 0.324541486012991 0.4065716448655622 0 +113 0.07563202650358743 0.4070850503712184 0 +114 -0.3339119951447812 0.4121734463348341 0 +115 -0.09201428431920435 0.4124129772614709 0 +116 -0.01180070048916525 0.1674642745142007 0 +117 0.2193780858772615 0.2241296897461216 0 +118 -0.2416407840271332 0.2410241998348156 0 +119 0.2990957701151897 0.29855856263609 0 +120 0.2223758008676434 0.4436443083908998 0 +121 -0.2356425947502926 0.4468018204506074 0 +122 -0.3131792956270039 0.3095669645533252 0 +123 0.09967546028892153 0.3148806833669398 0 +124 -0.1184427524766596 0.3183726346464444 0 +125 -0.008093344874798035 0.4558508202432131 0 +126 0.4032167088347607 0.4519639813861305 0 +127 0.06937118358063861 0.1441659822533401 0 +128 -0.4002626322533028 0.4615894007816893 0 +129 -0.09665482234202258 0.144033914896221 0 +130 0.1402435753877111 0.4655567158327233 0 +131 -0.1548934646812773 0.4669251805240694 0 +132 0.1745753020806016 0.163813027362911 0 +133 0.02897950961398314 0.2369351905556829 0 +134 0.003354096229198298 0.3775614896635333 0 +135 0.3683693295434359 0.3422908872656255 0 +136 -0.1980518741330708 0.1768605338794096 0 +137 -0.06382438314857168 0.2423164770832472 0 +138 0.2914632700279634 0.4747255340464032 0 +139 0.1612790492680848 0.2696466445251161 0 +140 -0.3803041732620199 0.3557382360148033 0 +141 0.2648879504324043 0.363800339721242 0 +142 0.1462044069552536 0.3922371125369243 0 +143 -0.1815962890744263 0.2769641035850443 0 +144 -0.301631180746838 0.4773160773813773 0 +145 -0.07504944590892793 0.4824228984362947 0 +146 -0.1630709987320259 0.3942939329180792 0 +147 0.2917295765859801 0.2289485757901872 0 +148 0.06913465691766513 0.4837988795546873 0 +149 -0.2787850708742034 0.3684469006891219 0 +150 -0.05980966434146976 0.3526285953524161 0 +151 -0.306620622417532 0.2432539043280366 0 +152 0.4052314301002875 0.405987604507889 0 +153 0.005077833844681695 0.114319383425068 0 +154 0.2403156806056027 0.2928621088938413 0 +155 -0.4161277787241805 0.3984273761910028 0 +156 -0.2582686862921693 0.3030655372542683 0 +157 0.1243590543871349 0.1377164830005601 0 +158 0.03129083890539174 0.3307097571933929 0 +159 0.2318702664956497 0.1716503883539014 0 +160 0.3489232304948642 0.4904397445050183 0 +161 -0.1552245197395196 0.1361178931713556 0 +162 0.3546398939464878 0.2766610791941949 0 +163 -0.05544737687452236 0.1067820028232205 0 +164 0.1914672628197577 0.4926731984268204 0 +165 -0.2737884614488015 0.1807299392913903 0 +166 -0.3548614877054251 0.4936396560739741 0 +167 -0.2044703640879302 0.4937897639166191 0 +168 -0.06630056455616884 0.2964060959279643 0 +169 0.4421032423046113 0.4947186103636957 0 +170 0.07945973633154368 0.261983605144952 0 +171 0.2003859640860983 0.3915426663849307 0 +172 0.2732926077748499 0.4214362891352958 0 +173 -0.371498862707666 0.2759742494358789 0 +174 0.1674430872397316 0.2198257262659796 0 +175 -0.03627007952766946 0.4147500353758257 0 +176 0.02735238095022039 0.5036797611761729 0 +177 0.3140527631741167 0.3527041033257431 0 +178 -0.4468447453609086 0.485186750265691 0 +179 -0.2137422031336154 0.394777210689763 0 +180 -0.2851845884598839 0.4258284209351755 0 +181 -0.06889495248577712 0.1879352566346309 0 +182 0.241359553770711 0.4979948530762517 0 +183 -0.1891528093320655 0.2302159276774366 0 +184 0.02777227602836244 0.1876175650483655 0 +185 0.1496059325977905 0.3399109787823761 0 +186 -0.1114072310706673 0.2675957613872883 0 +187 0.1004570049009593 0.3673864559135394 0 +188 -0.2543117753532098 0.4970931351000933 0 +189 -0.1667530273723027 0.3447260254456604 0 +190 -0.326972317888319 0.3598257991467358 0 +191 -0.1100982083475947 0.3695715757974858 0 +192 -0.1185708021044796 0.5047712745288343 0 +193 -0.02941953603323009 0.503940035234691 0 +194 0.02907937137290616 0.4197749866986372 0 +195 -0.01753974569809826 0.2165683560093989 0 +196 0.4439452748483976 0.433504139555073 0 +197 0.0277178652975977 0.2820704148120675 0 +198 0.3359923710227309 0.2356579264164653 0 +199 0.1756494444739442 0.4333847507093272 0 +200 0.3673317891759262 0.4235243573371964 0 +201 0.1017377854145348 0.1762100592310798 0 +202 -0.4450220407661828 0.4396769380422502 0 +203 0.2815726051317162 0.1860136926341269 0 +204 0.1105267143033833 0.5037301973139618 0 +205 0.05757325127151662 0.0978960014684036 0 +206 -0.1193230393143431 0.1844935393379391 0 +207 -0.1870293889015087 0.4368043420057113 0 +208 0.3907417755193248 0.5030363681752912 0 +209 -0.3749670592881499 0.4251372716655361 0 +210 -0.2356225695306847 0.1561083701096387 0 +211 -0.1059620371437915 0.4548904508513228 0 +212 0.4107493575122539 0.3578334456834684 0 +213 -0.1110742835373288 0.1028779134318098 0 +214 0.05244082199947542 0.3738691890389321 0 +215 -0.3975596684208971 0.3149273776094628 0 +216 0.1184941291811034 0.4238802984251188 0 +217 0.3288890669436891 0.4507858599143252 0 +218 0.06967761747481055 0.2038686428488731 0 +219 0.1982620479953664 0.1308101206812429 0 +220 -0.2005700406018551 0.1329656962422475 0 +221 -0.3977987956789985 0.5064078943662431 0 +222 0.1917262833156118 0.3007730067533394 0 +223 0.3921527937307541 0.3118606045519516 0 +224 -0.1648493502554533 0.5111628719170078 0 +225 -0.3389720916500273 0.4533313052911815 0 +226 -0.4200226021135264 0.3501639992795211 0 +227 -0.01899712660445669 0.3417105517259422 0 +228 -0.2077974896461788 0.3102432745771253 0 +229 0.2584295519499785 0.2545688927705138 0 +230 0.2682249419737601 0.3250371358655978 0 +231 -0.3124651686122938 0.2089633043933055 0 +232 0.3580905708489841 0.3823643416166114 0 +233 0.1238182220635272 0.2602443086291486 0 +234 0.3128513997432721 0.5138458816257248 0 +235 0.04060689587032788 0.4620639250985911 0 +236 -0.05180534377204095 0.1508710042993477 0 +237 -0.3193239435988624 0.5159735157797003 0 +238 0.1537473528313092 0.5090495980408375 0 +239 -0.07008044985833065 0.4468920441821373 0 +240 -0.026018263064329 0.2617703187820156 0 +241 0.1371845619812296 0.1850587574478819 0 +242 -0.02164711176291285 0.08406911043651009 0 +243 -0.09656813580532106 0.2196630411697524 0 +244 0.1613853793157926 0.1148446249021088 0 +245 -0.230557103565537 0.2026163307598432 0 +246 -0.3429802247290291 0.2387686276816876 0 +247 -0.2735602925875388 0.2651881681037174 0 +248 -0.3538422058256955 0.3166191725060327 0 +249 0.1983166854233111 0.1948547079322284 0 +250 0.338020068541857 0.3174888172560016 0 +251 -0.2847428114020354 0.3334032915559064 0 +252 -0.1306340875652953 0.4199775586273106 0 +253 0.2430378020655228 0.3932500324798395 0 +254 -0.159360261263618 0.1899598015764556 0 +255 0.1409420805855188 0.299809633978901 0 +256 -0.3664566648614901 0.3902604258705824 0 +257 0.02306312221128798 0.08053936361388506 0 +258 0.2766437175103308 0.5192432235211231 0 +259 0.07406360642317401 0.5199055137132259 0 +260 0.315873473233897 0.201487329995932 0 +261 -0.1514209796878465 0.2594423208359149 0 +262 -0.2570323150100888 0.3960389484760207 0 +263 0.2340429247633811 0.3414625156743125 0 +264 0.2571228081102919 0.4594502461622126 0 +265 -0.07499999999984314 0.5208357786184671 0 +266 0.02959876884720683 0.1422926413418756 0 +267 0.2031223750187628 0.2594051164465031 0 +268 0.1061641361134195 0.09850712120078124 0 +269 -0.1553019278656079 0.3035430000089227 0 +270 0.08556141498564572 0.4472151778202449 0 +271 -0.2498107532304026 0.3447078011151668 0 +272 0.06259509898165284 0.3060683709821248 0 +273 0.3079774327240734 0.2588790073614902 0 +274 -0.2731054054431393 0.2236293077851847 0 +275 0.441264416588948 0.3918464879686415 0 +276 -0.2184207398950682 0.2679798337053588 0 +277 0.1736544012544477 0.3673798712550967 0 +278 0.2902661315161378 0.3881028633152414 0 +279 0.2561344563932663 0.2139690037148108 0 +280 -0.2701821834023651 0.4591776446922826 0 +281 -0.3269469173157441 0.2748951811014464 0 +282 -0.2271955348599936 0.5227207247541781 0 +283 -0.1886290992714705 0.3710930004576891 0 +284 -0.03559165492463761 0.3748767952833069 0 +285 0.06929184501528016 0.3385828912989858 0 +286 -0.04180303110106959 0.4675857195929832 0 +287 -0.06637119329666961 0.386982234826536 0 +288 0.367820514193713 0.4639500622635924 0 +289 -0.08508457843681064 0.3284624454893579 0 +290 -0.3019109991004751 0.3931926511838181 0 +291 0.4674227995907099 0.5193065053354683 0 +292 -0.2850533799395462 0.5180765456522342 0 +11 2 0 161 +293 0.3249999999994945 -0.07500000000037334 0 +294 -0.02500000000131217 -0.0749999999998034 0 +295 -0.2762557122470346 -0.07335747568371341 0 +296 0.140927771481751 -0.07341321072321773 0 +297 -0.1441709400256675 -0.1319763767042312 0 +298 -0.3802163758887959 -0.1312697847879755 0 +299 -0.3707605829622759 -0.01269315908636637 0 +300 -0.1140799469975067 -0.01045213698631014 0 +301 0.2374999999992321 -0.008928571428739188 0 +302 0.239364587879715 -0.1449604228290726 0 +303 0.05683783668245218 -0.008801746536017062 0 +304 0.06583658483893899 -0.1455653842333428 0 +305 0.4058294203883725 0.003589839066427533 0 +306 0.4003176320690112 -0.1530607430235512 0 +307 -0.2131563825668621 0.005510011661840052 0 +308 -0.06574628780476942 -0.1609089153461749 0 +309 -0.2292510822997157 -0.163487062701872 0 +310 -0.01832301366378245 0.006118770606085574 0 +311 0.1458627224865235 0.0168324318136457 0 +312 0.1546894720466331 -0.1681829946989501 0 +313 -0.2998132363646029 0.02266120706859185 0 +314 0.3165426249152323 0.01811371638284878 0 +315 0.318095053984824 -0.1661592608461226 0 +316 -0.3123680408613039 -0.1712912911155654 0 +317 0.4261585365853502 -0.07500000000005624 0 +318 -0.4266925381952101 -0.07026217120839512 0 +319 -0.1970322777185727 -0.06782600942511278 0 +320 -0.4320993185776716 0.0346675590610601 0 +321 0.00492987484335233 -0.1831612038955783 0 +322 -0.4314583333333957 -0.1814583333334978 0 +323 -0.09435967819805501 -0.07642117668779931 0 +324 -0.1343261681383029 -0.183845300734195 0 +325 0.04369533527558502 -0.07500000000002563 0 +326 0.2563046647225432 -0.07500000000007111 0 +327 -0.3349648284605178 -0.07859782167149161 0 +328 -0.1637412547543429 0.0351617204242195 0 +329 -0.08255931020141968 0.04275203207335465 0 +330 -0.3702405164610144 -0.1809766195890308 0 +331 0.022575426266819 0.05161118379743541 0 +332 0.439216502731838 0.04429705662903685 0 +333 0.439119270412986 -0.1931114074733934 0 +334 0.2072241187967513 0.02949383807627853 0 +335 0.2039310462924435 -0.1876630333200203 0 +336 -0.3510057112607998 0.03465013901818974 0 +337 0.09815398314724696 0.04260141803923802 0 +338 0.09797445496060024 -0.1916683154896894 0 +339 -0.4492463717601048 -0.1253673555886963 0 +340 0.2652364975922613 0.04440354390811885 0 +341 0.2638009276412374 -0.1932347881627331 0 +342 -0.187506185179097 -0.1933787542274065 0 +343 -0.001119162789662051 -0.1158805233532229 0 +344 -0.0625746317432373 -0.02885420713300732 0 +345 0.3661557686616563 0.04880199593503241 0 +346 0.36570850892365 -0.1999844062607133 0 +347 0.2031523323606292 -0.05441602316641909 0 +348 0.09358051248362761 -0.09481039509601544 0 +349 0.1824065017231087 -0.1224384099942984 0 +350 -0.1452332833423528 -0.05718131276243767 0 +351 0.01780693008609372 -0.03999348697789763 0 +352 0.1152986824354105 -0.03450762086480957 0 +353 0.4494149872321017 -0.1212303837234156 0 +354 0.4524753165016159 -0.02967898329634851 0 +355 0.2906523323612117 -0.03117880226484608 0 +356 0.2844255108021391 -0.1135049715785271 0 +357 -0.2502100884663495 0.04241024194873935 0 +358 0.3783986585742644 -0.04619491847853411 0 +359 -0.4501147516717464 -0.01918773315224366 0 +360 0.348430867163 -0.1240691328370529 0 +361 -0.3089378795849138 -0.1183602247142645 0 +362 -0.2625377292513435 -0.01193770861134554 0 +363 -0.2719656490672786 -0.2020505814655382 0 +364 -0.3190922547876296 -0.03038514645452571 0 +365 -0.2011880350881818 -0.1186802430453308 0 +366 -0.1671481950778641 -0.01242998255354861 0 +367 -0.08666973111362178 -0.205646082517727 0 +368 0.04744735209038599 -0.1960831933374882 0 +369 -0.05621899124359844 -0.1099098804451636 0 +370 -0.1007791885774214 -0.1249481545394985 0 +371 0.1898820016336944 -0.01005297727902277 0 +372 0.3341173222252349 -0.02540162162023579 0 +373 -0.2545277041095648 -0.1203151887467775 0 +374 -0.2049619686057672 0.05479017640626124 0 +375 -0.03171629155463005 0.05150881976605286 0 +376 0.3879534469372036 -0.09909252967716388 0 +377 0.1414224860497475 -0.1161839709791613 0 +378 0.1338111374855427 0.06572298254391738 0 +379 0.1364671391125252 -0.2112910739870971 0 +380 0.3248537089306673 -0.2124068254013839 0 +381 0.3248537089308958 0.06240682540109627 0 +382 -0.1256815842421477 0.0611869758964487 0 +383 -0.03021860521877257 -0.2056068229455497 0 +384 -0.3929885467522186 0.05788346913814332 0 +385 -0.3818886385532236 -0.09052207022668085 0 +386 -0.3220820294845056 -0.2116461645562594 0 +387 -0.3142780707348832 0.06338505821889698 0 +388 -0.2277445833095234 -0.2117832796789633 0 +389 0.1107007020759095 -0.1432822120994318 0 +390 -0.1768269806045011 -0.1556325999393932 0 +391 0.2237019123365921 -0.09898397129084736 0 +392 -0.02009477569053965 -0.1547056303909107 0 +393 0.161062916458809 -0.03806316871277764 0 +394 -0.335167946343827 0.003558260136472377 0 +395 0.07223919230762489 -0.04359482012557379 0 +396 -0.4599985525069109 0.05999855250695798 0 +397 0.4664431753448862 -0.07509279255330195 0 +398 -0.343404293198881 -0.1422918915503993 0 +399 -0.4668715641122466 -0.07446986984293799 0 +400 -0.4604062001595937 -0.2104062001595361 0 +401 0.4046190340108923 0.05967078944665156 0 +402 0.4046190340109789 -0.2096707894466728 0 +403 -0.4053845452928674 -0.02889635899479814 0 +404 -0.2750107656807784 -0.1562188115270072 0 +405 -0.2305678346876351 -0.083134138244025 0 +406 0.1765583357985406 0.06597806545887611 0 +407 0.1761694691664686 -0.2177676153468385 0 +408 0.2774828637169844 0.005602471649345589 0 +409 0.2759002528978573 -0.1546403713943836 0 +410 -0.4052919077668761 -0.2109388137977913 0 +411 0.2282360732359924 0.06604803389212706 0 +412 0.2312035150989076 -0.2153698752499626 0 +413 -0.2238354442748217 -0.03549177431919057 0 +414 0.06274006251134574 0.06291719483494408 0 +415 0.4149491582882755 -0.03723616540989799 0 +416 0.3600491603314121 0.005307369499360825 0 +417 -0.4685053825153637 0.02233621826892482 0 +418 -0.4690027838002592 -0.1728617496926251 0 +419 0.0987585099425011 0.002828429020932451 0 +420 0.02686643505682871 0.004456393497077993 0 +421 -0.4123717661690934 -0.1463428917904274 0 +422 -0.4003922427529299 0.0110706526639975 0 +423 -0.1250000000010403 -0.2214815877879124 0 +424 0.4630477848227073 0.01387860461642073 0 +425 0.4638761431913198 -0.1608533972700924 0 +426 -0.3731305647086207 -0.05189278794037629 0 +427 -0.0619908517320216 0.008570932119441562 0 +428 -0.02631524049744779 -0.03770030720663005 0 +429 -0.2797345912752921 0.0683060955971492 0 +430 0.01201371173015106 -0.2146979071018821 0 +431 0.4693641240713498 0.06306574888427212 0 +432 0.4680842701547646 -0.2145106106168785 0 +433 0.3581380155351213 -0.16081838574186 0 +434 -0.165363539682183 -0.09122443647606679 0 +435 0.2022667162102881 -0.1525301795835149 0 +436 0.07385545176215709 -0.2219378772067944 0 +437 0.02550482335310124 -0.1518964321424977 0 +438 -0.4270172835672771 0.07050991614123228 0 +439 -0.4100831381132856 -0.112752854720435 0 +440 -0.1720236727659839 0.07101540252779293 0 +441 -0.1306420582546562 0.02324372177083282 0 +442 -0.06017825070573003 -0.06140099152802066 0 +443 -0.129981325965136 -0.0963502914340067 0 +444 0.01255456453922094 -0.08452274858821118 0 +445 0.04278282525557234 -0.1109714333529252 0 +446 -0.09625844788148934 -0.04217114490828659 0 +447 0.06154776135442938 0.026720852746579 0 +448 0.2374999999991322 -0.04464029426307854 0 +449 -0.3521970742577616 -0.110039933743747 0 +450 -0.1597166817299364 -0.2178819288723885 0 +451 0.177945462428461 -0.08404931907798031 0 +452 -0.1013296760148425 -0.1635039782492083 0 +453 0.2906523323610561 -0.06655049226162012 0 +$EndNodes +$Elements +7 868 +1 1 1 16 +9 3 9 +10 9 10 +11 10 11 +12 11 12 +13 12 13 +14 13 14 +15 14 15 +16 15 16 +17 16 17 +18 17 18 +19 18 19 +20 19 20 +21 20 21 +22 21 22 +23 22 23 +24 23 1 +2 1 1 16 +25 1 24 +26 24 25 +27 25 26 +28 26 27 +29 27 28 +30 28 29 +31 29 30 +32 30 31 +33 31 32 +34 32 33 +35 33 34 +36 34 35 +37 35 36 +38 36 37 +39 37 38 +40 38 2 +3 1 1 10 +41 2 39 +42 39 40 +43 40 41 +44 41 42 +45 42 43 +46 43 44 +47 44 45 +48 45 46 +49 46 47 +50 47 8 +6 1 1 20 +78 4 73 +79 73 74 +80 74 75 +81 75 76 +82 76 77 +83 77 78 +84 78 79 +85 79 80 +86 80 81 +87 81 82 +88 82 83 +89 83 84 +90 84 85 +91 85 86 +92 86 87 +93 87 88 +94 88 89 +95 89 90 +96 90 91 +97 91 5 +13 1 1 10 +105 8 98 +106 98 99 +107 99 100 +108 100 101 +109 101 102 +110 102 103 +111 103 104 +112 104 105 +113 105 106 +114 106 3 +9 2 2 422 +115 153 116 236 +116 44 164 182 +117 39 40 169 +118 22 23 163 +119 22 163 213 +120 119 154 229 +121 222 185 255 +122 106 9 178 +123 169 40 208 +124 108 185 222 +125 133 197 240 +126 153 205 266 +127 11 12 155 +128 9 106 3 +129 163 153 236 +130 150 168 227 +131 20 21 161 +132 43 44 182 +133 106 178 221 +134 168 107 227 +135 138 172 217 +136 189 228 269 +137 127 201 218 +138 105 106 221 +139 184 127 218 +140 126 152 196 +141 189 109 228 +142 28 29 159 +143 119 229 273 +144 205 127 266 +145 166 128 225 +146 37 38 169 +147 205 153 257 +148 172 112 217 +149 180 144 225 +150 128 209 225 +151 178 10 202 +152 175 134 194 +153 206 181 243 +154 37 169 196 +155 125 175 194 +156 156 122 247 +157 129 181 206 +158 161 21 213 +159 181 116 195 +160 137 181 195 +161 9 10 178 +162 130 204 270 +163 41 160 208 +164 166 104 221 +165 114 180 225 +166 195 133 240 +167 152 126 200 +168 154 119 230 +169 136 161 254 +170 120 164 199 +171 32 33 162 +172 116 153 266 +173 164 130 199 +174 107 158 227 +175 157 127 268 +176 164 120 182 +177 158 134 227 +178 212 152 232 +179 122 156 251 +180 135 212 232 +181 157 132 241 +182 163 129 213 +183 194 113 270 +184 104 105 221 +185 235 194 270 +186 137 168 186 +187 109 179 271 +188 28 159 219 +189 179 262 271 +190 187 142 216 +191 123 170 233 +192 247 122 281 +193 155 12 226 +194 116 184 195 +195 184 133 195 +196 161 206 254 +197 17 18 165 +198 113 187 216 +199 131 207 252 +200 229 154 267 +201 164 44 238 +202 11 155 202 +203 110 174 233 +204 204 148 270 +205 140 155 226 +206 155 128 202 +207 123 233 255 +208 40 41 208 +209 117 229 267 +210 155 140 256 +211 121 167 188 +212 129 161 213 +213 168 150 289 +214 127 157 201 +215 128 155 209 +216 132 157 244 +217 159 29 203 +218 207 146 252 +219 144 180 280 +220 169 126 196 +221 167 121 207 +222 131 167 207 +223 171 108 263 +224 197 107 240 +225 158 107 197 +226 253 171 263 +227 159 132 219 +228 127 205 268 +229 174 139 233 +230 47 8 176 +231 172 138 264 +232 44 45 238 +233 128 178 202 +234 156 247 276 +235 168 124 186 +236 160 41 234 +237 222 154 263 +238 32 162 198 +239 134 158 214 +240 26 244 268 +241 14 15 173 +242 16 231 246 +243 139 222 255 +244 176 8 193 +245 110 170 218 +246 170 133 218 +247 45 46 204 +248 228 143 269 +249 142 199 216 +250 244 157 268 +251 191 150 287 +252 47 176 259 +253 161 129 206 +254 150 191 289 +255 133 170 197 +256 247 118 276 +257 162 119 273 +258 108 222 263 +259 119 177 230 +260 137 186 243 +261 132 159 249 +262 177 141 230 +263 185 142 187 +264 125 176 193 +265 117 159 279 +266 186 124 269 +267 159 117 249 +268 101 167 224 +269 199 130 216 +270 261 186 269 +271 20 161 220 +272 160 138 217 +273 201 157 241 +274 161 136 220 +275 162 33 223 +276 175 125 286 +277 8 98 193 +278 123 185 187 +279 138 160 234 +280 99 100 192 +281 162 223 250 +282 223 135 250 +283 142 171 199 +284 119 162 250 +285 126 169 208 +286 239 175 286 +287 176 148 259 +288 165 18 210 +289 168 137 240 +290 107 168 240 +291 146 189 191 +292 129 163 236 +293 111 206 243 +294 174 110 241 +295 231 151 246 +296 45 204 238 +297 183 111 261 +298 112 177 232 +299 189 124 191 +300 171 120 199 +301 163 23 242 +302 153 163 242 +303 196 152 275 +304 182 120 264 +305 130 164 238 +306 169 38 291 +307 179 146 207 +308 29 30 203 +309 120 171 253 +310 172 120 253 +311 39 169 291 +312 124 168 289 +313 216 130 270 +314 17 165 231 +315 187 113 214 +316 128 166 221 +317 104 166 237 +318 167 101 282 +319 177 135 232 +320 100 101 224 +321 144 166 225 +322 179 121 262 +323 36 37 196 +324 121 180 262 +325 41 42 234 +326 167 131 224 +327 111 183 254 +328 166 144 237 +329 31 32 198 +330 184 116 266 +331 146 191 252 +332 192 131 211 +333 145 192 211 +334 190 122 251 +335 149 190 251 +336 194 134 214 +337 113 194 214 +338 154 222 267 +339 14 173 215 +340 121 179 207 +341 178 128 221 +342 10 11 202 +343 143 183 261 +344 34 35 212 +345 186 111 243 +346 121 188 280 +347 115 191 287 +348 191 115 252 +349 170 110 233 +350 156 228 271 +351 191 124 289 +352 177 112 278 +353 126 208 288 +354 208 160 288 +355 192 100 224 +356 198 162 273 +357 170 123 272 +358 108 171 277 +359 212 35 275 +360 24 25 205 +361 171 142 277 +362 228 156 276 +363 112 172 278 +364 181 129 236 +365 132 174 241 +366 190 114 256 +367 151 247 281 +368 201 110 218 +369 174 117 267 +370 120 172 264 +371 139 174 267 +372 182 138 258 +373 228 109 271 +374 112 200 217 +375 152 200 232 +376 136 183 245 +377 133 184 218 +378 183 118 245 +379 13 14 215 +380 18 19 210 +381 176 125 235 +382 209 155 256 +383 229 147 273 +384 173 15 246 +385 21 22 213 +386 190 140 248 +387 122 190 248 +388 211 131 252 +389 16 17 231 +390 219 132 244 +391 117 174 249 +392 115 175 239 +393 248 173 281 +394 122 248 281 +395 134 175 284 +396 148 176 235 +397 174 132 249 +398 175 115 287 +399 19 20 220 +400 138 182 264 +401 15 16 246 +402 140 190 256 +403 152 212 275 +404 141 177 278 +405 135 177 250 +406 177 119 250 +407 192 145 265 +408 116 181 236 +409 145 193 265 +410 99 192 265 +411 193 98 265 +412 118 183 276 +413 165 210 245 +414 33 34 223 +415 183 143 276 +416 114 190 290 +417 141 253 263 +418 180 114 290 +419 179 109 283 +420 181 137 243 +421 26 27 244 +422 146 179 283 +423 103 104 237 +424 27 28 219 +425 209 114 225 +426 131 192 224 +427 197 170 272 +428 159 203 279 +429 165 245 274 +430 245 118 274 +431 183 136 254 +432 210 136 245 +433 25 26 268 +434 148 235 270 +435 188 102 292 +436 180 121 280 +437 253 141 278 +438 172 253 278 +439 36 196 275 +440 13 215 226 +441 43 182 258 +442 187 214 285 +443 188 144 280 +444 125 193 286 +445 12 13 226 +446 127 184 266 +447 185 123 255 +448 102 103 292 +449 185 108 277 +450 262 149 271 +451 215 140 226 +452 142 185 277 +453 125 194 235 +454 211 115 239 +455 154 230 263 +456 145 211 239 +457 23 1 242 +458 111 186 261 +459 144 188 292 +460 24 205 257 +461 204 130 238 +462 2 39 291 +463 31 198 260 +464 34 212 223 +465 136 210 220 +466 198 147 260 +467 123 187 285 +468 149 262 290 +469 38 2 291 +470 262 180 290 +471 210 19 220 +472 190 149 290 +473 124 189 269 +474 188 167 282 +475 137 195 240 +476 102 188 282 +477 109 189 283 +478 189 146 283 +479 143 261 269 +480 233 139 255 +481 200 112 232 +482 212 135 223 +483 1 24 257 +484 42 43 258 +485 46 47 259 +486 30 31 260 +487 98 99 265 +488 193 145 286 +489 110 201 241 +490 138 234 258 +491 203 30 260 +492 147 203 260 +493 35 36 275 +494 101 102 282 +495 217 200 288 +496 148 204 259 +497 204 46 259 +498 206 111 254 +499 203 147 279 +500 158 197 272 +501 147 198 273 +502 215 173 248 +503 242 1 257 +504 153 242 257 +505 205 25 268 +506 214 158 285 +507 140 215 248 +508 115 211 252 +509 27 219 244 +510 114 209 256 +511 200 126 288 +512 251 156 271 +513 113 216 270 +514 237 144 292 +515 222 139 267 +516 147 229 279 +517 160 217 288 +518 230 141 263 +519 234 42 258 +520 231 165 274 +521 151 231 274 +522 143 228 276 +523 247 151 274 +524 118 247 274 +525 229 117 279 +526 150 227 284 +527 227 134 284 +528 173 246 281 +529 246 151 281 +530 272 123 285 +531 158 272 285 +532 149 251 271 +533 145 239 286 +534 103 237 292 +535 284 175 287 +536 150 284 287 +11 2 2 374 +537 83 84 329 +538 307 357 362 +539 357 313 362 +540 297 324 452 +541 358 293 376 +542 76 340 381 +543 341 63 380 +544 302 356 391 +545 356 326 391 +546 334 301 340 +547 302 335 341 +548 295 327 361 +549 293 360 376 +550 327 295 364 +551 343 294 369 +552 348 296 352 +553 76 77 340 +554 62 63 341 +555 370 297 452 +556 83 329 375 +557 324 297 390 +558 343 369 392 +559 81 82 331 +560 50 330 410 +561 49 50 410 +562 319 350 366 +563 332 305 424 +564 306 333 425 +565 303 351 395 +566 351 325 395 +567 364 299 426 +568 319 366 413 +569 305 354 424 +570 353 306 425 +571 299 364 394 +572 74 75 345 +573 64 65 346 +574 308 383 392 +575 87 357 374 +576 331 82 375 +577 338 304 368 +578 359 403 422 +579 320 359 422 +580 353 317 376 +581 306 353 376 +582 340 314 381 +583 315 341 380 +584 322 330 421 +585 330 298 421 +586 74 345 401 +587 346 65 402 +588 86 87 374 +589 336 89 384 +590 348 352 395 +591 375 329 427 +592 311 334 406 +593 335 312 407 +594 305 332 401 +595 333 306 402 +596 328 307 366 +597 295 405 413 +598 307 328 374 +599 293 356 360 +600 356 315 360 +601 310 375 427 +602 95 96 339 +603 362 295 413 +604 341 335 412 +605 334 340 411 +606 310 331 375 +607 332 73 401 +608 66 333 402 +609 51 363 386 +610 335 302 435 +611 345 305 401 +612 306 346 402 +613 330 316 398 +614 73 332 431 +615 333 66 432 +616 327 364 426 +617 329 300 427 +618 298 330 398 +619 51 52 363 +620 378 311 406 +621 312 379 407 +622 338 59 379 +623 80 337 378 +624 383 321 392 +625 93 94 359 +626 330 322 410 +627 309 342 390 +628 308 367 383 +629 89 336 387 +630 358 305 416 +631 70 71 354 +632 68 69 353 +633 309 365 373 +634 89 90 384 +635 81 331 414 +636 316 330 386 +637 300 329 441 +638 59 338 436 +639 330 50 386 +640 329 84 382 +641 372 358 416 +642 334 311 371 +643 301 334 371 +644 295 362 364 +645 359 94 399 +646 312 335 435 +647 368 304 437 +648 331 310 420 +649 342 309 388 +650 363 316 386 +651 400 48 410 +652 366 307 413 +653 373 365 405 +654 73 74 401 +655 65 66 402 +656 300 344 427 +657 365 309 390 +658 95 339 399 +659 339 318 399 +660 342 324 390 +661 367 324 423 +662 350 300 366 +663 337 80 414 +664 362 313 364 +665 313 357 429 +666 318 359 399 +667 348 304 389 +668 355 293 372 +669 324 367 452 +670 337 311 378 +671 312 338 379 +672 70 354 397 +673 353 69 397 +674 48 49 410 +675 293 358 372 +676 299 336 422 +677 336 313 387 +678 351 294 444 +679 322 339 418 +680 339 96 418 +681 311 337 419 +682 294 343 444 +683 348 325 445 +684 304 348 445 +685 59 60 379 +686 79 80 378 +687 338 312 389 +688 349 312 435 +689 336 299 394 +690 313 336 394 +691 357 307 374 +692 304 338 389 +693 339 322 421 +694 309 363 388 +695 354 317 397 +696 317 353 397 +697 355 301 448 +698 347 301 371 +699 88 89 387 +700 344 300 446 +701 325 351 444 +702 301 347 448 +703 325 348 395 +704 347 393 451 +705 393 296 451 +706 312 349 377 +707 317 358 376 +708 342 53 450 +709 347 391 448 +710 377 349 451 +711 296 377 451 +712 318 339 439 +713 391 326 448 +714 361 316 404 +715 62 341 412 +716 340 77 411 +717 302 341 409 +718 340 301 408 +719 314 355 372 +720 387 313 429 +721 53 342 388 +722 310 344 428 +723 350 323 446 +724 300 350 446 +725 314 340 408 +726 341 315 409 +727 310 351 420 +728 326 355 448 +729 296 348 377 +730 5 92 396 +731 91 5 396 +732 294 351 428 +733 345 75 381 +734 64 346 380 +735 314 345 381 +736 346 315 380 +737 67 425 432 +738 424 72 431 +739 324 342 450 +740 344 310 427 +741 351 310 428 +742 346 306 433 +743 6 48 400 +744 97 6 400 +745 295 361 373 +746 301 355 408 +747 356 302 409 +748 323 350 443 +749 369 323 370 +750 55 56 367 +751 293 355 453 +752 355 326 453 +753 356 293 453 +754 326 356 453 +755 352 311 419 +756 311 352 393 +757 315 346 433 +758 305 345 416 +759 345 314 416 +760 360 306 376 +761 350 319 434 +762 407 61 412 +763 78 406 411 +764 403 318 426 +765 318 385 426 +766 366 300 441 +767 351 303 420 +768 328 366 441 +769 352 296 393 +770 68 353 425 +771 354 71 424 +772 335 407 412 +773 406 334 411 +774 373 361 404 +775 303 395 419 +776 395 352 419 +777 425 333 432 +778 332 424 431 +779 354 305 415 +780 317 354 415 +781 56 57 383 +782 355 314 408 +783 315 356 409 +784 377 348 389 +785 93 359 417 +786 359 320 417 +787 308 369 370 +788 67 68 425 +789 71 72 424 +790 80 81 414 +791 84 85 382 +792 391 349 435 +793 302 391 435 +794 405 319 413 +795 322 400 410 +796 357 87 429 +797 360 315 433 +798 305 358 415 +799 358 317 415 +800 82 83 375 +801 347 371 393 +802 363 52 388 +803 359 318 403 +804 55 367 423 +805 57 58 430 +806 316 361 398 +807 53 54 450 +808 361 327 449 +809 58 368 430 +810 306 360 433 +811 369 308 392 +812 75 76 381 +813 63 64 380 +814 321 368 437 +815 365 319 405 +816 343 392 437 +817 392 321 437 +818 363 309 404 +819 316 363 404 +820 367 56 383 +821 367 308 452 +822 370 323 443 +823 382 85 440 +824 307 362 413 +825 52 53 388 +826 50 51 386 +827 364 313 394 +828 403 299 422 +829 297 365 390 +830 61 62 412 +831 77 78 411 +832 72 4 431 +833 4 73 431 +834 66 7 432 +835 7 67 432 +836 371 311 393 +837 328 382 440 +838 365 297 434 +839 319 365 434 +840 69 70 397 +841 94 95 399 +842 336 384 422 +843 383 57 430 +844 391 347 451 +845 92 93 417 +846 312 377 389 +847 297 370 443 +848 96 97 418 +849 369 294 442 +850 323 369 442 +851 60 61 407 +852 78 79 406 +853 314 372 416 +854 85 86 440 +855 87 88 429 +856 54 55 423 +857 58 59 436 +858 329 382 441 +859 90 91 438 +860 368 58 436 +861 368 321 430 +862 309 373 404 +863 338 368 436 +864 295 373 405 +865 86 374 440 +866 321 383 430 +867 79 378 406 +868 379 60 407 +869 308 370 452 +870 320 384 438 +871 385 318 439 +872 374 328 440 +873 398 361 449 +874 384 320 422 +875 396 92 417 +876 349 391 451 +877 385 327 426 +878 298 385 439 +879 97 400 418 +880 384 90 438 +881 299 403 426 +882 385 298 449 +883 327 385 449 +884 396 320 438 +885 382 328 441 +886 419 337 447 +887 303 419 447 +888 88 387 429 +889 343 437 445 +890 437 304 445 +891 423 324 450 +892 320 396 417 +893 323 442 446 +894 442 344 446 +895 294 428 442 +896 428 344 442 +897 400 322 418 +898 434 297 443 +899 350 434 443 +900 331 420 447 +901 420 303 447 +902 444 343 445 +903 325 444 445 +904 91 396 438 +905 298 398 449 +906 421 298 439 +907 337 414 447 +908 414 331 447 +909 339 421 439 +910 54 423 450 +$EndElements diff --git a/packages/contact_mechanics.cmake b/packages/contact_mechanics.cmake new file mode 100644 index 000000000..f15c8affd --- /dev/null +++ b/packages/contact_mechanics.cmake @@ -0,0 +1,73 @@ +#=============================================================================== +# @file contact_mechanics.cmake +# +# @author Mohit Pundir +# +# @date creation: Sun Oct 21 2018 +# @date last modification: Sun Oct 21 2018 +# +# @brief package description for contact mechanics +# +# @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 . +# +#=============================================================================== +package_declare(contact_mechanics + DEPENDS model_couplers cohesive_element + DESCRIPTION "Use Contact Mechanics package of Akantu") + +package_declare_sources(contact_mechanics + model/contact_mechanics/contact_mechanics_model.hh + model/contact_mechanics/contact_mechanics_model.cc + model/contact_mechanics/contact_detector.hh + model/contact_mechanics/contact_detector.cc + model/contact_mechanics/contact_detector_inline_impl.cc + model/contact_mechanics/contact_element.hh + model/contact_mechanics/geometry_utils.hh + model/contact_mechanics/geometry_utils.cc + model/contact_mechanics/geometry_utils_inline_impl.cc + + model/contact_mechanics/resolution.hh + model/contact_mechanics/resolution.cc + model/contact_mechanics/resolution_utils.hh + model/contact_mechanics/resolution_utils.cc + model/contact_mechanics/resolutions/resolution_penalty.hh + model/contact_mechanics/resolutions/resolution_penalty.cc + model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh + model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc + + model/contact_mechanics/surface_selector.hh + model/contact_mechanics/surface_selector.cc + + model/model_couplers/coupler_solid_contact.hh + model/model_couplers/coupler_solid_contact_tmpl.hh + model/model_couplers/coupler_solid_contact.cc + model/model_couplers/coupler_solid_cohesive_contact.hh + model/model_couplers/coupler_solid_cohesive_contact.cc + model/model_couplers/cohesive_contact_solvercallback.hh + model/model_couplers/cohesive_contact_solvercallback.cc + ) + +package_declare_documentation_files(contact_mechanics + manual-contactmechanicsmodel.tex + manual-contact-detector.tex + ) + +package_declare_documentation(contact_mechanics + "This package activates the contact mechanics model") diff --git a/packages/core.cmake b/packages/core.cmake index 78762cc20..ca7272e73 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,522 +1,528 @@ #=============================================================================== # @file core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for core # # @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 . # #=============================================================================== package_declare(core NOT_OPTIONAL DESCRIPTION "core package for Akantu" FEATURES_PUBLIC cxx_strong_enums cxx_defaulted_functions cxx_deleted_functions cxx_auto_type cxx_decltype_auto FEATURES_PRIVATE cxx_lambdas cxx_nullptr cxx_range_for cxx_delegating_constructors DEPENDS INTERFACE akantu_iterators Boost) package_declare_sources(core common/aka_array.cc common/aka_array.hh + common/aka_array_filter.hh common/aka_array_tmpl.hh common/aka_array_printer.hh common/aka_bbox.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.hh common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.hh common/aka_csr.hh common/aka_element_classes_info_inline_impl.hh common/aka_enum_macros.hh common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_factory.hh common/aka_fwd.hh common/aka_grid_dynamic.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_named_argument.hh common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_typelist.hh common/aka_types.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper_tmpl.hh common/aka_voigthelper.cc common/aka_warning.hh common/aka_warning_restore.hh fe_engine/element_class.hh + fe_engine/element_class_helper.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.hh fe_engine/element_classes/element_class_hexahedron_20_inline_impl.hh fe_engine/element_classes/element_class_pentahedron_6_inline_impl.hh fe_engine/element_classes/element_class_pentahedron_15_inline_impl.hh fe_engine/element_classes/element_class_point_1_inline_impl.hh fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh fe_engine/element_classes/element_class_quadrangle_8_inline_impl.hh fe_engine/element_classes/element_class_segment_2_inline_impl.hh fe_engine/element_classes/element_class_segment_3_inline_impl.hh fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.hh fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.hh fe_engine/element_classes/element_class_triangle_3_inline_impl.hh fe_engine/element_classes/element_class_triangle_6_inline_impl.hh fe_engine/element_type_conversion.hh fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.hh fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl_field.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_element_property.hh fe_engine/geometrical_element_property.cc fe_engine/gauss_integration.cc fe_engine/gauss_integration_tmpl.hh fe_engine/integrator.hh fe_engine/integrator_gauss.hh fe_engine/integrator_gauss_inline_impl.hh fe_engine/interpolation_element_tmpl.hh fe_engine/integration_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions.cc fe_engine/shape_functions_inline_impl.hh fe_engine/shape_lagrange_base.cc fe_engine/shape_lagrange_base.hh fe_engine/shape_lagrange_base_inline_impl.hh fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.hh fe_engine/element.hh io/dumper/dumpable.hh io/dumper/dumpable.cc io/dumper/dumpable_dummy.hh io/dumper/dumpable_inline_impl.hh io/dumper/dumper_field.hh io/dumper/dumper_material_padders.hh io/dumper/dumper_filtered_connectivity.hh io/dumper/dumper_element_partition.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh #io/model_io.cc #io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parser.cc io/parser/parser_real.cc io/parser/parser_random.cc io/parser/parser_types.cc io/parser/parser_input_files.cc io/parser/parser.hh io/parser/parser_tmpl.hh io/parser/parser_grammar_tmpl.hh io/parser/cppargparse/cppargparse.hh io/parser/cppargparse/cppargparse.cc io/parser/cppargparse/cppargparse_tmpl.hh io/parser/parameter_registry.cc io/parser/parameter_registry.hh io/parser/parameter_registry_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.hh mesh/element_type_map.cc mesh/element_type_map.hh mesh/element_type_map_tmpl.hh mesh/element_type_map_filter.hh mesh/group_manager.cc mesh/group_manager.hh mesh/group_manager_inline_impl.hh mesh/mesh.cc mesh/mesh.hh mesh/mesh_periodic.cc mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_global_data_updater.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.hh mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.hh mesh/mesh_iterators.hh mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_utils_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_distribution.cc mesh_utils/mesh_utils_distribution.hh mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.hh mesh_utils/global_ids_updater.hh mesh_utils/global_ids_updater.cc mesh_utils/global_ids_updater_inline_impl.hh model/common/boundary_condition/boundary_condition.hh model/common/boundary_condition/boundary_condition_functor.hh model/common/boundary_condition/boundary_condition_functor_inline_impl.hh model/common/boundary_condition/boundary_condition_tmpl.hh model/common/non_local_toolbox/neighborhood_base.hh model/common/non_local_toolbox/neighborhood_base.cc model/common/non_local_toolbox/neighborhood_base_inline_impl.hh model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.hh model/common/non_local_toolbox/non_local_manager.hh model/common/non_local_toolbox/non_local_manager.cc model/common/non_local_toolbox/non_local_manager_inline_impl.hh model/common/non_local_toolbox/non_local_manager_callback.hh model/common/non_local_toolbox/non_local_neighborhood_base.hh model/common/non_local_toolbox/non_local_neighborhood_base.cc model/common/non_local_toolbox/non_local_neighborhood.hh model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh model/common/non_local_toolbox/non_local_neighborhood_inline_impl.hh model/common/non_local_toolbox/base_weight_function.hh model/common/non_local_toolbox/base_weight_function_inline_impl.hh model/common/model_solver.cc model/common/model_solver.hh model/common/solver_callback.hh model/common/solver_callback.cc model/common/dof_manager/dof_manager.cc model/common/dof_manager/dof_manager.hh model/common/dof_manager/dof_manager_default.cc model/common/dof_manager/dof_manager_default.hh model/common/dof_manager/dof_manager_default_inline_impl.hh model/common/dof_manager/dof_manager_inline_impl.hh model/common/non_linear_solver/non_linear_solver.cc model/common/non_linear_solver/non_linear_solver.hh model/common/non_linear_solver/non_linear_solver_default.hh model/common/non_linear_solver/non_linear_solver_lumped.cc model/common/non_linear_solver/non_linear_solver_lumped.hh model/common/time_step_solvers/time_step_solver.hh model/common/time_step_solvers/time_step_solver.cc model/common/time_step_solvers/time_step_solver_default.cc model/common/time_step_solvers/time_step_solver_default.hh model/common/time_step_solvers/time_step_solver_default_explicit.hh model/common/integration_scheme/generalized_trapezoidal.cc model/common/integration_scheme/generalized_trapezoidal.hh model/common/integration_scheme/integration_scheme.cc model/common/integration_scheme/integration_scheme.hh model/common/integration_scheme/integration_scheme_1st_order.cc model/common/integration_scheme/integration_scheme_1st_order.hh model/common/integration_scheme/integration_scheme_2nd_order.cc model/common/integration_scheme/integration_scheme_2nd_order.hh model/common/integration_scheme/newmark-beta.cc model/common/integration_scheme/newmark-beta.hh model/common/integration_scheme/pseudo_time.cc model/common/integration_scheme/pseudo_time.hh model/model.cc model/model.hh model/model_inline_impl.hh model/model_options.hh solver/solver_vector.hh solver/solver_vector_default.cc solver/solver_vector_default.hh solver/solver_vector_default_tmpl.hh solver/solver_vector_distributed.cc solver/solver_vector_distributed.hh solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_aij.cc solver/sparse_matrix_aij.hh solver/sparse_matrix_aij_inline_impl.hh solver/sparse_matrix_inline_impl.hh solver/sparse_solver.cc solver/sparse_solver.hh solver/sparse_solver_inline_impl.hh solver/terms_to_assemble.hh synchronizer/communication_buffer_inline_impl.hh synchronizer/communication_descriptor.hh synchronizer/communication_descriptor_tmpl.hh synchronizer/communication_request.hh synchronizer/communication_tag.hh synchronizer/communications.hh synchronizer/communications_tmpl.hh synchronizer/communicator.cc synchronizer/communicator.hh synchronizer/communicator_dummy_inline_impl.hh synchronizer/communicator_event_handler.hh synchronizer/communicator_inline_impl.hh synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.hh synchronizer/element_info_per_processor.cc synchronizer/element_info_per_processor.hh synchronizer/element_info_per_processor_tmpl.hh synchronizer/element_synchronizer.cc synchronizer/element_synchronizer.hh synchronizer/facet_synchronizer.cc synchronizer/facet_synchronizer.hh synchronizer/facet_synchronizer_inline_impl.hh synchronizer/grid_synchronizer.cc synchronizer/grid_synchronizer.hh synchronizer/grid_synchronizer_tmpl.hh synchronizer/master_element_info_per_processor.cc synchronizer/node_info_per_processor.cc synchronizer/node_info_per_processor.hh synchronizer/node_synchronizer.cc synchronizer/node_synchronizer.hh synchronizer/node_synchronizer_inline_impl.hh synchronizer/periodic_node_synchronizer.cc synchronizer/periodic_node_synchronizer.hh synchronizer/slave_element_info_per_processor.cc synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_impl.hh synchronizer/synchronizer_impl_tmpl.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh synchronizer/synchronizer_tmpl.hh synchronizer/communication_buffer.hh ) set(AKANTU_SPIRIT_SOURCES io/mesh_io/mesh_io_abaqus.cc io/parser/parser_real.cc io/parser/parser_random.cc io/parser/parser_types.cc io/parser/parser_input_files.cc PARENT_SCOPE ) package_declare_elements(core ELEMENT_TYPES _point_1 _segment_2 _segment_3 _triangle_3 _triangle_6 _quadrangle_4 _quadrangle_8 _tetrahedron_4 _tetrahedron_10 _pentahedron_6 _pentahedron_15 _hexahedron_8 _hexahedron_20 KIND regular GEOMETRICAL_TYPES _gt_point _gt_segment_2 _gt_segment_3 _gt_triangle_3 _gt_triangle_6 _gt_quadrangle_4 _gt_quadrangle_8 _gt_tetrahedron_4 _gt_tetrahedron_10 _gt_hexahedron_8 _gt_hexahedron_20 _gt_pentahedron_6 _gt_pentahedron_15 INTERPOLATION_TYPES _itp_lagrange_point_1 _itp_lagrange_segment_2 _itp_lagrange_segment_3 _itp_lagrange_triangle_3 _itp_lagrange_triangle_6 _itp_lagrange_quadrangle_4 _itp_serendip_quadrangle_8 _itp_lagrange_tetrahedron_4 _itp_lagrange_tetrahedron_10 _itp_lagrange_hexahedron_8 _itp_serendip_hexahedron_20 _itp_lagrange_pentahedron_6 _itp_lagrange_pentahedron_15 GEOMETRICAL_SHAPES _gst_point _gst_triangle _gst_square _gst_prism GAUSS_INTEGRATION_TYPES _git_point _git_segment _git_triangle _git_tetrahedron _git_pentahedron INTERPOLATION_KIND _itk_lagrangian FE_ENGINE_LISTS gradient_on_integration_points interpolate_on_integration_points interpolate compute_normals_on_integration_points inverse_map contains compute_shapes compute_shapes_derivatives + get_N + compute_dnds + compute_d2nds2 + compute_jmat get_shapes_derivatives lagrange_base assemble_fields ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-authors.tex manual-changelog.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-feengine.tex manual-elements.tex manual-appendix-elements.tex manual-appendix-packages.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg # figures/doc_wheel.pdf # figures/doc_wheel.svg figures/hot-point-1.png figures/hot-point-2.png figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/vectors.pdf figures/vectors.svg figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) package_declare_documentation(core "This package is the core engine of \\akantu. It depends on:" "\\begin{itemize}" "\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel})." "\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system." "\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries." "\\item The \\href{http://www.zlib.net/}{zlib} compression library." "\\end{itemize}" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install cmake libboost-dev zlib1g-dev g++" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{itemize}" "\\item Install Xcode" "\\item Install the command line tools." "\\item Install the MacPorts project which allows to automatically" "download and install opensource packages." "\\end{itemize}" "Then the following commands should be typed in a terminal:" "\\begin{command}" " > sudo port install cmake gcc48 boost" "\\end{command}" ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) find_program(PATCH_COMMAND patch) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) package_declare_extra_files_to_package(core SOURCES common/aka_element_classes_info.hh.in common/aka_config.hh.in ) if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang" AND (NOT CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.9)) package_set_compile_flags(core CXX "-Wno-undefined-var-template") endif() if(DEFINED AKANTU_CXX11_FLAGS) package_declare(core_cxx11 NOT_OPTIONAL DESCRIPTION "C++ 11 additions for Akantu core" COMPILE_FLAGS CXX "${AKANTU_CXX11_FLAGS}") if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "4.6") set(AKANTU_CORE_CXX11 OFF CACHE BOOL "C++ 11 additions for Akantu core - not supported by the selected compiler" FORCE) endif() endif() package_declare_documentation(core_cxx11 "This option activates some features of the C++11 standard. This is usable with GCC>=4.7 or Intel>=13.") else() if(CMAKE_VERSION VERSION_LESS 3.1) message(FATAL_ERROR "Since version 3.0 Akantu requires at least c++11 capable compiler") endif() endif() diff --git a/packages/damage_non_local.cmake b/packages/damage_non_local.cmake index bb3de2867..7e0ad5e09 100644 --- a/packages/damage_non_local.cmake +++ b/packages/damage_non_local.cmake @@ -1,66 +1,67 @@ #=============================================================================== # @file damage_non_local.cmake # # @author Nicolas Richart # # @date creation: Fri Jun 15 2012 # @date last modification: Mon Jan 18 2016 # # @brief package description for non-local materials # # @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 . # #=============================================================================== - package_declare(damage_non_local DESCRIPTION "Package for Non-local damage constitutives laws Akantu" DEPENDS lapack) package_declare_sources(damage_non_local model/solid_mechanics/materials/material_damage/material_damage_non_local.hh model/solid_mechanics/materials/material_damage/material_marigo_non_local.cc model/solid_mechanics/materials/material_damage/material_marigo_non_local.hh model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh - + model/solid_mechanics/materials/material_damage/material_von_mises_mazars_non_local.cc + model/solid_mechanics/materials/material_damage/material_von_mises_mazars_non_local.hh + model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.hh model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.hh model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.hh model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh model/solid_mechanics/materials/weight_functions/stress_based_weight_function.cc model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.hh ) package_declare_material_infos(damage_non_local LIST AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST INCLUDE material_non_local_includes.hh ) package_declare_documentation_files(damage_non_local manual-constitutive-laws-non_local.tex manual-appendix-materials-non-local.tex) package_declare_documentation(damage_non_local "This package activates the non local damage feature of AKANTU" "") diff --git a/packages/model_couplers.cmake b/packages/model_couplers.cmake index 47c375b8c..2b14e53a7 100644 --- a/packages/model_couplers.cmake +++ b/packages/model_couplers.cmake @@ -1,47 +1,32 @@ #=============================================================================== # @file model_couplers.cmake # # @author Mohit Pundir # # @date creation: Sun Sep 30 2018 # @date last modification: Sun Sep 28 2018 # # @brief package description for model couplers # # @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 . # #=============================================================================== - package_declare(model_couplers ADVANCED DESCRIPTION "Use Model Couplers package of Akantu") - -package_declare_sources(model_couplers - model/model_couplers/coupler_solid_phasefield.hh - model/model_couplers/coupler_solid_phasefield.cc - ) - -package_declare_documentation_files(model_couplers - # - ) - -package_declare_documentation(model_couplers - "This package activates the modle couplers within Akantu. " - "It has no additional dependencies." - ) diff --git a/packages/phase_field.cmake b/packages/phase_field.cmake index 01ab95a37..f9dbd8560 100644 --- a/packages/phase_field.cmake +++ b/packages/phase_field.cmake @@ -1,57 +1,60 @@ #=============================================================================== # @file phase_field.cmake # # @author Mohit Pundir # # @date creation: Sun Sep 30 2018 # @date last modification: Sun Sep 30 2018 # # @brief package description for phase field model # # @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 . # #=============================================================================== package_declare(phase_field DEPENDS model_couplers DESCRIPTION "Use Phase Field package of Akantu") package_declare_sources(phase_field model/phase_field/phasefield.cc model/phase_field/phasefield.hh model/phase_field/phasefield_inline_impl.cc model/phase_field/phasefield_selector.hh model/phase_field/phasefield_selector_tmpl.hh model/phase_field/phasefields/phasefield_exponential.hh model/phase_field/phasefields/phasefield_exponential.cc model/phase_field/phase_field_model.cc model/phase_field/phase_field_model.hh model/phase_field/phase_field_model_inline_impl.cc + + model/model_couplers/coupler_solid_phasefield.hh + model/model_couplers/coupler_solid_phasefield.cc ) package_declare_documentation_files(phase_field # ) package_declare_documentation(phase_field "This package activates the phase field model within Akantu. " ) diff --git a/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake index ae2d6ed5c..e2c363359 100644 --- a/packages/solid_mechanics.cmake +++ b/packages/solid_mechanics.cmake @@ -1,137 +1,143 @@ #=============================================================================== # @file solid_mechanics.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for core # # @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 . # #=============================================================================== package_declare(solid_mechanics DEFAULT ON DESCRIPTION "Solid mechanics model" DEPENDS core lapack ) package_declare_sources(solid_mechanics model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.hh model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.hh model/solid_mechanics/solid_mechanics_model_io.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.hh model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.cc model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.hh model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_phasefield.cc model/solid_mechanics/materials/material_damage/material_phasefield.hh model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.hh + model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc + model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh + model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.hh + model/solid_mechanics/materials/material_damage/material_von_mises_mazars.cc + model/solid_mechanics/materials/material_damage/material_von_mises_mazars.hh + model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_non_local_tmpl.hh model/solid_mechanics/materials/material_non_local_includes.hh ) package_declare_material_infos(solid_mechanics LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(solid_mechanics manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-appendix-materials.tex figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf ) package_declare_extra_files_to_package(solid_mechanics SOURCES model/solid_mechanics/material_list.hh.in ) diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index c4ef70381..f90ffb6bb 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,133 +1,147 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Dec 12 2014 # @date last modification: Mon Jan 18 2016 # # @brief CMake file for the python wrapping of akantu # # @section LICENSE # # Copyright (©) 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(NOT SKBUILD) package_get_all_include_directories( AKANTU_LIBRARY_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 ) endif() set(PYAKANTU_SRCS py_aka_common.cc py_aka_error.cc py_akantu.cc py_boundary_conditions.cc py_fe_engine.cc py_group_manager.cc py_mesh.cc py_model.cc py_parser.cc py_solver.cc ) package_is_activated(iohelper _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_dumpable.cc ) endif() package_is_activated(solid_mechanics _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model.cc py_material.cc py_material_selector.cc ) endif() package_is_activated(cohesive_element _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model_cohesive.cc py_fragment_manager.cc ) endif() package_is_activated(heat_transfer _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_heat_transfer_model.cc ) endif() + +package_is_activated(contact_mechanics _is_activated) +if(_is_activated) + list(APPEND PYAKANTU_SRCS + py_contact_mechanics_model.cc + ) +endif() + +package_is_activated(model_couplers _is_activated) +if(_is_activated) + list(APPEND PYAKANTU_SRCS + py_model_couplers.cc + ) +endif() + package_is_activated(phase_field _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_phase_field_model.cc ) endif() - package_is_activated(structural_mechanics _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_structural_mechanics_model.cc ) endif() pybind11_add_module(py11_akantu ${PYAKANTU_SRCS}) # to avoid compilation warnings from pybind11 target_include_directories(py11_akantu SYSTEM BEFORE PRIVATE ${PYBIND11_INCLUDE_DIR} PRIVATE ${pybind11_INCLUDE_DIR} PRIVATE ${PYTHON_INCLUDE_DIRS}) target_link_libraries(py11_akantu PUBLIC akantu) set_target_properties(py11_akantu PROPERTIES DEBUG_POSTFIX "" LIBRARY_OUTPUT_DIRECTORY akantu) file(COPY akantu DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) if(NOT SKBUILD) set(_python_install_dir ${CMAKE_INSTALL_LIBDIR}/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages) else() set(_python_install_dir .) endif() install(TARGETS py11_akantu LIBRARY DESTINATION ${_python_install_dir}) if(NOT SKBUILD) install(DIRECTORY akantu DESTINATION ${_python_install_dir} FILES_MATCHING PATTERN "*.py") endif() diff --git a/python/py_aka_common.cc b/python/py_aka_common.cc index 07fbcc881..22b86b372 100644 --- a/python/py_aka_common.cc +++ b/python/py_aka_common.cc @@ -1,112 +1,114 @@ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; namespace akantu { /* -------------------------------------------------------------------------- */ #define PY_AKANTU_PP_VALUE(s, data, elem) \ .value(BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data, elem)) #define PY_AKANTU_REGISTER_ENUM_(type_name, list, prefix, mod) \ py::enum_(mod, BOOST_PP_STRINGIZE(type_name)) \ BOOST_PP_SEQ_FOR_EACH(PY_AKANTU_PP_VALUE, prefix, list) \ .export_values() #define PY_AKANTU_REGISTER_CLASS_ENUM(type_name, list, mod) \ PY_AKANTU_REGISTER_ENUM_(type_name, list, type_name::_, mod) #define PY_AKANTU_REGISTER_ENUM(type_name, list, mod) \ PY_AKANTU_REGISTER_ENUM_(type_name, list, , mod) /* -------------------------------------------------------------------------- */ void register_initialize(py::module & mod) { mod.def("__initialize", []() { int nb_args = 0; char ** null = nullptr; initialize(nb_args, null); }); } void register_enums(py::module & mod) { py::enum_(mod, "SpatialDirection") .value("_x", _x) .value("_y", _y) .value("_z", _z) .export_values(); py::enum_(mod, "AnalysisMethod") .value("_static", _static) .value("_implicit_dynamic", _implicit_dynamic) .value("_explicit_lumped_mass", _explicit_lumped_mass) .value("_explicit_lumped_capacity", _explicit_lumped_capacity) .value("_explicit_consistent_mass", _explicit_consistent_mass) + .value("_explicit_contact", _explicit_contact) + .value("_implicit_contact", _implicit_contact) .export_values(); PY_AKANTU_REGISTER_CLASS_ENUM(ModelType, AKANTU_MODEL_TYPES, mod); PY_AKANTU_REGISTER_CLASS_ENUM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES, mod); PY_AKANTU_REGISTER_CLASS_ENUM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE, mod); PY_AKANTU_REGISTER_CLASS_ENUM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE, mod); PY_AKANTU_REGISTER_CLASS_ENUM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA, mod); py::enum_(mod, "CohesiveMethod") .value("_intrinsic", _intrinsic) .value("_extrinsic", _extrinsic) .export_values(); py::enum_(mod, "GhostType") .value("_not_ghost", _not_ghost) .value("_ghost", _ghost) .value("_casper", _casper) .export_values(); py::enum_(mod, "MeshIOType") .value("_miot_auto", _miot_auto) .value("_miot_gmsh", _miot_gmsh) .value("_miot_gmsh_struct", _miot_gmsh_struct) .value("_miot_diana", _miot_diana) .value("_miot_abaqus", _miot_abaqus) .export_values(); py::enum_(mod, "MatrixType") .value("_unsymmetric", _unsymmetric) .value("_symmetric", _symmetric) .export_values(); PY_AKANTU_REGISTER_ENUM(ElementType, AKANTU_ALL_ELEMENT_TYPE(_not_defined), mod); PY_AKANTU_REGISTER_ENUM(ElementKind, AKANTU_ELEMENT_KIND(_ek_not_defined), mod); } /* -------------------------------------------------------------------------- */ #define AKANTU_PP_STR_TO_TYPE2(s, data, elem) ({BOOST_PP_STRINGIZE(elem), elem}) void register_functions(py::module & mod) { mod.def("getElementTypes", []() { std::map element_types{ BOOST_PP_SEQ_FOR_EACH_I( AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(AKANTU_ek_regular_ELEMENT_TYPE), BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE2, akantu, AKANTU_ek_regular_ELEMENT_TYPE))}; return element_types; }); } #undef AKANTU_PP_STR_TO_TYPE2 } // namespace akantu diff --git a/python/py_akantu.cc b/python/py_akantu.cc index 614751c8c..8713b02d3 100644 --- a/python/py_akantu.cc +++ b/python/py_akantu.cc @@ -1,127 +1,137 @@ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_common.hh" #include "py_aka_error.hh" #include "py_boundary_conditions.hh" #include "py_fe_engine.hh" #include "py_group_manager.hh" #include "py_mesh.hh" #include "py_model.hh" #include "py_parser.hh" #include "py_solver.hh" #if defined(AKANTU_USE_IOHELPER) #include "py_dumpable.hh" #endif #if defined(AKANTU_SOLID_MECHANICS) #include "py_material.hh" #include "py_material_selector.hh" #include "py_solid_mechanics_model.hh" #endif #if defined(AKANTU_HEAT_TRANSFER) #include "py_heat_transfer_model.hh" #endif #if defined(AKANTU_COHESIVE_ELEMENT) #include "py_fragment_manager.hh" #include "py_solid_mechanics_model_cohesive.hh" #endif +#if defined(AKANTU_CONTACT_MECHANICS) +#include "py_contact_mechanics_model.hh" +#include "py_model_couplers.hh" +#endif + #if defined(AKANTU_PHASE_FIELD) #include "py_phase_field_model.hh" #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "py_structural_mechanics_model.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; namespace akantu { void register_all(pybind11::module & mod) { register_initialize(mod); register_enums(mod); register_error(mod); register_functions(mod); register_parser(mod); register_solvers(mod); register_group_manager(mod); #if defined(AKANTU_USE_IOHELPER) register_dumpable(mod); #endif register_mesh(mod); register_fe_engine(mod); register_boundary_conditions(mod); register_model(mod); #if defined(AKANTU_HEAT_TRANSFER) register_heat_transfer_model(mod); #endif #if defined(AKANTU_SOLID_MECHANICS) register_solid_mechanics_model(mod); register_material(mod); register_material_selector(mod); #endif #if defined(AKANTU_COHESIVE_ELEMENT) register_solid_mechanics_model_cohesive(mod); register_fragment_manager(mod); #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) register_structural_mechanics_model(mod); #endif +#if defined(AKANTU_CONTACT_MECHANICS) + register_contact_mechanics_model(mod); + register_model_couplers(mod); +#endif + #if defined(AKANTU_PHASE_FIELD) register_phase_field_model(mod); register_phase_field_coupler(mod); #endif } } // namespace akantu /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ PYBIND11_MODULE(py11_akantu, mod) { mod.doc() = "Akantu python interface"; static py::exception akantu_exception(mod, "Exception"); py::register_exception_translator([](std::exception_ptr ptr) { try { if (ptr) { std::rethrow_exception(ptr); } } catch (akantu::debug::Exception & e) { if (akantu::debug::debugger.printBacktrace()) { akantu::debug::printBacktrace(); } akantu_exception(e.info().c_str()); } }); akantu::register_all(mod); mod.def("has_mpi", []() { #if defined(AKANTU_USE_MPI) return true; #else return false; #endif }); } // Module akantu diff --git a/python/py_contact_mechanics_model.cc b/python/py_contact_mechanics_model.cc new file mode 100644 index 000000000..5f191a31e --- /dev/null +++ b/python/py_contact_mechanics_model.cc @@ -0,0 +1,130 @@ +/* -------------------------------------------------------------------------- */ +#include "py_aka_array.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +namespace py = pybind11; +/* -------------------------------------------------------------------------- */ + +namespace akantu { +/* -------------------------------------------------------------------------- */ +#define def_function_nocopy(func_name) \ + def( \ + #func_name, \ + [](ContactMechanicsModel & self) -> decltype(auto) { \ + return self.func_name(); \ + }, \ + py::return_value_policy::reference) + +#define def_function(func_name) \ + def(#func_name, [](ContactMechanicsModel & self) -> decltype(auto) { \ + return self.func_name(); \ + }) + +/* -------------------------------------------------------------------------- */ +void register_contact_mechanics_model(py::module & mod) { + py::class_(mod, "ContactDetector", + py::multiple_inheritance()) + .def(py::init(), py::arg("mesh"), + py::arg("id") = "contact_detector") + .def(py::init, const ID &>(), py::arg("mesh"), + py::arg("positions"), py::arg("id") = "contact_detector") + .def("setSurfaceSelector", &ContactDetector::setSurfaceSelector); + + py::class_>( + mod, "SurfaceSelector", py::multiple_inheritance()) + .def(py::init(), py::arg("mesh")); + + py::class_>( + mod, "PhysicalSurfaceSelector") + .def(py::init(), py::arg("mesh")); + + py::class_>( + mod, "CohesiveSurfaceSelector") + .def(py::init(), py::arg("mesh")); + + py::class_>(mod, "AllSurfaceSelector") + .def(py::init(), py::arg("mesh")); + + py::class_(mod, "ContactMechanicsModelOptions") + .def(py::init(), + py::arg("analysis_method") = _explicit_contact); + + /* ------------------------------------------------------------------------ */ + py::class_(mod, "ContactMechanicsModel", + py::multiple_inheritance()) + .def(py::init, + const ModelType>(), + py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, + py::arg("id") = "contact_mechanics_model", + py::arg("dof_manager") = nullptr, + py::arg("model_type") = ModelType::_contact_mechanics_model) + .def( + "initFull", + [](ContactMechanicsModel & self, + const ContactMechanicsModelOptions & options) { + self.initFull(options); + }, + py::arg("options") = ContactMechanicsModelOptions()) + .def( + "initFull", + [](ContactMechanicsModel & self, + const AnalysisMethod & analysis_method) { + self.initFull(_analysis_method = analysis_method); + }, + py::arg("_analysis_method")) + .def_function(search) + .def_function(assembleStiffnessMatrix) + .def_function(assembleInternalForces) + .def_function_nocopy(getExternalForce) + .def_function_nocopy(getNormalForce) + .def_function_nocopy(getTangentialForce) + .def_function_nocopy(getInternalForce) + .def_function_nocopy(getGaps) + .def_function_nocopy(getNormals) + .def_function_nocopy(getNodalArea) + .def_function_nocopy(getContactDetector); + + py::class_(mod, "ContactElement") + .def(py::init<>()) + .def_readwrite("master", &ContactElement::master) + .def_readwrite("slave", &ContactElement::slave); + + py::class_(mod, "GeometryUtils") + .def_static( + "normal", + py::overload_cast &, const Element &, + Vector &, bool>(&GeometryUtils::normal), + py::arg("mesh"), py::arg("positions"), py::arg("element"), + py::arg("normal"), py::arg("outward") = true) + .def_static( + "covariantBasis", + py::overload_cast &, const Element &, + const Vector &, Vector &, + Matrix &>(&GeometryUtils::covariantBasis), + py::arg("mesh"), py::arg("positions"), py::arg("element"), + py::arg("normal"), py::arg("natural_projection"), py::arg("basis")) + .def_static("curvature", &GeometryUtils::curvature) + .def_static("contravariantBasis", &GeometryUtils::contravariantBasis, + py::arg("covariant_basis"), py::arg("basis")) + .def_static("realProjection", + py::overload_cast &, + const Vector &, const Element &, + const Vector &, Vector &>( + &GeometryUtils::realProjection), + py::arg("mesh"), py::arg("positions"), py::arg("slave"), + py::arg("element"), py::arg("normal"), py::arg("projection")) + .def_static("isBoundaryElement", &GeometryUtils::isBoundaryElement); +} + +} // namespace akantu diff --git a/python/py_contact_mechanics_model.hh b/python/py_contact_mechanics_model.hh new file mode 100644 index 000000000..cec86c46a --- /dev/null +++ b/python/py_contact_mechanics_model.hh @@ -0,0 +1,10 @@ +#include + +#ifndef __AKANTU_PY_CONTACT_MECHANICS_MODEL_HH__ +#define __AKANTU_PY_CONTACT_MECHANICS_MODEL_HH__ + +namespace akantu { +void register_contact_mechanics_model(pybind11::module & mod); +} // namespace akantu + +#endif // __AKANTU_PY_CONTACT_MECHANICS_MODEL_HH__ diff --git a/python/py_fe_engine.cc b/python/py_fe_engine.cc index 0f0119fda..021539fc7 100644 --- a/python/py_fe_engine.cc +++ b/python/py_fe_engine.cc @@ -1,122 +1,120 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" #include "py_aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { -__attribute__((visibility("default"))) void -register_fe_engine(py::module & mod) { - +void register_fe_engine(py::module & mod) { py::class_(mod, "Element") .def(py::init([](ElementType type, UInt id) { return new Element{type, id, _not_ghost}; })) .def(py::init([](ElementType type, UInt id, GhostType ghost_type) { return new Element{type, id, ghost_type}; })) .def("__lt__", [](Element & self, const Element & other) { return (self < other); }) .def("__repr__", [](Element & self) { return std::to_string(self); }); mod.attr("ElementNull") = ElementNull; py::class_(mod, "FEEngine") .def( "getNbIntegrationPoints", [](FEEngine & fem, const ElementType & type, const GhostType & ghost_type) { return fem.getNbIntegrationPoints(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "gradientOnIntegrationPoints", [](FEEngine & fem, const Array & u, Array & nablauq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, const Array & u, Array & uq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } self.interpolateOnIntegrationPoints(u, uq, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, const Array & u, ElementTypeMapArray & uq, const ElementTypeMapArray * filter_elements) { self.interpolateOnIntegrationPoints(u, uq, filter_elements); }, py::arg("u"), py::arg("uq"), py::arg("filter_elements") = nullptr) .def( "computeIntegrationPointsCoordinates", [](FEEngine & self, ElementTypeMapArray & coordinates, const ElementTypeMapArray * filter_elements) -> decltype(auto) { return self.computeIntegrationPointsCoordinates(coordinates, filter_elements); }, py::arg("coordinates"), py::arg("filter_elements") = nullptr) .def( "assembleFieldLumped", [](FEEngine & fem, const std::function &, const Element &)> & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, GhostType ghost_type) { fem.assembleFieldLumped(field_funct, matrix_id, dof_id, dof_manager, type, ghost_type); }, py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"), py::arg("dof_manager"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "assembleFieldMatrix", [](FEEngine & fem, const std::function &, const Element &)> & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, GhostType ghost_type = _not_ghost) { fem.assembleFieldMatrix(field_funct, matrix_id, dof_id, dof_manager, type, ghost_type); }, py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"), py::arg("dof_manager"), py::arg("type"), py::arg("ghost_type") = _not_ghost); py::class_(mod, "IntegrationPoint"); } } // namespace akantu diff --git a/python/py_material.cc b/python/py_material.cc index b2cc09d94..fdc0528e0 100644 --- a/python/py_material.cc +++ b/python/py_material.cc @@ -1,177 +1,178 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ +#include #include #if defined(AKANTU_COHESIVE_ELEMENT) #include #endif #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { template class PyMaterial : public _Material { public: /* Inherit the constructors */ using _Material::_Material; ~PyMaterial() override = default; void initMaterial() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, initMaterial, ); }; void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, _Material, computeStress, el_type, ghost_type); } void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type, tangent_matrix, ghost_type); } void computePotentialEnergy(ElementType el_type) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type); } Real getPushWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element); } Real getShearWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element); } template void registerInternal(const std::string & name, UInt nb_component) { auto && internal = std::make_shared>(name, *this); AKANTU_DEBUG_INFO("alloc internal " << name << " " << &this->internals[name]); internal->initialize(nb_component); this->internals[name] = internal; } protected: std::map> internals; }; /* ------------------------------------------------------------------------ */ template void register_internal_field(py::module & mod, const std::string & name) { py::class_, ElementTypeMapArray, std::shared_ptr>>( mod, ("InternalField" + name).c_str()); } /* ------------------------------------------------------------------------ */ template void register_material_classes(py::module & mod, const std::string & name) { py::class_<_Material, Material, Parsable, PyMaterial<_Material>>( mod, name.c_str(), py::multiple_inheritance()) .def(py::init()); } } // namespace /* -------------------------------------------------------------------------- */ void register_material(py::module & mod) { py::class_(mod, "MaterialFactory") .def_static( "getInstance", []() -> MaterialFactory & { return Material::getFactory(); }, py::return_value_policy::reference) .def("registerAllocator", [](MaterialFactory & self, const std::string id, py::function func) { self.registerAllocator( id, [func, id](UInt dim, const ID & /*unused*/, SolidMechanicsModel & model, const ID & option) -> std::unique_ptr { py::object obj = func(dim, id, model, option); auto & ptr = py::cast(obj); obj.release(); return std::unique_ptr(&ptr); }); }) .def("getPossibleAllocators", &MaterialFactory::getPossibleAllocators); register_internal_field(mod, "Real"); register_internal_field(mod, "UInt"); py::class_>( mod, "Material", py::multiple_inheritance()) .def(py::init()) .def( "getGradU", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getGradU(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getStress", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getStress(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getPotentialEnergy", [](Material & self, ElementType el_type) -> decltype(auto) { return self.getPotentialEnergy(el_type); }, py::return_value_policy::reference) .def("initMaterial", &Material::initMaterial) .def("getModel", &Material::getModel) .def("registerInternalReal", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def("registerInternalUInt", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def( "getInternalReal", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getInternalUInt", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getElementFilter", [](Material & self) -> decltype(auto) { return self.getElementFilter(); }, py::return_value_policy::reference); register_material_classes>(mod, "MaterialElastic2D"); register_material_classes>(mod, "MaterialElastic3D"); } } // namespace akantu diff --git a/python/py_material.hh b/python/py_material.hh index 186538181..910a95843 100644 --- a/python/py_material.hh +++ b/python/py_material.hh @@ -1,12 +1,12 @@ #include #ifndef AKANTU_PY_MATERIAL_HH_ #define AKANTU_PY_MATERIAL_HH_ namespace akantu { void register_material(pybind11::module & mod); - +void register_material_selector(pybind11::module & mod); } // namespace akantu #endif // AKANTU_PY_MATERIAL_HH_ diff --git a/python/py_mesh.cc b/python/py_mesh.cc index f73aeb8fd..25ea2e8c4 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,150 +1,160 @@ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { /* ------------------------------------------------------------------------ */ template void register_element_type_map_array(py::module & mod, const std::string & name) { py::class_, std::shared_ptr>>( mod, ("ElementTypeMapArray" + name).c_str()) .def( "__call__", [](ElementTypeMapArray & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "elementTypes", [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, ElementKind _kind) -> std::vector { auto types = self.elementTypes(_dim, _ghost_type, _kind); std::vector _types; for (auto && t : types) { _types.push_back(t); } return _types; }, py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular); } } // namespace /* -------------------------------------------------------------------------- */ void register_mesh(py::module & mod) { register_element_type_map_array(mod, "Real"); register_element_type_map_array(mod, "UInt"); //register_element_type_map_array(mod, "String"); py::class_(mod, "MeshData") .def( "getElementalDataUInt", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference) .def( "getElementalDataReal", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference); py::class_(mod, "Mesh", py::multiple_inheritance()) .def(py::init(), py::arg("spatial_dimension"), py::arg("id") = "mesh") .def("read", &Mesh::read, py::arg("filename"), py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file") .def( "getNodes", [](Mesh & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("getNbNodes", &Mesh::getNbNodes) .def( "getConnectivity", [](Mesh & self, ElementType type) -> decltype(auto) { return self.getConnectivity(type); }, py::return_value_policy::reference) .def( "addConnectivityType", [](Mesh & self, ElementType type, GhostType ghost_type) -> void { self.addConnectivityType(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("distribute", [](Mesh & self) { self.distribute(); }) + .def("fillNodesToElements", &Mesh::fillNodesToElements, + py::arg("dimension") = _all_dimensions) + .def("getAssociatedElements", + [](Mesh & self, const UInt & node, py::list list) { + Array elements; + self.getAssociatedElements(node, elements); + for (auto && element : elements) { + list.append(element); + } + }) .def("makePeriodic", [](Mesh & self, const SpatialDirection & direction) { self.makePeriodic(direction); }) .def( "getNbElement", [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type, ElementKind kind) { return self.getNbElement(spatial_dimension, ghost_type, kind); }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined) .def( "getNbElement", [](Mesh & self, ElementType type, GhostType ghost_type) { return self.getNbElement(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def_static( "getSpatialDimension", [](ElementType & type) { return Mesh::getSpatialDimension(type); }) .def( "getDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> decltype(auto) { return _this.getData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "hasDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> bool { return _this.hasData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost); /* ------------------------------------------------------------------------ */ py::class_(mod, "MeshUtils") .def_static("buildFacets", &MeshUtils::buildFacets); py::class_(mod, "MeshAccessor") .def(py::init(), py::arg("mesh")) .def( "resizeConnectivity", [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt) -> void { self.resizeConnectivity(new_size, type, gt); }, py::arg("new_size"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "resizeNodes", [](MeshAccessor & self, UInt new_size) -> void { self.resizeNodes(new_size); }, py::arg("new_size")) .def("makeReady", &MeshAccessor::makeReady); } } // namespace akantu diff --git a/python/py_model.cc b/python/py_model.cc index 1373cf9dc..ab06474f6 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,102 +1,107 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_model(py::module & mod) { py::class_(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def( "hasMatrix", [](DOFManager & self, const ID & name) -> bool { return self.hasMatrix(name); }, py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual); py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_(mod, "ModelSolver", py::multiple_inheritance()) .def("getNonLinearSolver", (NonLinearSolver & (ModelSolver::*)(const ID &)) & ModelSolver::getNonLinearSolver, py::arg("solver_id") = "", py::return_value_policy::reference) .def("solveStep", [](ModelSolver & self) { self.solveStep(); }) .def("solveStep", [](ModelSolver & self, const ID & solver_id) { self.solveStep(solver_id); }); py::class_(mod, "Model", py::multiple_inheritance()) .def("setBaseName", &Model::setBaseName) .def("setDirectory", &Model::setDirectory) .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("addDumpFieldVector", &Model::addDumpFieldVector) .def("addDumpField", &Model::addDumpField) .def("setBaseNameToDumper", &Model::setBaseNameToDumper) .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper) .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper) .def("dump", py::overload_cast<>(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) .def("initNewSolver", &Model::initNewSolver) - .def("getNewSolver", [](Model & self, const std::string id, const TimeStepSolverType & time, - const NonLinearSolverType & type) { - self.getNewSolver(id, time, type); - }, py::return_value_policy::reference) - .def("setIntegrationScheme", [](Model & self, const std::string id, const std::string primal, - const IntegrationSchemeType & scheme) { - self.setIntegrationScheme(id, primal, scheme); - }) + .def( + "getNewSolver", + [](Model & self, const std::string id, + const TimeStepSolverType & time, + const NonLinearSolverType & type) { + self.getNewSolver(id, time, type); + }, + py::return_value_policy::reference) + .def("setIntegrationScheme", + [](Model & self, const std::string id, const std::string primal, + const IntegrationSchemeType & scheme) { + self.setIntegrationScheme(id, primal, scheme); + }) .def("getDOFManager", &Model::getDOFManager, py::return_value_policy::reference) .def("assembleMatrix", &Model::assembleMatrix); } } // namespace akantu diff --git a/python/py_model_couplers.cc b/python/py_model_couplers.cc new file mode 100644 index 000000000..b461d19d8 --- /dev/null +++ b/python/py_model_couplers.cc @@ -0,0 +1,90 @@ +/* -------------------------------------------------------------------------- */ +#include "py_aka_array.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +namespace py = pybind11; +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +namespace { + template + auto register_coupler_solid_contact(py::module & mod, + const std::string & name) + -> py::class_ { + return py::class_(mod, name.c_str(), + py::multiple_inheritance()) + .def(py::init, + const ModelType>(), + py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, + py::arg("id") = "coupler_solid_contact", + py::arg("dof_manager") = nullptr, + py::arg("model_type") = ModelType::_coupler_solid_contact) + .def("applyBC", + [](CouplerSolidContact_ & self, + BC::Dirichlet::DirichletFunctor & func, + const std::string & element_group) { + self.applyBC(func, element_group); + }) + .def("applyBC", + [](CouplerSolidContact_ & self, BC::Neumann::NeumannFunctor & func, + const std::string & element_group) { + self.applyBC(func, element_group); + }) + + .def("setTimeStep", &CouplerSolidContact_::setTimeStep, + py::arg("time_step"), py::arg("solver_id") = "") + .def("getContactMechanicsModel", + &CouplerSolidContact_::getContactMechanicsModel, + py::return_value_policy::reference); + } +} // namespace + +/* -------------------------------------------------------------------------- */ +void register_model_couplers(py::module & mod) { + register_coupler_solid_contact(mod, + "CouplerSolidContact") + .def( + "getSolidMechanicsModel", + [](CouplerSolidContact & self) -> decltype(auto) { + return self.getSolidMechanicsModel(); + }, + py::return_value_policy::reference) + .def( + "initFull", + [](CouplerSolidContact & self, + const AnalysisMethod & analysis_method) { + self.initFull(_analysis_method = analysis_method); + }, + py::arg("_analysis_method") = _explicit_lumped_mass); + + register_coupler_solid_contact( + mod, "CouplerSolidCohesiveContact") + .def( + "initFull", + [](CouplerSolidCohesiveContact & self, + const AnalysisMethod & analysis_method, bool is_extrinsic) { + self.initFull(_analysis_method = analysis_method, + _is_extrinsic = is_extrinsic); + }, + py::arg("_analysis_method") = _explicit_lumped_mass, + py::arg("_is_extrinsic") = false) + .def("checkCohesiveStress", + [](CouplerSolidCohesiveContact & self) { + return self.checkCohesiveStress(); + }) + .def( + "getSolidMechanicsModelCohesive", + [](CouplerSolidCohesiveContact & self) -> decltype(auto) { + return self.getSolidMechanicsModelCohesive(); + }, + py::return_value_policy::reference); +} + +} // namespace akantu diff --git a/python/py_model_couplers.hh b/python/py_model_couplers.hh new file mode 100644 index 000000000..bb517af36 --- /dev/null +++ b/python/py_model_couplers.hh @@ -0,0 +1,10 @@ +#include + +#ifndef __AKANTU_PY_MODEL_COUPLERS_HH__ +#define __AKANTU_PY_MODEL_COUPLERS_HH__ + +namespace akantu { +void register_model_couplers(pybind11::module & mod); +} // namespace akantu + +#endif // __AKANTU_PY_MODEL_COUPLERS_HH__ diff --git a/python/py_solid_mechanics_model.cc b/python/py_solid_mechanics_model.cc index b9e486598..9166c4319 100644 --- a/python/py_solid_mechanics_model.cc +++ b/python/py_solid_mechanics_model.cc @@ -1,121 +1,122 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ - void register_solid_mechanics_model(py::module & mod) { py::class_(mod, "SolidMechanicsModelOptions") .def(py::init(), py::arg("_analysis_method") = _explicit_lumped_mass); py::class_(mod, "SolidMechanicsModel", py::multiple_inheritance()) - .def(py::init(), + .def(py::init, + const ModelType>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model", + py::arg("dof_manager") = nullptr, py::arg("model_type") = ModelType::_solid_mechanics_model) .def( "initFull", [](SolidMechanicsModel & self, const SolidMechanicsModelOptions & options) { self.initFull(options); }, py::arg("option") = SolidMechanicsModelOptions()) .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_deprecated("applyDirichletBC", "Deprecated: use applyBC") .def("applyBC", [](SolidMechanicsModel & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](SolidMechanicsModel & self, BC::Neumann::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &SolidMechanicsModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getEnergy", py::overload_cast( &SolidMechanicsModel::getEnergy), py::arg("energy_id")) .def("getEnergy", py::overload_cast( &SolidMechanicsModel::getEnergy), py::arg("energy_id"), py::arg("group_id")) .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function(assembleMass) .def_function(assembleMassLumped) .def_function(getStableTimeStep) .def_function_nocopy(getExternalForce) .def_function_nocopy(getDisplacement) .def_function_nocopy(getPreviousDisplacement) .def_function_nocopy(getCurrentPosition) .def_function_nocopy(getIncrement) .def_function_nocopy(getInternalForce) .def_function_nocopy(getMass) .def_function_nocopy(getVelocity) .def_function_nocopy(getAcceleration) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getMesh) .def("getMaterial", py::overload_cast(&SolidMechanicsModel::getMaterial), py::return_value_policy::reference) .def("getMaterial", py::overload_cast( &SolidMechanicsModel::getMaterial), py::return_value_policy::reference) .def("getMaterialIndex", &SolidMechanicsModel::getMaterialIndex) // .def( // "setMaterialSelector", // [](SolidMechanicsModel & self, MaterialSelector & // material_selector) { // self.setMaterialSelector(material_selector.shared_from_this()); // }) .def("setMaterialSelector", [](SolidMechanicsModel & self, std::shared_ptr material_selector) { - std::cout << (*material_selector)(ElementNull) << std::endl; + std::cout << (*material_selector)(ElementNull) << std::endl; self.setMaterialSelector(material_selector); }) .def("getMaterialSelector", &SolidMechanicsModel::getMaterialSelector); } } // namespace akantu diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh index 6390db524..9f93dc68c 100644 --- a/src/common/aka_array.hh +++ b/src/common/aka_array.hh @@ -1,442 +1,443 @@ /** * @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