diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3eb326080..0446abb71 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,348 +1,348 @@ 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_PYTHON_INTERFACE:BOOL=TRUE -DAKANTU_EXAMPLES:BOOL=TRUE -DAKANTU_BUILD_ALL_EXAMPLES:BOOL=TRUE - -DAKANTU_TEST_EXAMPLES:BOOL=FALSE + -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_BASEDIRE: ${CI_PROJECT_DIR}/build + 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/examples/cohesive_element/cohesive_extrinsic/CMakeLists.txt b/examples/cohesive_element/cohesive_extrinsic/CMakeLists.txt index 9c563e853..b3de41926 100644 --- a/examples/cohesive_element/cohesive_extrinsic/CMakeLists.txt +++ b/examples/cohesive_element/cohesive_extrinsic/CMakeLists.txt @@ -1,39 +1,39 @@ #=============================================================================== # @file CMakeLists.txt # # @author Seyedeh Mohadeseh Taheri Mousavi # # @date creation: Mon Jan 18 2016 # # @brief configuration for extrinsic cohesive elements # # @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(cohesive_extrinsic_mesh triangle.geo 2 2) +#add_mesh(cohesive_extrinsic_mesh triangle.geo 2 2) register_example(cohesive_extrinsic SOURCES cohesive_extrinsic.cc - DEPENDS cohesive_extrinsic_mesh - FILES_TO_COPY material.dat + #DEPENDS cohesive_extrinsic_mesh + FILES_TO_COPY material.dat triangle.msh ) diff --git a/examples/cohesive_element/cohesive_extrinsic/triangle.msh b/examples/cohesive_element/cohesive_extrinsic/triangle.msh new file mode 100644 index 000000000..9acaf5dc7 --- /dev/null +++ b/examples/cohesive_element/cohesive_extrinsic/triangle.msh @@ -0,0 +1,66 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +41 +1 1 1 0 +2 -1 1 0 +3 -1 -1 0 +4 1 -1 0 +5 2.750244476601438e-12 1 0 +6 0.5000000000011042 1 0 +7 -0.4999999999986249 1 0 +8 -1 2.750244476601438e-12 0 +9 -1 0.5000000000011042 0 +10 -1 -0.4999999999986249 0 +11 -2.750244476601438e-12 -1 0 +12 -0.5000000000011042 -1 0 +13 0.4999999999986249 -1 0 +14 1 -2.750244476601438e-12 0 +15 1 -0.5000000000011042 0 +16 1 0.4999999999986249 0 +17 0 0 0 +18 0.5000000000013751 0.4999999999986249 0 +19 0.4999999999986249 -0.5000000000013751 0 +20 -0.4999999999986249 0.5000000000013751 0 +21 -0.5000000000013751 -0.4999999999986249 0 +22 1.375122238300719e-12 0.5 0 +23 0.2500000000020627 0.7499999999993124 0 +24 -0.2499999999979373 0.7500000000006876 0 +25 -0.7499999999993124 0.2500000000020627 0 +26 -0.7500000000006876 -0.2499999999979373 0 +27 -0.5 1.375122238300719e-12 0 +28 -0.2500000000020627 -0.7499999999993124 0 +29 0.2499999999979373 -0.7500000000006876 0 +30 -1.375122238300719e-12 -0.5 0 +31 0.5 -1.375122238300719e-12 0 +32 0.7499999999993124 -0.2500000000020627 0 +33 0.7500000000006876 0.2499999999979373 0 +34 0.2499999999993124 -0.2500000000006876 0 +35 0.2500000000006876 0.2499999999993124 0 +36 -0.2499999999993124 0.2500000000006876 0 +37 -0.2500000000006876 -0.2499999999993124 0 +38 0.7500000000006876 0.7499999999993124 0 +39 -0.7499999999993124 0.7500000000006876 0 +40 -0.7500000000006876 -0.7499999999993124 0 +41 0.7499999999993124 -0.7500000000006876 0 +$EndNodes +$Elements +16 +1 9 2 7 6 20 18 5 22 23 24 +2 9 2 7 6 20 8 21 25 26 27 +3 9 2 7 6 21 11 19 28 29 30 +4 9 2 7 6 18 19 14 31 32 33 +5 9 2 7 6 17 19 18 34 31 35 +6 9 2 7 6 17 18 20 35 22 36 +7 9 2 7 6 17 21 19 37 30 34 +8 9 2 7 6 17 20 21 36 27 37 +9 9 2 7 6 1 18 14 38 33 16 +10 9 2 7 6 2 20 5 39 24 7 +11 9 2 7 6 3 21 8 40 26 10 +12 9 2 7 6 4 19 11 41 29 13 +13 9 2 7 6 1 5 18 6 23 38 +14 9 2 7 6 2 8 20 9 25 39 +15 9 2 7 6 3 11 21 12 28 40 +16 9 2 7 6 4 14 19 15 32 41 +$EndElements diff --git a/examples/cohesive_element/cohesive_intrinsic/CMakeLists.txt b/examples/cohesive_element/cohesive_intrinsic/CMakeLists.txt index afbdccbd5..a149423a4 100644 --- a/examples/cohesive_element/cohesive_intrinsic/CMakeLists.txt +++ b/examples/cohesive_element/cohesive_intrinsic/CMakeLists.txt @@ -1,38 +1,37 @@ #=============================================================================== # @file CMakeLists.txt # # @author Seyedeh Mohadeseh Taheri Mousavi # # @date creation: Mon Jan 18 2016 # # @brief Intrinsic cohesive element configuration # # @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( cohesive_intrinsic_mesh triangle.geo 2 2) +#add_mesh( cohesive_intrinsic_mesh triangle.geo 2 2) register_example(cohesive_intrinsic SOURCES cohesive_intrinsic.cc - DEPENDS cohesive_intrinsic_mesh - FILES_TO_COPY material.dat) - + #DEPENDS cohesive_intrinsic_mesh + FILES_TO_COPY material.dat triangle.msh) diff --git a/examples/cohesive_element/cohesive_intrinsic/triangle.msh b/examples/cohesive_element/cohesive_intrinsic/triangle.msh new file mode 100644 index 000000000..9acaf5dc7 --- /dev/null +++ b/examples/cohesive_element/cohesive_intrinsic/triangle.msh @@ -0,0 +1,66 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +41 +1 1 1 0 +2 -1 1 0 +3 -1 -1 0 +4 1 -1 0 +5 2.750244476601438e-12 1 0 +6 0.5000000000011042 1 0 +7 -0.4999999999986249 1 0 +8 -1 2.750244476601438e-12 0 +9 -1 0.5000000000011042 0 +10 -1 -0.4999999999986249 0 +11 -2.750244476601438e-12 -1 0 +12 -0.5000000000011042 -1 0 +13 0.4999999999986249 -1 0 +14 1 -2.750244476601438e-12 0 +15 1 -0.5000000000011042 0 +16 1 0.4999999999986249 0 +17 0 0 0 +18 0.5000000000013751 0.4999999999986249 0 +19 0.4999999999986249 -0.5000000000013751 0 +20 -0.4999999999986249 0.5000000000013751 0 +21 -0.5000000000013751 -0.4999999999986249 0 +22 1.375122238300719e-12 0.5 0 +23 0.2500000000020627 0.7499999999993124 0 +24 -0.2499999999979373 0.7500000000006876 0 +25 -0.7499999999993124 0.2500000000020627 0 +26 -0.7500000000006876 -0.2499999999979373 0 +27 -0.5 1.375122238300719e-12 0 +28 -0.2500000000020627 -0.7499999999993124 0 +29 0.2499999999979373 -0.7500000000006876 0 +30 -1.375122238300719e-12 -0.5 0 +31 0.5 -1.375122238300719e-12 0 +32 0.7499999999993124 -0.2500000000020627 0 +33 0.7500000000006876 0.2499999999979373 0 +34 0.2499999999993124 -0.2500000000006876 0 +35 0.2500000000006876 0.2499999999993124 0 +36 -0.2499999999993124 0.2500000000006876 0 +37 -0.2500000000006876 -0.2499999999993124 0 +38 0.7500000000006876 0.7499999999993124 0 +39 -0.7499999999993124 0.7500000000006876 0 +40 -0.7500000000006876 -0.7499999999993124 0 +41 0.7499999999993124 -0.7500000000006876 0 +$EndNodes +$Elements +16 +1 9 2 7 6 20 18 5 22 23 24 +2 9 2 7 6 20 8 21 25 26 27 +3 9 2 7 6 21 11 19 28 29 30 +4 9 2 7 6 18 19 14 31 32 33 +5 9 2 7 6 17 19 18 34 31 35 +6 9 2 7 6 17 18 20 35 22 36 +7 9 2 7 6 17 21 19 37 30 34 +8 9 2 7 6 17 20 21 36 27 37 +9 9 2 7 6 1 18 14 38 33 16 +10 9 2 7 6 2 20 5 39 24 7 +11 9 2 7 6 3 21 8 40 26 10 +12 9 2 7 6 4 19 11 41 29 13 +13 9 2 7 6 1 5 18 6 23 38 +14 9 2 7 6 2 8 20 9 25 39 +15 9 2 7 6 3 11 21 12 28 40 +16 9 2 7 6 4 14 19 15 32 41 +$EndElements diff --git a/examples/python/cohesive/plate.py b/examples/python/cohesive/plate.py index d8090a1cc..214f08c52 100644 --- a/examples/python/cohesive/plate.py +++ b/examples/python/cohesive/plate.py @@ -1,84 +1,84 @@ #!/usr/bin/env python3 import akantu as aka import numpy as np def solve(material_file, mesh_file, traction): aka.parseInput(material_file) spatial_dimension = 2 # ------------------------------------------------------------------------- # Initialization # ------------------------------------------------------------------------- mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) model = aka.SolidMechanicsModelCohesive(mesh) model.initFull(_analysis_method=aka._static, _is_extrinsic=True) model.initNewSolver(aka._explicit_lumped_mass) model.setBaseName('plate') model.addDumpFieldVector('displacement') model.addDumpFieldVector('external_force') model.addDumpField('strain') model.addDumpField('stress') model.addDumpField('blocked_dofs') model.setBaseNameToDumper('cohesive elements', 'cohesive') model.addDumpFieldVectorToDumper('cohesive elements', 'displacement') model.addDumpFieldToDumper('cohesive elements', 'damage') model.addDumpFieldVectorToDumper('cohesive elements', 'traction') model.addDumpFieldVectorToDumper('cohesive elements', 'opening') # ------------------------------------------------------------------------- # Boundary conditions # ------------------------------------------------------------------------- model.applyBC(aka.FixedValue(0.0, aka._x), 'XBlocked') model.applyBC(aka.FixedValue(0.0, aka._y), 'YBlocked') trac = np.zeros(spatial_dimension) trac[int(aka._y)] = traction print('Solve for traction ', traction) model.getExternalForce()[:] = 0 model.applyBC(aka.FromTraction(trac), 'Traction') solver = model.getNonLinearSolver('static') solver.set('max_iterations', 100) solver.set('threshold', 1e-10) - solver.set('convergence_type', aka._scc_residual) + solver.set("convergence_type", aka.SolveConvergenceCriteria.residual) model.solveStep('static') model.dump() model.dump('cohesive elements') model.setTimeStep(model.getStableTimeStep()*0.1) maxsteps = 100 for i in range(0, maxsteps): print('{0}/{1}'.format(i, maxsteps)) model.checkCohesiveStress() model.solveStep('explicit_lumped') if i % 10 == 0: model.dump() model.dump('cohesive elements') # ----------------------------------------------------------------------------- # main # ----------------------------------------------------------------------------- def main(): mesh_file = 'plate.msh' material_file = 'material.dat' traction = .095 solve(material_file, mesh_file, traction) # ----------------------------------------------------------------------------- if __name__ == '__main__': main() diff --git a/examples/python/custom-material/bi-material.py b/examples/python/custom-material/bi-material.py index d58c12058..823f7dee8 100644 --- a/examples/python/custom-material/bi-material.py +++ b/examples/python/custom-material/bi-material.py @@ -1,178 +1,180 @@ import akantu as aka import numpy as np # ------------------------------------------------------------------------------ class LocalElastic(aka.Material): def __init__(self, model, _id): super().__init__(model, _id) super().registerParamReal('E', aka._pat_readable | aka._pat_parsable, 'Youngs modulus') super().registerParamReal('nu', aka._pat_readable | aka._pat_parsable, 'Poisson ratio') # change it to have the initialize wrapped - super().registerInternal('factor', 1) - super().registerInternal('quad_coordinates', 2) + super().registerInternalReal('factor', 1) + super().registerInternalReal('quad_coordinates', 2) def initMaterial(self): nu = self.getReal('nu') E = self.getReal('E') self.mu = E / (2 * (1 + nu)) self.lame_lambda = nu * E / ( (1. + nu) * (1. - 2. * nu)) # Second Lame coefficient (shear modulus) self.lame_mu = E / (2. * (1. + nu)) super().initMaterial() - quad_coords = self.internals["quad_coordinates"] - factor = self.internals["factor"] + quad_coords = self.getInternalReal("quad_coordinates") + factor = self.getInternalReal("factor") model = self.getModel() model.getFEEngine().computeIntegrationPointsCoordinates( - quad_coords, self.element_filter) + quad_coords, self.getElementFilter()) for elem_type in factor.elementTypes(): factor = factor(elem_type) coords = quad_coords(elem_type) factor[:] = 1. factor[coords[:, 1] < 0.5] = .5 # declares all the parameters that are needed def getPushWaveSpeed(self, params): return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho) # compute small deformation tensor @staticmethod def computeEpsilon(grad_u): return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) # constitutive law def computeStress(self, el_type, ghost_type): grad_u = self.getGradU(el_type, ghost_type) sigma = self.getStress(el_type, ghost_type) n_quads = grad_u.shape[0] grad_u = grad_u.reshape((n_quads, 2, 2)) - factor = self.internals['factor'](el_type, ghost_type).reshape(n_quads) + factor = self.getInternalReal('factor')( + el_type, ghost_type).reshape(n_quads) epsilon = self.computeEpsilon(grad_u) sigma = sigma.reshape((n_quads, 2, 2)) trace = np.einsum('aii->a', grad_u) sigma[:, :, :] = ( np.einsum('a,ij->aij', trace, self.lame_lambda * np.eye(2)) + 2. * self.lame_mu * epsilon) sigma[:, :, :] = np.einsum('aij, a->aij', sigma, factor) # constitutive law tangent modulii def computeTangentModuli(self, el_type, tangent_matrix, ghost_type): n_quads = tangent_matrix.shape[0] tangent = tangent_matrix.reshape(n_quads, 3, 3) - factor = self.internals['factor'](el_type, ghost_type).reshape(n_quads) + factor = self.getInternalReal('factor')( + el_type, ghost_type).reshape(n_quads) Miiii = self.lame_lambda + 2 * self.lame_mu Miijj = self.lame_lambda Mijij = self.lame_mu tangent[:, 0, 0] = Miiii tangent[:, 1, 1] = Miiii tangent[:, 0, 1] = Miijj tangent[:, 1, 0] = Miijj tangent[:, 2, 2] = Mijij tangent[:, :, :] = np.einsum('aij, a->aij', tangent, factor) # computes the energy density def computePotentialEnergy(self, el_type): sigma = self.getStress(el_type) grad_u = self.getGradU(el_type) nquads = sigma.shape[0] stress = sigma.reshape(nquads, 2, 2) grad_u = grad_u.reshape((nquads, 2, 2)) epsilon = self.computeEpsilon(grad_u) energy_density = self.getPotentialEnergy(el_type) energy_density[:, 0] = 0.5 * np.einsum('aij,aij->a', stress, epsilon) # ------------------------------------------------------------------------------ # applies manually the boundary conditions def applyBC(model): nbNodes = model.getMesh().getNbNodes() position = model.getMesh().getNodes() displacement = model.getDisplacement() blocked_dofs = model.getBlockedDOFs() width = 1. height = 1. epsilon = 1e-8 for node in range(0, nbNodes): if((np.abs(position[node, 0]) < epsilon) or # left side (np.abs(position[node, 0] - width) < epsilon)): # right side blocked_dofs[node, 0] = True displacement[node, 0] = 0 * position[node, 0] + 0. if(np.abs(position[node, 1]) < epsilon): # lower side blocked_dofs[node, 1] = True displacement[node, 1] = - 1. if(np.abs(position[node, 1] - height) < epsilon): # upper side blocked_dofs[node, 1] = True displacement[node, 1] = 1. # register the material to the material factory def allocator(dim, option, model, id): return LocalElastic(model, id) mat_factory = aka.MaterialFactory.getInstance() mat_factory.registerAllocator("local_elastic", allocator) # main parameters spatial_dimension = 2 mesh_file = 'square.msh' # read mesh mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) # parse input file aka.parseInput('material.dat') # init the SolidMechanicsModel model = aka.SolidMechanicsModel(mesh) model.initFull(_analysis_method=aka._static) # configure the solver solver = model.getNonLinearSolver() solver.set("max_iterations", 2) solver.set("threshold", 1e-3) solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) # prepare the dumper model.setBaseName("bimaterial") model.addDumpFieldVector("displacement") model.addDumpFieldVector("internal_force") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") # model.addDumpField("factor") model.addDumpField("blocked_dofs") # Boundary conditions applyBC(model) # solve the problem model.solveStep() # dump paraview files model.dump() epot = model.getEnergy('potential') print('Potential energy: ' + str(epot)) diff --git a/examples/python/eigen_modes/eigen_modes.py b/examples/python/eigen_modes/eigen_modes.py index acee26aa6..5f7b72047 100644 --- a/examples/python/eigen_modes/eigen_modes.py +++ b/examples/python/eigen_modes/eigen_modes.py @@ -1,257 +1,264 @@ -import sys +#!/usr/bin/env python + import subprocess import argparse import akantu as aka import numpy as np -from image_saver import ImageSaver -import matplotlib.pyplot as plt from scipy.sparse.linalg import eigsh from scipy.sparse import csr_matrix +try: + import matplotlib.pyplot as plt + from image_saver import ImageSaver + has_matplotlib = True +except ImportError: + has_matplotlib = False # ----------------------------------------------------------------------------- # parser # ----------------------------------------------------------------------------- parser = argparse.ArgumentParser(description='Eigen mode exo') parser.add_argument('-m', '--mode_number', type=int, help='precise the mode to study', default=2) parser.add_argument('-wL', '--wave_width', type=float, help='precise the width of the wave for ' 'the initial displacement', default=5) parser.add_argument('-L', '--Lbar', type=float, help='precise the length of the bar', default=10) parser.add_argument('-t', '--time_step', type=float, help='precise the timestep', default=None) parser.add_argument('-n', '--max_steps', type=int, help='precise the number of timesteps', default=500) parser.add_argument('-mh', '--mesh_h', type=float, help='characteristic mesh size', default=.2) parser.add_argument('-p', '--plot', action='store_true', help='plot the results') args = parser.parse_args() mode = args.mode_number wave_width = args.wave_width time_step = args.time_step max_steps = args.max_steps mesh_h = args.mesh_h Lbar = args.Lbar plot = args.plot # ----------------------------------------------------------------------------- # Mesh Generation # ----------------------------------------------------------------------------- geo_content = """ // Mesh size h = {0}; """.format(mesh_h) geo_content += """ h1 = h; h2 = h; // Dimensions of the bar Lx = 10; Ly = 1; // ------------------------------------------ // Geometry // ------------------------------------------ Point(101) = { 0.0, -Ly/2, 0.0, h1}; Point(102) = { Lx, -Ly/2, 0.0, h2}; Point(103) = { Lx, 0., 0.0, h2}; Point(104) = { Lx, Ly/2., 0.0, h2}; Point(105) = { 0.0, Ly/2., 0.0, h1}; Point(106) = { 0.0, 0., 0.0, h1}; Line(101) = {101,102}; Line(102) = {102,103}; Line(103) = {103,104}; Line(104) = {104,105}; Line(105) = {105,106}; Line(106) = {106,101}; Line(107) = {106,103}; Line Loop(108) = {101, 102, -107, 106}; Plane Surface(109) = {108}; Line Loop(110) = {103, 104, 105, 107}; Plane Surface(111) = {110}; Physical Surface(112) = {109, 111}; Transfinite Surface "*"; Recombine Surface "*"; Physical Surface(113) = {111, 109}; Physical Line("XBlocked") = {103, 102}; Physical Line("ImposedVelocity") = {105, 106}; Physical Line("YBlocked") = {104, 101}; """ mesh_file = 'bar' with open(mesh_file + '.geo', 'w') as f: f.write(geo_content) subprocess.call(['gmsh', '-2', mesh_file + '.geo']) mesh_file = mesh_file + '.msh' # ----------------------------------------------------------------------------- # Initialization # ----------------------------------------------------------------------------- spatial_dimension = 2 aka.parseInput('material.dat') mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) model = aka.SolidMechanicsModel(mesh) model.initFull(aka._implicit_dynamic) model.setBaseName("waves-{0}".format(mode)) model.addDumpFieldVector("displacement") model.addDumpFieldVector("acceleration") model.addDumpFieldVector("velocity") model.addDumpField("blocked_dofs") # ----------------------------------------------------------------------------- # Boundary conditions # ----------------------------------------------------------------------------- internal_force = model.getInternalForce() displacement = model.getDisplacement() acceleration = model.getAcceleration() velocity = model.getVelocity() blocked_dofs = model.getBlockedDOFs() nbNodes = mesh.getNbNodes() position = mesh.getNodes() model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked") model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked") # ------------------------------------------------------------------------ # timestep value computation # ------------------------------------------------------------------------ time_factor = 0.8 stable_time_step = model.getStableTimeStep() * time_factor if time_step: print("Required Time Step = {0}".format(time_step)) if stable_time_step * time_factor < time_step: print("Stable Time Step = {0}".format(stable_time_step)) raise RuntimeError("required time_step too large") print("Required Time Step = {0}".format(time_step)) else: print("Stable Time Step = {0}".format(stable_time_step)) time_step = stable_time_step * time_factor model.setTimeStep(time_step) -disp_sav = ImageSaver(mesh, displacement, 0, Lbar) -velo_sav = ImageSaver(mesh, velocity, 0, Lbar) - # ------------------------------------------------------------------------ # compute the eigen modes # ------------------------------------------------------------------------ model.assembleStiffnessMatrix() model.assembleMass() stiff = model.getDOFManager().getMatrix('K') stiff = aka.AkantuSparseMatrix(stiff).toarray() mass = model.getDOFManager().getMatrix('M') mass = aka.AkantuSparseMatrix(mass).toarray() # select the non blocked DOFs by index in the mask mask = np.equal(blocked_dofs.flatten(), False) Mass_star = mass[mask, :] Mass_star = csr_matrix(Mass_star[:, mask].copy()) K_star = stiff[mask, :] K_star = csr_matrix(K_star[:, mask].copy()) print('getting the eigen values') vals, vects = eigsh(K_star, M=Mass_star, which='SM', k=20) # ----------------------------------------------------------------------------- # import the initial conditions in displacement # ----------------------------------------------------------------------------- displacement.reshape(nbNodes*2)[mask] = vects[:, mode] with open('modes.txt', 'a') as f: f.write('{0} {1}\n'.format(mode, vals[mode])) model.dump() # ----------------------------------------------------------------------------- # prepare the storage of the dynamical evolution # ----------------------------------------------------------------------------- e_p = np.zeros(max_steps + 1) e_k = np.zeros(max_steps + 1) e_t = np.zeros(max_steps + 1) time = np.zeros(max_steps + 1) norm = np.zeros(max_steps + 1) epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') e_p[0] = epot e_k[0] = ekin e_t[0] = epot + ekin time[0] = 0 +if has_matplotlib: + disp_sav = ImageSaver(mesh, displacement, 0, Lbar) + velo_sav = ImageSaver(mesh, velocity, 0, Lbar) + + # ----------------------------------------------------------------------------- # loop for evolution of motion dynamics # ----------------------------------------------------------------------------- for step in range(1, max_steps + 1): model.solveStep() # outputs epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') print(step, '/', max_steps, epot, ekin, epot + ekin) e_p[step] = epot e_k[step] = ekin e_t[step] = epot + ekin time[step] = (step + 1) * time_step - disp_sav.storeStep() - velo_sav.storeStep() + if has_matplotlib: + disp_sav.storeStep() + velo_sav.storeStep() + if step % 10 == 0: model.dump() -if not plot: - sys.exit(0) - -# ----------------------------------------------------------------------------- -# plot figures for global evolution -# ----------------------------------------------------------------------------- -# energy norms -plt.figure(1) -plt.plot(time, e_t, 'r', time, e_p, 'b', time, e_k, 'g') - -# space-time diagram for diplacements -plt.figure(2) -plt.imshow(disp_sav.getImage(), extent=(0, Lbar, max_steps * time_step, 0)) -plt.xlabel("Space ") -plt.ylabel("Time ") - -# space-time diagram for velocities -plt.figure(3) -plt.imshow(velo_sav.getImage(), extent=(0, Lbar, max_steps * time_step, 0)) -plt.xlabel("Velocity") -plt.ylabel("Time") -plt.show() +if plot and has_matplotlib: + # -------------------------------------------------------------------------- + # plot figures for global evolution + # -------------------------------------------------------------------------- + # energy norms + plt.figure(1) + plt.plot(time, e_t, 'r', time, e_p, 'b', time, e_k, 'g') + + # space-time diagram for diplacements + plt.figure(2) + plt.imshow(disp_sav.getImage(), extent=(0, Lbar, max_steps * time_step, 0)) + plt.xlabel("Space ") + plt.ylabel("Time ") + + # space-time diagram for velocities + plt.figure(3) + plt.imshow(velo_sav.getImage(), extent=(0, Lbar, max_steps * time_step, 0)) + plt.xlabel("Velocity") + plt.ylabel("Time") + plt.show() diff --git a/examples/python/structural_mechanics/CMakeLists.txt b/examples/python/structural_mechanics/CMakeLists.txt index f8a6b2751..7496e58f4 100644 --- a/examples/python/structural_mechanics/CMakeLists.txt +++ b/examples/python/structural_mechanics/CMakeLists.txt @@ -1,14 +1,14 @@ register_example(structural_mechanics.py SCRIPT structural_mechanics.py PYTHON ) register_example(structural_mechanics_softening.py - SCRIPT structural_mechanics.py + SCRIPT structural_mechanics_softening.py PYTHON ) register_example(structural_mechanics_dynamics.py - SCRIPT structural_mechanics.py + SCRIPT structural_mechanics_dynamics.py PYTHON ) diff --git a/examples/python/structural_mechanics/structural_mechanics.py b/examples/python/structural_mechanics/structural_mechanics.py index 8599a50a0..d0706beda 100644 --- a/examples/python/structural_mechanics/structural_mechanics.py +++ b/examples/python/structural_mechanics/structural_mechanics.py @@ -1,135 +1,141 @@ #!/usr/bin/env python # coding: utf-8 # # Test of Structural Mechanics # In this example a beam, consisting of two elements, three nodes, is created. # The left most node is fixed and a force is applied at the right most node. import akantu as aka import numpy -import matplotlib.pyplot as plt import numpy as np +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except ImportError: + has_matplotlib = False # ### Creating the Mesh # Create a mesh for the two dimensional case beam = aka.Mesh(2) # We now create the connectivity array for the beam. -beam.addConnectivity(aka._bernoulli_beam_2) +beam.addConnectivityType(aka._bernoulli_beam_2) # We need a `MeshAccessor` in order to change the size of the mesh entities. beamAcc = aka.MeshAccessor(beam) -# Now we create the array to store the nodes and the connectivities and give them their size. +# Now we create the array to store the nodes and the connectivities and give +# them their size. beamAcc.resizeConnectivity(2, aka._bernoulli_beam_2) beamAcc.resizeNodes(3) Nodes = beam.getNodes() Nodes[0, :] = [0., 0.] Nodes[1, :] = [1., 0.] Nodes[2, :] = [2., 0.] # #### Setting the Connections Conn = beam.getConnectivity(aka._bernoulli_beam_2) Conn[0, :] = [0, 1] Conn[1, :] = [1, 2] # #### Ready # We have to make the mesh ready. beamAcc.makeReady() # ### Creating the Model model = aka.StructuralMechanicsModel(beam) # #### Setting up the Modell # ##### Creating and Inserting the Materials mat1 = aka.StructuralMaterial() mat1.E = 1e9 mat1.rho = 1. -mat1.I = 1. +mat1.I = 1. # noqa: E741 mat1.Iz = 1. mat1.Iy = 1. mat1.A = 1. mat1.GJ = 1. model.addMaterial(mat1) mat2 = aka.StructuralMaterial() mat2.E = 1e9 mat2.rho = 1. -mat2.I = 1. +mat2.I = 1. # noqa: E741 mat2.Iz = 1. mat2.Iy = 1. mat2.A = 1. mat2.GJ = 1. model.addMaterial(mat2) # ##### Initializing the Model model.initFull(aka._implicit_dynamic) # ##### Assigning the Materials materials = model.getElementMaterial(aka._bernoulli_beam_2) materials[0][0] = 0 materials[1][0] = 1 # ##### Setting Boundaries # Neumann # Apply a force of `10` at the last (right most) node. forces = model.getExternalForce() forces[:] = 0 forces[2, 0] = 100. # Dirichlets # Block all dofs of the first node, since it is fixed. # All other nodes have no restrictions boundary = model.getBlockedDOFs() boundary[0, :] = True boundary[1, :] = False boundary[2, :] = False # ### Solving the System # Set up the system deltaT = 1e-10 model.setTimeStep(deltaT) solver = model.getNonLinearSolver() solver.set("max_iterations", 100) solver.set("threshold", 1e-8) solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) # Perform N time steps. # At each step records the displacement of all three nodes in x direction. N = 1000000 disp1 = np.zeros(N) disp2 = np.zeros(N) disp0 = np.zeros(N) times = np.zeros(N) for i in range(N): model.solveStep() disp = model.getDisplacement() disp0[i] = disp[0, 0] disp1[i] = disp[1, 0] disp2[i] = disp[2, 0] times[i] = deltaT * i disps = [disp0, disp1, disp2] maxMin = [-1.0, 1.0] for d in disps: maxMin[0] = max(np.max(d), maxMin[0]) maxMin[1] = min(np.min(d), maxMin[1]) -plt.plot(disp1, times, color='g', label = "middle node") -plt.plot(disp2, times, color='b', label = "right node") +if has_matplotlib: + plt.plot(disp1, times, color='g', label="middle node") + plt.plot(disp2, times, color='b', label="right node") -plt.title("Displacement in $x$ of the nodes") -plt.ylabel("Time [S]") -plt.xlabel("displacement [m]") + plt.title("Displacement in $x$ of the nodes") + plt.ylabel("Time [S]") + plt.xlabel("displacement [m]") -plt.xlim((maxMin[1] * 1.3, maxMin[0] * 1.1)) + plt.xlim((maxMin[1] * 1.3, maxMin[0] * 1.1)) -plt.legend() + plt.legend() -plt.show() + plt.show() diff --git a/examples/python/structural_mechanics/structural_mechanics_dynamics.py b/examples/python/structural_mechanics/structural_mechanics_dynamics.py index b3d7cf850..4a9745a3a 100644 --- a/examples/python/structural_mechanics/structural_mechanics_dynamics.py +++ b/examples/python/structural_mechanics/structural_mechanics_dynamics.py @@ -1,166 +1,175 @@ #!/usr/bin/env python # coding: utf-8 import akantu as aka import numpy -import matplotlib.pyplot as plt import numpy as np +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except ImportError: + has_matplotlib = False # ### Creating the Mesh # Create a mesh for the two dimensional case el_type = aka._bernoulli_beam_2 beam = aka.Mesh(2) # We now create the connectivity array for the beam. beam.addConnectivityType(el_type) # We need a `MeshAccessor` in order to change the size of the mesh entities. beamAcc = aka.MeshAccessor(beam) -# Now we create the array to store the nodes and the connectivities and give them their size. +# Now we create the array to store the nodes and the connectivities and give +# them their size. nb_elem = 40 L = 2 beamAcc.resizeConnectivity(nb_elem, el_type) beamAcc.resizeNodes(nb_elem + 1) # #### Setting the Nodes Nodes = beam.getNodes() length = L / nb_elem Nodes[:, :] = 0. Nodes[:, 0] = np.arange(nb_elem+1) * length # #### Setting the Connections Conn = beam.getConnectivity(el_type) for e in range(nb_elem): Conn[e, :] = [e, e + 1] # #### Ready # We have to make the mesh ready. beamAcc.makeReady() # ### Creating the Model model = aka.StructuralMechanicsModel(beam) if el_type == aka._bernoulli_beam_3: normal = beam.getDataReal("extra_normal", el_type) for e in range(nb_elem): normal[e, :] = [0, 0, 1] # #### Setting up the Modell # ##### Creating and Inserting the Materials mat1 = aka.StructuralMaterial() mat1.E = 1e9 mat1.rho = 10. -mat1.I = 1. +mat1.I = 1. # noqa: E741 mat1.Iz = 1. mat1.Iy = 1. mat1.A = 1. mat1.GJ = 1. model.addMaterial(mat1, 'mat1') # ##### Initializing the Model model.initFull(aka.AnalysisMethod._implicit_dynamic) # ##### Assigning the Materials materials = model.getElementMaterial(el_type) materials[:, :] = 0 # ##### Setting Boundaries # Neumann F = 1e4 no_print = int(nb_elem / 2) # Apply a force of `10` at the last (right most) node. forces = model.getExternalForce() forces[:, :] = 0 forces[no_print, 1] = F # Dirichlets # Block all dofs of the first node, since it is fixed. # All other nodes have no restrictions boundary = model.getBlockedDOFs() boundary[:, :] = False boundary[0, 0] = True boundary[0, 1] = True if el_type == aka._bernoulli_beam_3: boundary[0, 2] = True boundary[nb_elem, 1] = True # ### Solving the System # Set up the system deltaT = 1e-6 model.setTimeStep(deltaT) solver = model.getNonLinearSolver() solver.set("max_iterations", 100) solver.set("threshold", 1e-8) solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) model.assembleMatrix("M") M_ = model.getDOFManager().getMatrix("M") M = aka.AkantuSparseMatrix(M_) model.assembleMatrix("K") K_ = model.getDOFManager().getMatrix("K") K = aka.AkantuSparseMatrix(K_) C_ = model.getDOFManager().getMatrix("C") C_.add(M_, 0.00001) C_.add(K_, 0.00001) -def analytical_solution(time, L, rho, E, A, I, F): + +def analytical_solution(time, L, rho, E, A, I, F): # noqa: E741 omega = np.pi**2 / L**2 * np.sqrt(E * I / rho) sum = 0. N = 110 for n in range(1, N, 2): sum += (1. - np.cos(n * n * omega * time)) / n**4 return 2. * F * L**3 / np.pi**4 / E / I * sum + # Perform N time steps. # At each step records the displacement of all three nodes in x direction. N = 900 mat1 = model.getMaterial('mat1') disp = model.getDisplacement() velo = model.getVelocity() disp[:, :] = 0. displs = np.zeros(N) ekin = np.zeros(N) epot = np.zeros(N) ework = np.zeros(N) _ework = 0. for i in range(1, N): model.solveStep() displs[i] = disp[no_print, 1] _ework += F * velo[no_print, 1] * deltaT ekin[i] = model.getEnergy("kinetic") epot[i] = model.getEnergy("potential") ework[i] = _ework def sol(x): return analytical_solution(x, L, mat1.rho, mat1.E, mat1.A, mat1.I, F) -times = np.arange(N) * deltaT -plt.plot(times, sol(times)) -plt.plot(times, displs) -plt.plot(times, displs - sol(times)) +if has_matplotlib: + times = np.arange(N) * deltaT + plt.plot(times, sol(times)) + plt.plot(times, displs) + plt.plot(times, displs - sol(times)) -# What I do not fully understand is why the middle node first go backwards until it goes forward. -# I could imagine that there is some vibration, because everything is in rest. -np.max(displs - sol(times)) -plt.plot(times, ekin+epot) -plt.plot(times, ework) + # What I do not fully understand is why the middle node first go backwards + # until it goes forward. I could imagine that there is some vibration, + # because everything is in rest. + np.max(displs - sol(times)) + plt.plot(times, ekin+epot) + plt.plot(times, ework) diff --git a/examples/python/structural_mechanics/structural_mechanics_softening.py b/examples/python/structural_mechanics/structural_mechanics_softening.py index b196c8a71..602852693 100644 --- a/examples/python/structural_mechanics/structural_mechanics_softening.py +++ b/examples/python/structural_mechanics/structural_mechanics_softening.py @@ -1,188 +1,189 @@ #!/usr/bin/env python # coding: utf-8 # # Test of Structural Mechanics # In this test there is a beam consisting of three parts, all have the same materials. # The left most node is fixed. # On the right most node a force is applied in x direction. # # After a certain time, the material of the middle _element_ is waekened, lower Young's modulus. # In each step the modulus is lowered by a coinstant factor. import akantu as aka import numpy -import matplotlib.pyplot as plt import numpy as np +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except ImportError: + has_matplotlib = False # ### Creating the Mesh # Create a mesh for the two dimensional case beam = aka.Mesh(2) # We now create the connectivity array for the beam. beam.addConnectivityType(aka._bernoulli_beam_2) # We need a `MeshAccessor` in order to change the size of the mesh entities. beamAcc = aka.MeshAccessor(beam) # Now we create the array to store the nodes and the connectivities and give them their size. beamAcc.resizeConnectivity(3, aka._bernoulli_beam_2) beamAcc.resizeNodes(4) # #### Setting the Nodes Nodes = beam.getNodes() Nodes[0, :] = [0., 0.] Nodes[1, :] = [1., 0.] Nodes[2, :] = [2., 0.] Nodes[3, :] = [3., 0.] -# #### Setting the Connections +# Setting the Connections Conn = beam.getConnectivity(aka._bernoulli_beam_2) Conn[0, :] = [0, 1] Conn[1, :] = [1, 2] Conn[2, :] = [2, 3] -#### Ready +# Ready # We have to make the mesh ready. beamAcc.makeReady() -# ### Creating the Model +# Creating the Model model = aka.StructuralMechanicsModel(beam) -# #### Setting up the Modell -# ##### Creating and Inserting the Materials +# Setting up the Modell +# Creating and Inserting the Materials mat1 = aka.StructuralMaterial() mat1.E = 1e9 mat1.rho = 1. -mat1.I = 1. -mat1.Iz = 1. -mat1.Iy = 1. +mat1.I = 1. # noqa: E741 mat1.A = 1. mat1.GJ = 1. mat1ID = model.addMaterial(mat1, 'mat1') mat2 = aka.StructuralMaterial() mat2.E = 1e9 mat2.rho = 1. -mat2.I = 1. -mat2.Iz = 1. -mat2.Iy = 1. +mat2.I = 1. # noqa: E741 mat2.A = 1. mat2.GJ = 1. mat2ID = model.addMaterial(mat2, 'mat2') mat3 = aka.StructuralMaterial() mat3.E = mat2.E / 100000 mat3.rho = 1. -mat3.I = 1. -mat3.Iz = 1. -mat3.Iy = 1. +mat3.I = 1. # noqa: E741 mat3.A = mat2.A / 100 mat3.GJ = 1. mat3ID = model.addMaterial(mat3, 'mat3') # ##### Initializing the Model model.initFull(aka.AnalysisMethod._implicit_dynamic) # ##### Assigning the Materials materials = model.getElementMaterial(aka._bernoulli_beam_2) materials[0][0] = mat1ID materials[1][0] = mat2ID materials[2][0] = mat1ID # ##### Setting Boundaries # Neumann # Apply a force of `10` at the last (right most) node. forces = model.getExternalForce() forces[:] = 0 forces[2, 0] = 100. # Dirichlets # Block all dofs of the first node, since it is fixed. # All other nodes have no restrictions boundary = model.getBlockedDOFs() boundary[0, :] = True boundary[1, :] = False boundary[2, :] = False boundary[3, :] = False # ### Solving the System # Set up the system deltaT = 1e-9 model.setTimeStep(deltaT) solver = model.getNonLinearSolver() solver.set("max_iterations", 100) solver.set("threshold", 1e-8) solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) # Perform N time steps. # At each step records the displacement of all three nodes in x direction. N = 10000 * 60 disp0 = np.zeros(N) disp1 = np.zeros(N) disp2 = np.zeros(N) disp3 = np.zeros(N) times = np.zeros(N) switchT = None switchEnd = None softDuration = 1000 SoftStart = (N // 2) - softDuration // 2 SoftEnd = SoftStart + softDuration if(softDuration > 0): - softFactor = (model.getMaterial('mat3ID').E / model.getMaterial('mat2ID').E) ** (1.0 / softDuration) + softFactor = (model.getMaterial('mat3').E + / model.getMaterial('mat2').E) ** (1.0 / softDuration) mat2 = model.getMaterial('mat2') for i in range(N): times[i] = deltaT * i if((SoftStart <= i <= SoftEnd) and (softDuration > 0)): if switchT is None: switchT = times[i] elif(i == SoftEnd): switchEnd = times[i] # mat2.E *= softFactor # model.solveStep() disp = model.getDisplacement() disp0[i] = disp[0, 0] disp1[i] = disp[1, 0] disp2[i] = disp[2, 0] disp3[i] = disp[3, 0] disps = [disp0, disp1, disp2, disp3] maxMin = [-1.0, 1.0] for d in disps: maxMin[0] = max(np.max(d), maxMin[0]) maxMin[1] = min(np.min(d), maxMin[1]) -#plt.plot(disp0, times, color='k', label = "left node (fix)") -plt.plot(disp1, times, color='g', label = "middle, left node") -plt.plot(disp2, times, color='g', linestyle = '--', label = "middle, right node") -plt.plot(disp3, times, color='b', label = "right node") +if has_matplotlib: + # plt.plot(disp0, times, color='k', label = "left node (fix)") + plt.plot(disp1, times, color='g', label="middle, left node") + plt.plot(disp2, times, color='g', linestyle='--', + label="middle, right node") + plt.plot(disp3, times, color='b', label="right node") -if(softDuration > 0): - plt.plot((maxMin[1], maxMin[0]), (switchT, switchT),) - plt.plot((maxMin[1], maxMin[0]), (switchEnd, switchEnd), ) + if(softDuration > 0): + plt.plot((maxMin[1], maxMin[0]), (switchT, switchT),) + plt.plot((maxMin[1], maxMin[0]), (switchEnd, switchEnd), ) -plt.title("Displacement in $x$ of the nodes") -plt.ylabel("Time [S]") -plt.xlabel("displacement [m]") -plt.xlim((maxMin[1] * 1.3, maxMin[0] * 1.1)) -plt.legend() -plt.show() + plt.title("Displacement in $x$ of the nodes") + plt.ylabel("Time [S]") + plt.xlabel("displacement [m]") + plt.xlim((maxMin[1] * 1.3, maxMin[0] * 1.1)) + plt.legend() + plt.show() # If the softening is disabled, then the displacement looks wierd. # Because the displacement first increases and then decreases. # In this case `softDuration > 0` holds. # -# However if the softening is enabled, it looks rather good. -# The left middle node will start to vibrate, because it is not pulled in the other direction. +# However if the softening is enabled, it looks rather good. The left middle +# node will start to vibrate, because it is not pulled in the other direction. diff --git a/examples/structural_mechanics/bernoulli_beam_2_example.cc b/examples/structural_mechanics/bernoulli_beam_2_example.cc index c2f250518..adc1d0be5 100644 --- a/examples/structural_mechanics/bernoulli_beam_2_example.cc +++ b/examples/structural_mechanics/bernoulli_beam_2_example.cc @@ -1,151 +1,153 @@ /** * @file bernoulli_beam_2_exemple.cc * * @author Fabian Barras * * @date creation: Mon Jan 18 2016 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #define TYPE _bernoulli_beam_2 using namespace akantu; // Linear load function static void lin_load(const Array & nodes, Array& load) { - for(auto &&data : zip(make_view(nodes, 2), make_view(load, 2))) { + for(auto &&data : zip(make_view(nodes, 2), make_view(load, 3))) { if (std::get<0>(data)[_y] <= 10) { std::get<1>(data)[_y] = -6000; } } } /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); // Defining the mesh Mesh beams(2); UInt nb_nodes = 3; UInt nb_nodes_1 = 1; UInt nb_nodes_2 = nb_nodes - nb_nodes_1 - 1; UInt nb_element = nb_nodes - 1; MeshAccessor mesh_accessor(beams); Array & nodes = mesh_accessor.getNodes(); nodes.resize(nb_nodes); beams.addConnectivityType(_bernoulli_beam_2); Array & connectivity = mesh_accessor.getConnectivity(_bernoulli_beam_2); connectivity.resize(nb_element); for (UInt i = 0; i < nb_nodes; ++i) { nodes(i, 1) = 0; } for (UInt i = 0; i < nb_nodes_1; ++i) { nodes(i, 0) = 10. * i / ((Real)nb_nodes_1); } nodes(nb_nodes_1, 0) = 10; for (UInt i = 0; i < nb_nodes_2; ++i) { nodes(nb_nodes_1 + i + 1, 0) = 10 + 8. * (i + 1) / ((Real)nb_nodes_2); } for (UInt i = 0; i < nb_element; ++i) { connectivity(i, 0) = i; connectivity(i, 1) = i + 1; } + mesh_accessor.makeReady(); + // Defining the materials StructuralMechanicsModel model(beams); StructuralMaterial mat1; mat1.E = 3e10; mat1.I = 0.0025; mat1.A = 0.01; model.addMaterial(mat1); StructuralMaterial mat2; mat2.E = 3e10; mat2.I = 0.00128; mat2.A = 0.01; model.addMaterial(mat2); // Defining the forces model.initFull(); const Real M = -3600; // Momentum at 3 Array & forces = model.getExternalForce(); Array & displacement = model.getDisplacement(); Array & boundary = model.getBlockedDOFs(); const Array & N_M = model.getStress(_bernoulli_beam_2); Array & element_material = model.getElementMaterial(_bernoulli_beam_2); forces.zero(); displacement.zero(); for (UInt i = 0; i < nb_nodes_2; ++i) { element_material(i + nb_nodes_1) = 1; } forces(nb_nodes - 1, 2) += M; - Array load(nodes.size(), 2); + Array load(nodes.size(), 3); lin_load(nodes, load); - model.computeForcesByGlobalTractionArray<_bernoulli_beam_2>(load); + model.computeForcesByGlobalTractionArray(load, _bernoulli_beam_2); // Defining the boundary conditions boundary(0, 0) = true; boundary(0, 1) = true; boundary(0, 2) = true; boundary(nb_nodes_1, 1) = true; boundary(nb_nodes - 1, 1) = true; model.addDumpFieldVector("displacement"); model.addDumpField("rotation"); model.addDumpFieldVector("force"); model.addDumpField("momentum"); model.solveStep(); // Post-Processing std::cout << " d1 = " << displacement(nb_nodes_1, 2) << std::endl; std::cout << " d2 = " << displacement(nb_nodes - 1, 2) << std::endl; std::cout << " M1 = " << N_M(0, 1) << std::endl; std::cout << " M2 = " << N_M(2 * (nb_nodes - 2), 1) << std::endl; model.dump(); finalize(); } diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc index a7eaafa2e..87262a9f1 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc @@ -1,603 +1,602 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * * Copyright (©) 2010-2011 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_RVE.hh" #include "element_group.hh" #include "material_damage_iterative.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "non_local_manager.hh" #include "parser.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, - const ID & id, - const MemoryID & memory_id) - : SolidMechanicsModel(mesh, dim, id, memory_id), volume(0.), + const ID & id) + : SolidMechanicsModel(mesh, dim, id), volume(0.), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// find the four corner nodes of the RVE findCornerNodes(); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); bottom_nodes.insert(bottom.begin(), bottom.end()); const auto & left = mesh.getElementGroup("left").getNodeGroup(); left_nodes.insert(left.begin(), left.end()); // /// enforce periodicity on the displacement fluctuations // auto surface_pair_1 = std::make_pair("top", "bottom"); // auto surface_pair_2 = std::make_pair("right", "left"); // SurfacePairList surface_pairs_list; // surface_pairs_list.push_back(surface_pair_1); // surface_pairs_list.push_back(surface_pair_2); // TODO: To Nicolas correct the PBCs // this->setPBC(surface_pairs_list); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); auto options_cp(options); options_cp.analysis_method = AnalysisMethod::_static; SolidMechanicsModel::initFullImpl(options_cp); // this->initMaterials(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); /// compute the volume of the RVE GhostType gt = _not_ghost; for (auto element_type : this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { Array Volume(this->mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume = fem.integrate(Volume, element_type); } std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; - base_name << this->id; // << this->memory_id - 1; + base_name << this->id; this->setBaseName(base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField("stress"); this->addDumpField("grad_u"); this->addDumpField("eigen_grad_u"); this->addDumpField("blocked_dofs"); this->addDumpField("material_index"); this->addDumpField("damage"); this->addDumpField("Sc"); this->addDumpField("external_force"); this->addDumpField("equivalent_stress"); this->addDumpField("internal_force"); this->addDumpField("delta_T"); this->dump(); this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement that will /// be applied Vector x(spatial_dimension); Vector appl_disp(spatial_dimension); /// fix top right node UInt node = this->corner_nodes(2); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); // (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.; // (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.; /// apply Hx at all the other corner nodes; H: displ. gradient node = this->corner_nodes(0); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(1); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(3); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyHomogeneousTemperature( const Real & temperature) { for (UInt m = 0; m < this->getNbMaterials(); ++m) { Material & mat = this->getMaterial(m); const ElementTypeMapArray & filter_map = mat.getElementFilter(); // Loop over all element types for (auto && type : filter_map.elementTypes(spatial_dimension)) { const Array & filter = filter_map(type); if (filter.size() == 0) continue; Array & delta_T = mat.getArray("delta_T", type); Array::scalar_iterator delta_T_it = delta_T.begin(); Array::scalar_iterator delta_T_end = delta_T.end(); for (; delta_T_it != delta_T_end; ++delta_T_it) { *delta_T_it = temperature; } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { AKANTU_DEBUG_IN(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, spatial_dimension))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(3) = node; } } for (UInt i = 0; i < corner_nodes.size(); ++i) { if (corner_nodes(i) == UInt(-1)) AKANTU_ERROR("The corner node " << i + 1 << " wasn't found"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } /// advance the damage MaterialDamageIterative<2> & mat_paste = dynamic_cast &>(*this->materials[1]); MaterialDamageIterative<2> & mat_aggregate = dynamic_cast &>(*this->materials[0]); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", SolveConvergenceCriteria::_solution); do { this->solveStep(); /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; if (max_eq_stress_aggregate > max_eq_stress_paste) nb_damaged_elements = mat_aggregate.updateDamage(); else if (max_eq_stress_aggregate < max_eq_stress_paste) nb_damaged_elements = mat_paste.updateDamage(); else nb_damaged_elements = (mat_paste.updateDamage() + mat_aggregate.updateDamage()); std::cout << "the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); if (this->nb_dumps % 10 == 0) { this->dump(); } this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); Real average = 0; GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & stress_vec = this->materials[m]->getStress(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_stress_vec(k, row_index * spatial_dimension + col_index); // 3 is the value // for the yy (in // 3D, the value is // 4) } } else if (field_type == "strain") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & gradu_vec = this->materials[m]->getGradU(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } } else { AKANTU_ERROR("Averaging not implemented for this field!!!"); } } return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { AKANTU_DEBUG_IN(); const UInt dim = 2; AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5) /// clear the eigenstrain Matrix zero_eigengradu(dim, dim, 0.); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigengradu; } /// storage for results of 3 different loading states UInt voigt_size = VoigtHelper::size; Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); /// save the damage state before filling up cracks // ElementTypeMapReal saved_damage("saved_damage"); // saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value = // 0); // this->fillCracks(saved_damage); /// virtual test 1: H(0, 0) = 0.01; this->performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.zero(); H(1, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.zero(); H(0, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 2); /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); C_macro.mul(stresses, eps_inverse); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); this->applyBoundaryConditions(H); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", SolveConvergenceCriteria::_solution); this->solveStep(); /// get average stress and strain eff_stresses(0, test_no) = this->averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = this->averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = this->averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = this->averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = this->averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * this->averageTensorField(1, 0, "strain"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU( Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = this->averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = this->averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = this->averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = this->averageTensorField(1, 0, "eigen_grad_u"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); if (use_RVE_mat_selector) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real box_size = std::abs(top - bottom); Real eps = box_size * 1e-6; auto tmp = std::make_shared(*this, box_size, "gel", this->nb_gel_pockets, eps); tmp->setFallback(material_selector); material_selector = tmp; } this->assignMaterialToElements(); // synchronize the element material arrays this->synchronize(SynchronizationTag::_material_id); for (auto & material : materials) { /// init internals properties const auto type = material->getID(); if (type.find("material_FE2") != std::string::npos) continue; material->initMaterial(); } this->synchronize(SynchronizationTag::_smm_init_mat); if (this->non_local_manager) { this->non_local_manager->initialize(); } // SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = this->getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; Real E = mat->get("E"); auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; sav_dam = dam; for (auto q : arange(dam.size())) { E_homogenized = (E_gel - E) * dam(q) + E; dam(q) = 1. - (E_homogenized / E); } } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::drainCracks( const ElementTypeMapReal & saved_damage) { for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; dam = sav_dam; } } } } } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh index 68c37ee23..771229c53 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh @@ -1,232 +1,231 @@ /** * @file solid_mechanics_model_RVE.hh * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 14:54:18 2016 * * @brief SMM for RVE computations in FE2 simulations * * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_ /* -------------------------------------------------------------------------- */ #include "aka_grid_dynamic.hh" #include "solid_mechanics_model.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelRVE : public SolidMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector = true, UInt nb_gel_pockets = 400, UInt spatial_dimension = _all_dimensions, - const ID & id = "solid_mechanics_model", - const MemoryID & memory_id = 0); + const ID & id = "solid_mechanics_model"); virtual ~SolidMechanicsModelRVE(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: void initFullImpl(const ModelOptions & option) override; /// initialize the materials void initMaterials() override; public: /// apply boundary contions based on macroscopic deformation gradient virtual void applyBoundaryConditions(const Matrix & displacement_gradient); /// apply homogeneous temperature field from the macroscale level to the RVEs virtual void applyHomogeneousTemperature(const Real & temperature); /// advance the reactions -> grow gel and apply homogenized properties void advanceASR(const Matrix & prestrain); /// compute average stress or strain in the model Real averageTensorField(UInt row_index, UInt col_index, const ID & field_type); /// compute effective stiffness of the RVE void homogenizeStiffness(Matrix & C_macro); /// compute average eigenstrain void homogenizeEigenGradU(Matrix & eigen_gradu_macro); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ inline void unpackData(CommunicationBuffer & buffer, const Array & index, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); AKANTU_GET_MACRO(Volume, volume, Real); private: /// find the corner nodes void findCornerNodes(); /// perform virtual testing void performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no); void fillCracks(ElementTypeMapReal & saved_damage); void drainCracks(const ElementTypeMapReal & saved_damage); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ /// volume of the RVE Real volume; /// corner nodes 1, 2, 3, 4 (see Leonardo's thesis, page 98) Array corner_nodes; /// bottom nodes std::unordered_set bottom_nodes; /// left nodes std::unordered_set left_nodes; /// standard mat selector or user one bool use_RVE_mat_selector; /// the number of gel pockets inside the RVE UInt nb_gel_pockets; /// dump counter UInt nb_dumps; }; inline void SolidMechanicsModelRVE::unpackData(CommunicationBuffer & buffer, const Array & index, const SynchronizationTag & tag) { SolidMechanicsModel::unpackData(buffer, index, tag); // if (tag == SynchronizationTag::_smm_uv) { // auto disp_it = displacement->begin(spatial_dimension); // // for (auto node : index) { // Vector current_disp(disp_it[node]); // // // if node is at the bottom, u_bottom = u_top +u_2 -u_3 // if (bottom_nodes.count(node)) { // current_disp += Vector(disp_it[corner_nodes(1)]); // current_disp -= Vector(disp_it[corner_nodes(2)]); // } // // if node is at the left, u_left = u_right +u_4 -u_3 // else if (left_nodes.count(node)) { // current_disp += Vector(disp_it[corner_nodes(3)]); // current_disp -= Vector(disp_it[corner_nodes(2)]); // } // } // } } /* -------------------------------------------------------------------------- */ /* ASR material selector */ /* -------------------------------------------------------------------------- */ class GelMaterialSelector : public MeshDataMaterialSelector { public: GelMaterialSelector(SolidMechanicsModel & model, const Real box_size, const std::string & gel_material, const UInt nb_gel_pockets, Real /*tolerance*/ = 0.) : MeshDataMaterialSelector("physical_names", model), model(model), gel_material(gel_material), nb_gel_pockets(nb_gel_pockets), nb_placed_gel_pockets(0), box_size(box_size) { Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = model.getSpatialDimension(); Element el{_triangle_3, 0, _not_ghost}; UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); Array barycenter(nb_element, spatial_dimension); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { el.element = std::get<0>(data); auto & bary = std::get<1>(data); mesh.getBarycenter(el, bary); } /// generate the gel pockets srand(0.); Vector center(spatial_dimension); UInt placed_gel_pockets = 0; std::set checked_baries; while (placed_gel_pockets != nb_gel_pockets) { /// get a random bary center UInt bary_id = rand() % nb_element; if (checked_baries.find(bary_id) != checked_baries.end()) continue; checked_baries.insert(bary_id); el.element = bary_id; if (MeshDataMaterialSelector::operator()(el) == 1) continue; /// element belongs to paste gel_pockets.push_back(el); placed_gel_pockets += 1; } } UInt operator()(const Element & elem) { UInt temp_index = MeshDataMaterialSelector::operator()(elem); if (temp_index == 1) return temp_index; std::vector::const_iterator iit = gel_pockets.begin(); std::vector::const_iterator eit = gel_pockets.end(); if (std::find(iit, eit, elem) != eit) { nb_placed_gel_pockets += 1; std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; return model.getMaterialIndex(gel_material); ; } return 0; } protected: SolidMechanicsModel & model; std::string gel_material; std::vector gel_pockets; UInt nb_gel_pockets; UInt nb_placed_gel_pockets; Real box_size; }; } // namespace akantu ///#include "material_selector_tmpl.hh" #endif /* AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_ */ diff --git a/extra_packages/igfem/src/integrator_gauss_igfem.hh b/extra_packages/igfem/src/integrator_gauss_igfem.hh index bee55b2f7..c32092989 100644 --- a/extra_packages/igfem/src/integrator_gauss_igfem.hh +++ b/extra_packages/igfem/src/integrator_gauss_igfem.hh @@ -1,124 +1,123 @@ /** * @file integrator_gauss_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Gauss integration facilities for IGFEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_INTEGRATOR_IGFEM_HH_ #define AKANTU_INTEGRATOR_IGFEM_HH_ /* -------------------------------------------------------------------------- */ #include "integrator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template class IntegratorGauss<_ek_igfem, IOF> : public Integrator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - IntegratorGauss(const Mesh & mesh, const ID & id = "integrator_gauss", - const MemoryID & memory_id = 0); + IntegratorGauss(const Mesh & mesh, const ID & id = "integrator_gauss"); virtual ~IntegratorGauss(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initIntegrator(const Array & nodes, ElementType type, GhostType ghost_type); /// precompute jacobians on elements of type "type" template void precomputeJacobiansOnQuadraturePoints(const Array & nodes, GhostType ghost_type); /// integrate f on the element "elem" of type "type" template inline void integrateOnElement(const Array & f, Real * intf, - UInt nb_degree_of_freedom, const UInt elem, + UInt nb_degree_of_freedom, UInt elem, GhostType ghost_type) const; /// integrate f for all elements of type "type" template void integrate(const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const; /// integrate one element scalar value on all elements of type "type" template Real integrate(const Vector & in_f, UInt index, GhostType ghost_type) const; /// integrate scalar field in_f template Real integrate(const Array & in_f, GhostType ghost_type, const Array & filter_elements) const; /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q * /// w_q @f$) template void integrateOnIntegrationPoints(const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const; /// return a vector with quadrature points natural coordinates template const Matrix & getIntegrationPoints(GhostType ghost_type) const; /// return the number of quadrature points for a given element type template inline UInt getNbIntegrationPoints(GhostType ghost_type = _not_ghost) const; /// compute the vector of quadrature points natural coordinates template void computeQuadraturePoints(GhostType ghost_type); /// check that the jacobians are not negative template void checkJacobians(GhostType ghost_type) const; public: /// compute the jacobians on quad points for a given element template void computeJacobianOnQuadPointsByElement(const Matrix & node_coords, Vector & jacobians); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: inline void integrate(Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const; private: ElementTypeMap> quadrature_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss_igfem_inline_impl.hh" } // namespace akantu #endif /*AKANTU_INTEGRATOR_IGFEM_HH_ */ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ diff --git a/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh b/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh index 9b92acc8e..376ccba79 100644 --- a/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh +++ b/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh @@ -1,451 +1,451 @@ /** * @file integrator_gauss_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Inline functions of gauss integrator for the case of IGFEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ } // namespace akantu #include "fe_engine.hh" #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif namespace akantu { /* -------------------------------------------------------------------------- */ #define INIT_INTEGRATOR(type) \ computeQuadraturePoints(ghost_type); \ precomputeJacobiansOnQuadraturePoints(nodes, ghost_type); \ checkJacobians(ghost_type); template inline void IntegratorGauss<_ek_igfem, IOF>::initIntegrator(const Array & nodes, ElementType type, GhostType ghost_type) { AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(INIT_INTEGRATOR); } #undef INIT_INTEGRATOR /* -------------------------------------------------------------------------- */ template template inline UInt IntegratorGauss<_ek_igfem, IOF>::getNbIntegrationPoints( GhostType) const { const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); return (nb_quad_points_sub_1 + nb_quad_points_sub_2); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrateOnElement( const Array & f, Real * intf, UInt nb_degree_of_freedom, const UInt elem, GhostType ghost_type) const { Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = getNbIntegrationPoints(); AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom, "The vector f do not have the good number of component."); Real * f_val = f.storage() + elem * f.getNbComponent(); Real * jac_val = jac_loc.storage() + elem * nb_quadrature_points; integrate(f_val, jac_val, intf, nb_degree_of_freedom, nb_quadrature_points); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss<_ek_igfem, IOF>::integrate( const Vector & in_f, UInt index, GhostType ghost_type) const { const Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = getNbIntegrationPoints(); AKANTU_DEBUG_ASSERT(in_f.size() == nb_quadrature_points, "The vector f do not have nb_quadrature_points entries."); Real * jac_val = jac_loc.storage() + index * nb_quadrature_points; Real intf; integrate(in_f.storage(), jac_val, &intf, 1, nb_quadrature_points); return intf; return 0.; } /* -------------------------------------------------------------------------- */ template inline void IntegratorGauss<_ek_igfem, IOF>::integrate(Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const { memset(inte, 0, nb_degree_of_freedom * sizeof(Real)); Real * cjac = jac; for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt dof = 0; dof < nb_degree_of_freedom; ++dof) { inte[dof] += *f * *cjac; ++f; } ++cjac; } } /* -------------------------------------------------------------------------- */ template template inline const Matrix & IntegratorGauss<_ek_igfem, IOF>::getIntegrationPoints( GhostType ghost_type) const { AKANTU_DEBUG_ASSERT( quadrature_points.exists(type, ghost_type), "Quadrature points for type " << quadrature_points.printType(type, ghost_type) << " have not been initialized." << " Did you use 'computeQuadraturePoints' function ?"); return quadrature_points(type, ghost_type); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::computeQuadraturePoints( GhostType ghost_type) { /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; /// store the quadrature points on the two subelements Matrix & quads_sub_1 = quadrature_points(sub_type_1, ghost_type); Matrix & quads_sub_2 = quadrature_points(sub_type_2, ghost_type); quads_sub_1 = GaussIntegrationElement::getQuadraturePoints(); quads_sub_2 = GaussIntegrationElement::getQuadraturePoints(); /// store all quad points for the current type UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt spatial_dimension = mesh.getSpatialDimension(); Matrix & quads = quadrature_points(type, ghost_type); quads = Matrix(spatial_dimension, nb_quad_points_sub_1 + nb_quad_points_sub_2); Matrix quads_1(quads.storage(), quads.rows(), nb_quad_points_sub_1); quads_1 = quads_sub_1; Matrix quads_2(quads.storage() + quads.rows() * nb_quad_points_sub_1, quads.rows(), nb_quad_points_sub_2); quads_2 = quads_sub_2; } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::computeJacobianOnQuadPointsByElement( const Matrix & node_coords, Vector & jacobians) { /// optimize: get the matrix from the ElementTypeMap Matrix quad = GaussIntegrationElement::getQuadraturePoints(); // jacobian ElementClass::computeJacobian(quad, node_coords, jacobians); } /* -------------------------------------------------------------------------- */ template inline IntegratorGauss<_ek_igfem, IOF>::IntegratorGauss( - const Mesh & mesh, const ID & id, const MemoryID & memory_id) - : Integrator(mesh, id, memory_id) { + const Mesh & mesh, const ID & id) + : Integrator(mesh, id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::checkJacobians( GhostType ghost_type) const { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points = nb_quad_points_sub_1 + nb_quad_points_sub_2; UInt nb_element; nb_element = mesh.getConnectivity(type, ghost_type).getSize(); Real * jacobians_val = jacobians(type, ghost_type).storage(); for (UInt i = 0; i < nb_element * nb_quadrature_points; ++i, ++jacobians_val) { if (*jacobians_val < 0) AKANTU_ERROR( "Negative jacobian computed," << " possible problem in the element node ordering (Quadrature Point " << i % nb_quadrature_points << ":" << i / nb_quadrature_points << ":" << type << ":" << ghost_type << ")"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::precomputeJacobiansOnQuadraturePoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// get the number of nodes for the subelements and the parent element UInt nb_nodes_sub_1 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_2 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_quadrature_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points = nb_quadrature_points_sub_1 + nb_quadrature_points_sub_2; UInt nb_element = mesh.getNbElement(type, ghost_type); Array * jacobians_tmp; if (!jacobians.exists(type, ghost_type)) jacobians_tmp = &jacobians.alloc(nb_element * nb_quadrature_points, 1, type, ghost_type); else { jacobians_tmp = &jacobians(type, ghost_type); jacobians_tmp->resize(nb_element * nb_quadrature_points); } Array::vector_iterator jacobians_it = jacobians_tmp->begin_reinterpret(nb_quadrature_points, nb_element); Vector weights_sub_1 = GaussIntegrationElement::getWeights(); Vector weights_sub_2 = GaussIntegrationElement::getWeights(); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Array::const_matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); // Matrix local_coord(spatial_dimension, nb_nodes_per_element); for (UInt elem = 0; elem < nb_element; ++elem, ++jacobians_it, ++x_it) { const Matrix & X = *x_it; Matrix sub_1_coords(spatial_dimension, nb_nodes_sub_1); Matrix sub_2_coords(spatial_dimension, nb_nodes_sub_2); ElementClass::getSubElementCoords(X, sub_1_coords, 0); ElementClass::getSubElementCoords(X, sub_2_coords, 1); Vector & J = *jacobians_it; /// initialize vectors to store the jacobians for each subelement Vector J_sub_1(nb_quadrature_points_sub_1); Vector J_sub_2(nb_quadrature_points_sub_2); computeJacobianOnQuadPointsByElement(sub_1_coords, J_sub_1); computeJacobianOnQuadPointsByElement(sub_2_coords, J_sub_2); J_sub_1 *= weights_sub_1; J_sub_2 *= weights_sub_2; /// copy results into the jacobian vector for this element for (UInt i = 0; i < nb_quadrature_points_sub_1; ++i) { J(i) = J_sub_1(i); } for (UInt i = 0; i < nb_quadrature_points_sub_2; ++i) { J(i + nb_quadrature_points_sub_1) = J_sub_2(i); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrate( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); const Matrix & quads = quadrature_points(type, ghost_type); UInt nb_points = quads.cols(); const Array & jac_loc = jacobians(type, ghost_type); Array::const_matrix_iterator J_it; Array::matrix_iterator inte_it; Array::const_matrix_iterator f_it; UInt nb_element; Array * filtered_J = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); const Array & cfiltered_J = *filtered_J; // \todo temporary patch J_it = cfiltered_J.begin_reinterpret(nb_points, 1, nb_element); } else { nb_element = mesh.getNbElement(type, ghost_type); J_it = jac_loc.begin_reinterpret(nb_points, 1, nb_element); } intf.resize(nb_element); f_it = in_f.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element); inte_it = intf.begin_reinterpret(nb_degree_of_freedom, 1, nb_element); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Matrix & f = *f_it; const Matrix & J = *J_it; Matrix & inte_f = *inte_it; inte_f.mul(f, J); } delete filtered_J; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss<_ek_igfem, IOF>::integrate( const Array & in_f, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); Array intfv(0, 1); integrate(in_f, intfv, 1, ghost_type, filter_elements); UInt nb_values = intfv.getSize(); if (nb_values == 0) return 0.; UInt nb_values_to_sum = nb_values >> 1; std::sort(intfv.begin(), intfv.end()); // as long as the half is not empty while (nb_values_to_sum) { UInt remaining = (nb_values - 2 * nb_values_to_sum); if (remaining) intfv(nb_values - 2) += intfv(nb_values - 1); // sum to consecutive values and store the sum in the first half for (UInt i = 0; i < nb_values_to_sum; ++i) { intfv(i) = intfv(2 * i) + intfv(2 * i + 1); } nb_values = nb_values_to_sum; nb_values_to_sum >>= 1; } AKANTU_DEBUG_OUT(); return intfv(0); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrateOnIntegrationPoints( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); UInt nb_element; const Matrix & quads = quadrature_points(type, ghost_type); UInt nb_points = quads.cols(); const Array & jac_loc = jacobians(type, ghost_type); Array::const_scalar_iterator J_it; Array::vector_iterator inte_it; Array::const_vector_iterator f_it; Array * filtered_J = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); J_it = filtered_J->begin(); } else { nb_element = mesh.getNbElement(type, ghost_type); J_it = jac_loc.begin(); } intf.resize(nb_element * nb_points); f_it = in_f.begin(nb_degree_of_freedom); inte_it = intf.begin(nb_degree_of_freedom); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Real & J = *J_it; const Vector & f = *f_it; Vector & inte_f = *inte_it; inte_f = f; inte_f *= J; } delete filtered_J; AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/igfem/src/non_local_manager_igfem.cc b/extra_packages/igfem/src/non_local_manager_igfem.cc index 5f1b2b851..c530c2e90 100644 --- a/extra_packages/igfem/src/non_local_manager_igfem.cc +++ b/extra_packages/igfem/src/non_local_manager_igfem.cc @@ -1,307 +1,306 @@ /** * @file non_local_manager_igfem.cc * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 15:32:10 2015 * * @brief Implementation of non-local manager igfem * * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_DAMAGE_NON_LOCAL #include "non_local_manager_igfem.hh" #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NonLocalManagerIGFEM::NonLocalManagerIGFEM(SolidMechanicsModelIGFEM & model, - const ID & id, - const MemoryID & memory_id) - : NonLocalManager(model, id, memory_id) { + const ID & id) + : NonLocalManager(model, id) { Mesh & mesh = this->model.getMesh(); /// initialize the element type map array /// it will be resized to nb_quad * nb_element during the computation of /// coords mesh.initElementTypeMapArray(quad_positions, spatial_dimension, spatial_dimension, false, _ek_igfem, true); } /* -------------------------------------------------------------------------- */ -NonLocalManagerIGFEM::~NonLocalManagerIGFEM() {} +NonLocalManagerIGFEM::~NonLocalManagerIGFEM() =default; /* -------------------------------------------------------------------------- */ void NonLocalManagerIGFEM::init() { /// store the number of current ghost elements for each type in the mesh ElementTypeMap nb_ghost_protected; Mesh & mesh = this->model.getMesh(); for (UInt k = _ek_regular; k <= _ek_igfem; ++k) { ElementKind el_kind = (ElementKind)k; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost, el_kind); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost, el_kind); for (; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); } /// exchange the missing ghosts for the non-local neighborhoods this->createNeighborhoodSynchronizers(); /// insert the ghost quadrature points of the non-local materials into the /// non-local neighborhoods for (UInt m = 0; m < this->non_local_materials.size(); ++m) { switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])) .insertQuadsInNeighborhoods(_ghost); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])) .insertQuadsInNeighborhoods(_ghost); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])) .insertQuadsInNeighborhoods(_ghost); break; } } FEEngine & fee_regular = this->model.getFEEngine(); FEEngine & fee_igfem = this->model.getFEEngine("IGFEMFEEngine"); this->updatePairLists(); /// cleanup the unneccessary ghost elements this->cleanupExtraGhostElements(nb_ghost_protected); this->initElementTypeMap(1, volumes, fee_regular, _ek_regular); this->initElementTypeMap(1, volumes, fee_igfem, _ek_igfem); this->setJacobians(fee_regular, _ek_regular); this->setJacobians(fee_igfem, _ek_igfem); this->initNonLocalVariables(); this->computeWeights(); } /* -------------------------------------------------------------------------- */ void NonLocalManagerIGFEM::computeAllNonLocalStresses() { /// update the flattened version of the internals std::map::iterator non_local_variable_it = non_local_variables.begin(); std::map::iterator non_local_variable_end = non_local_variables.end(); for (; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { non_local_variable_it->second->local.zero(); non_local_variable_it->second->non_local.zero(); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType)gt; this->flattenInternal(non_local_variable_it->second->local, ghost_type, _ek_regular); this->flattenInternal(non_local_variable_it->second->local, ghost_type, _ek_igfem); } } this->volumes.zero(); /// loop over all the neighborhoods and compute intiate the /// exchange of the non-local_variables std::set::const_iterator global_neighborhood_it = global_neighborhoods.begin(); NeighborhoodMap::iterator it; for (; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { it = neighborhoods.find(*global_neighborhood_it); if (it != neighborhoods.end()) it->second->getSynchronizerRegistry().asynchronousSynchronize( SynchronizationTag::_mnl_for_average); else dummy_synchronizers[*global_neighborhood_it]->asynchronousSynchronize( dummy_accessor, SynchronizationTag::_mnl_for_average); } this->averageInternals(_not_ghost); AKANTU_DEBUG_INFO("Wait distant non local stresses"); /// loop over all the neighborhoods and block until all non-local /// variables have been exchanged global_neighborhood_it = global_neighborhoods.begin(); it = neighborhoods.begin(); for (; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { it = neighborhoods.find(*global_neighborhood_it); if (it != neighborhoods.end()) it->second->getSynchronizerRegistry().waitEndSynchronize( SynchronizationTag::_mnl_for_average); else dummy_synchronizers[*global_neighborhood_it]->waitEndSynchronize( dummy_accessor, SynchronizationTag::_mnl_for_average); } this->averageInternals(_ghost); /// copy the results in the materials this->distributeInternals(_ek_regular); /// loop over all the materials and update the weights for (UInt m = 0; m < this->non_local_materials.size(); ++m) { switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])) .computeNonLocalStresses(_not_ghost); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])) .computeNonLocalStresses(_not_ghost); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])) .computeNonLocalStresses(_not_ghost); break; } } ++this->compute_stress_calls; } /* -------------------------------------------------------------------------- */ void NonLocalManagerIGFEM::cleanupExtraGhostElements( ElementTypeMap & nb_ghost_protected) { typedef std::set ElementSet; ElementSet relevant_ghost_elements; ElementSet to_keep_per_neighborhood; /// loop over all the neighborhoods and get their protected ghosts NeighborhoodMap::iterator neighborhood_it = neighborhoods.begin(); NeighborhoodMap::iterator neighborhood_end = neighborhoods.end(); for (; neighborhood_it != neighborhood_end; ++neighborhood_it) { neighborhood_it->second->cleanupExtraGhostElements( to_keep_per_neighborhood); ElementSet::const_iterator it = to_keep_per_neighborhood.begin(); for (; it != to_keep_per_neighborhood.end(); ++it) relevant_ghost_elements.insert(*it); to_keep_per_neighborhood.zero(); } /// remove all unneccessary ghosts from the mesh /// Create list of element to remove and new numbering for element to keep Mesh & mesh = this->model.getMesh(); ElementSet ghost_to_erase; RemovedElementsEvent remove_elem(mesh); Element element; for (UInt k = _ek_regular; k < _ek_igfem; ++k) { ElementKind el_kind = (ElementKind)k; element.kind = _ek_igfem; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost, el_kind); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost, el_kind); element.ghost_type = _ghost; for (; it != last_type; ++it) { element.type = *it; UInt nb_ghost_elem = mesh.getNbElement(*it, _ghost); UInt nb_ghost_elem_protected = 0; try { nb_ghost_elem_protected = nb_ghost_protected(*it, _ghost); } catch (...) { } if (!remove_elem.getNewNumbering().exists(*it, _ghost)) remove_elem.getNewNumbering().alloc(nb_ghost_elem, 1, *it, _ghost); else remove_elem.getNewNumbering(*it, _ghost).resize(nb_ghost_elem); Array & new_numbering = remove_elem.getNewNumbering(*it, _ghost); for (UInt g = 0; g < nb_ghost_elem; ++g) { element.element = g; if (element.element >= nb_ghost_elem_protected && relevant_ghost_elements.find(element) == relevant_ghost_elements.end()) { remove_elem.getList().push_back(element); new_numbering(element.element) = UInt(-1); } } /// renumber remaining ghosts UInt ng = 0; for (UInt g = 0; g < nb_ghost_elem; ++g) { if (new_numbering(g) != UInt(-1)) { new_numbering(g) = ng; ++ng; } } } } for (UInt k = _ek_regular; k < _ek_igfem; ++k) { ElementKind el_kind = (ElementKind)k; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost, el_kind); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost, el_kind); for (; it != last_type; ++it) { UInt nb_elem = mesh.getNbElement(*it, _not_ghost); if (!remove_elem.getNewNumbering().exists(*it, _not_ghost)) remove_elem.getNewNumbering().alloc(nb_elem, 1, *it, _not_ghost); Array & new_numbering = remove_elem.getNewNumbering(*it, _not_ghost); for (UInt e = 0; e < nb_elem; ++e) { new_numbering(e) = e; } } } mesh.sendEvent(remove_elem); } /* -------------------------------------------------------------------------- */ void NonLocalManagerIGFEM::onElementsAdded(__attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { FEEngine & fee = this->model.getFEEngine("IGFEMFEEngine"); this->resizeElementTypeMap(1, volumes, fee, _ek_igfem); this->resizeElementTypeMap(spatial_dimension, quad_positions, fee, _ek_igfem); NonLocalManager::onElementsAdded(element_list, event); } /* -------------------------------------------------------------------------- */ void NonLocalManagerIGFEM::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { FEEngine & fee = this->model.getFEEngine("IGFEMFEEngine"); this->removeIntegrationPointsFromMap(event.getNewNumbering(), spatial_dimension, quad_positions, fee, _ek_igfem); this->removeIntegrationPointsFromMap(event.getNewNumbering(), 1, volumes, fee, _ek_igfem); NonLocalManager::onElementsRemoved(element_list, new_numbering, event); } } // namespace akantu #endif /* AKANTU_DAMAGE_NON_LOCAL */ diff --git a/extra_packages/igfem/src/non_local_manager_igfem.hh b/extra_packages/igfem/src/non_local_manager_igfem.hh index 1937a2c86..ae2838523 100644 --- a/extra_packages/igfem/src/non_local_manager_igfem.hh +++ b/extra_packages/igfem/src/non_local_manager_igfem.hh @@ -1,97 +1,96 @@ /** * @file non_local_manager_igfem.hh * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 14:21:33 2015 * * @brief Class that manages all the non-local neighborhoods for IGFEM * simulations * * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_DAMAGE_NON_LOCAL #ifndef AKANTU_NON_LOCAL_MANAGER_IGFEM_HH_ #define AKANTU_NON_LOCAL_MANAGER_IGFEM_HH_ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_igfem.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class NonLocalManagerIGFEM : public NonLocalManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLocalManagerIGFEM(SolidMechanicsModelIGFEM & model, - const ID & id = "non_local_manager_igfem", - const MemoryID & memory_id = 0); + const ID & id = "non_local_manager_igfem"); virtual ~NonLocalManagerIGFEM(); /* -------------------------------------------------------------------------- */ /* Methods */ /* -------------------------------------------------------------------------- */ public: /// initialize the non-local manager: compute pair lists and weights for all /// neighborhoods virtual void init(); /// average the internals and compute the non-local stresses virtual void computeAllNonLocalStresses(); /* -------------------------------------------------------------------------- */ /* MeshEventHandler inherited members */ /* -------------------------------------------------------------------------- */ virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); virtual void onElementsAdded(const Array & element_list, const NewElementsEvent & event); private: /// cleanup unneccessary ghosts virtual void cleanupExtraGhostElements(ElementTypeMap & nb_ghost_protected); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #endif /* AKANTU_NON_LOCAL_MANAGER_IGFEM_HH_ */ #endif /* AKANTU_DAMAGE_NON_LOCAL */ diff --git a/extra_packages/igfem/src/shape_igfem.cc b/extra_packages/igfem/src/shape_igfem.cc index d698a9cdb..f75a5407e 100644 --- a/extra_packages/igfem/src/shape_igfem.cc +++ b/extra_packages/igfem/src/shape_igfem.cc @@ -1,96 +1,94 @@ /** * @file shape_igfem_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief ShapeIGFEM inline implementation * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ -#include "aka_memory.hh" //#include "mesh.hh" #include "shape_igfem.hh" /* -------------------------------------------------------------------------- */ #if defined(AKANTU_IGFEM) namespace akantu { /* -------------------------------------------------------------------------- */ -ShapeLagrange<_ek_igfem>::ShapeLagrange(const Mesh & mesh, const ID & id, - const MemoryID & memory_id) - : ShapeFunctions(mesh, id, memory_id), - shapes("shapes_generic", id, memory_id), - shapes_derivatives("shapes_derivatives_generic", id, memory_id), - igfem_integration_points("igfem_integration_points", id, memory_id), - shapes_at_enrichments("shapes_at_enrichments", id, memory_id) { +ShapeLagrange<_ek_igfem>::ShapeLagrange(const Mesh & mesh, const ID & id) + : ShapeFunctions(mesh, id), + shapes("shapes_generic", id), + shapes_derivatives("shapes_derivatives_generic", id), + igfem_integration_points("igfem_integration_points", id), + shapes_at_enrichments("shapes_at_enrichments", id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /*-------------------------------------------------------------------------- */ void ShapeLagrange<_ek_igfem>::extractValuesAtStandardNodes( const Array & nodal_values, Array & extracted_values, GhostType ghost_type) const { AKANTU_DEBUG_ASSERT(nodal_values.getNbComponent() == extracted_values.getNbComponent(), "The arrays are not of the same size!!!!!"); extracted_values.zero(); UInt spatial_dimension = mesh.getSpatialDimension(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_igfem); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type, _ek_igfem); for (; it != end; ++it) { ElementType type = *it; UInt nb_elements = mesh.getNbElement(type, ghost_type); UInt nb_parent_nodes = 0; UInt nb_nodes_per_element = 0; #define GET_NODES_INFO(type) \ const ElementType parent_type = \ ElementClassProperty::parent_element_type; \ nb_parent_nodes = \ ElementClass::getNbNodesPerInterpolationElement(); \ nb_nodes_per_element = \ ElementClass::getNbNodesPerInterpolationElement(); AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_NODES_INFO); #undef GET_NODES_INFO UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); for (UInt e = 0; e < nb_elements; ++e) { /// copy the value at standard nodes UInt offset = e * nb_nodes_per_element; for (UInt n = 0; n < nb_parent_nodes; ++n) { UInt node = conn_val[offset + n]; for (UInt i = 0; i < nodal_values.getNbComponent(); ++i) extracted_values(node, i) = nodal_values(node, i); } } } } /* -------------------------------------------------------------------------- */ void ShapeLagrange<_ek_igfem>::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Shapes Lagrange [" << std::endl; ShapeFunctions::printself(stream, indent + 1); shapes.printself(stream, indent + 1); shapes_derivatives.printself(stream, indent + 1); stream << space << "]" << std::endl; } } // namespace akantu #endif diff --git a/extra_packages/igfem/src/shape_igfem.hh b/extra_packages/igfem/src/shape_igfem.hh index a3db5d3dd..537f6b401 100644 --- a/extra_packages/igfem/src/shape_igfem.hh +++ b/extra_packages/igfem/src/shape_igfem.hh @@ -1,200 +1,199 @@ /** * @file shape_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief shape functions for interface-enriched generalized FEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "shape_functions.hh" #ifndef AKANTU_SHAPE_IGFEM_HH_ #define AKANTU_SHAPE_IGFEM_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template <> class ShapeLagrange<_ek_igfem> : public ShapeFunctions { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - ShapeLagrange(const Mesh & mesh, const ID & id = "shape_igfem", - const MemoryID & memory_id = 0); + ShapeLagrange(const Mesh & mesh, const ID & id = "shape_igfem"); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initShapeFunctions(const Array & nodes, const Matrix & integration_points, const Matrix & integration_points_1, const Matrix & integration_points_2, ElementType type, GhostType ghost_type); inline void interpolateEnrichmentsAllTypes(const Array & src, Array & dst, ElementType type, GhostType ghost_type) const; template inline void precomputeShapesOnEnrichedNodes(const Array & nodes, GhostType ghost_type); template void interpolateAtEnrichedNodes(const Array & src, Array & dst, GhostType ghost_type) const; /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// pre compute all shape derivatives on the element integration points from /// natural coordinates template void precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// interpolate nodal values on the integration points template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// interpolate on physical point template void interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, Vector & interpolated, GhostType ghost_type) const; /// compute the gradient of u on the integration points template void gradientOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// multiply a field by shape functions @f$ fts_{ij} = f_i * \varphi_j @f$ template void fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const; /// find natural coords in parent element from real coords provided an element template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, GhostType ghost_type = _not_ghost) const; /// find natural coords in sub-element from real coords provided an element template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, UInt sub_element, GhostType ghost_type = _not_ghost) const; /// return true if the coordinates provided are inside the element, false /// otherwise template bool contains(const Vector & real_coords, UInt elem, GhostType ghost_type) const; /// compute the shape on a provided point template void computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, GhostType ghost_type) const; /// compute the shape derivatives on a provided point template void computeShapeDerivatives(const Matrix & real_coords, UInt elem, Tensor3 & shapes, GhostType ghost_type) const; /// interpolate a field on a given physical point template void interpolateOnPhysicalPoint(const Vector & real_coords, UInt elem, const Array & field, Vector & interpolated, GhostType ghost_type) const; /// function to extract values at standard nodes and zero-out enriched values /// of a nodal field void extractValuesAtStandardNodes(const Array & nodal_values, Array & extracted_values, GhostType ghost_type) const; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute the shape derivatives on integration points for a given element template inline void computeShapeDerivativesOnCPointsByElement(const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get a the shapes vector inline const Array & getShapes(ElementType el_type, GhostType ghost_type = _not_ghost) const; /// get a the shapes derivatives vector inline const Array & getShapesDerivatives(ElementType el_type, GhostType ghost_type = _not_ghost) const; /// get a the shapes vector inline const Array & getShapesAtEnrichedNodes(ElementType el_type, GhostType ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// shape functions for all elements ElementTypeMapArray shapes; /// shape functions derivatives for all elements ElementTypeMapArray shapes_derivatives; /// additional integration points for the IGFEM formulation ElementTypeMapArray igfem_integration_points; /// values of shape functions for all elements on the enriched nodes ElementTypeMapArray shapes_at_enrichments; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "shape_igfem_inline_impl.hh" /// standard output stream operator // template // inline std::ostream & operator <<(std::ostream & stream, const // ShapeIGFEM & _this) // { // _this.printself(stream); // return stream; // } #endif /* AKANTU_SHAPE_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/solid_mechanics_model_igfem.cc b/extra_packages/igfem/src/solid_mechanics_model_igfem.cc index f7dfab8e3..8422cb2d2 100644 --- a/extra_packages/igfem/src/solid_mechanics_model_igfem.cc +++ b/extra_packages/igfem/src/solid_mechanics_model_igfem.cc @@ -1,586 +1,585 @@ /** * @file solid_mechanics_model_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief solid mechanics model for IGFEM analysis * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_igfem.hh" #include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.hh" #include "igfem_helper.hh" #include "material_igfem.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_igfem_element_partition.hh" #include "dumper_igfem_elemental_field.hh" #include "dumper_igfem_material_internal_field.hh" #include "dumper_material_padders.hh" #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { const SolidMechanicsModelIGFEMOptions default_solid_mechanics_model_igfem_options(_static, false); SolidMechanicsModelIGFEM::SolidMechanicsModelIGFEM(Mesh & mesh, UInt dim, - const ID & id, - const MemoryID & memory_id) - : SolidMechanicsModel(mesh, dim, id, memory_id), IGFEMEnrichment(mesh), + const ID & id) + : SolidMechanicsModel(mesh, dim, id), IGFEMEnrichment(mesh), global_ids_updater(NULL) { AKANTU_DEBUG_IN(); delete material_selector; material_selector = new DefaultMaterialIGFEMSelector(*this); this->registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("igfem elements", id); this->mesh.addDumpMeshToDumper("igfem elements", mesh, spatial_dimension, _not_ghost, _ek_igfem); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelIGFEM::~SolidMechanicsModelIGFEM() { AKANTU_DEBUG_IN(); if (global_ids_updater) delete global_ids_updater; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::initFull(const ModelOptions & options) { AKANTU_DEBUG_IN(); /// intialize the IGFEM enrichment this->initialize(); SolidMechanicsModel::initFull(options); // set the initial condition to 0 real_force->clear(); real_displacement->clear(); real_residual->clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModelIGFEM::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_rdisp; sstr_rdisp << id << ":real_displacement"; std::stringstream sstr_rforc; sstr_rforc << id << ":real_force"; std::stringstream sstr_rresi; sstr_rresi << id << ":real_residual"; real_displacement = &(alloc(sstr_rdisp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); real_force = &(alloc(sstr_rforc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); real_residual = &(alloc(sstr_rresi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); SolidMechanicsModel::initArrays(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { SolidMechanicsModel::initParallel(partition, data_accessor); this->intersector_sphere.setDistributedSynchronizer(synch_parallel); if (mesh.isDistributed()) global_ids_updater = new GlobalIdsUpdater(mesh, synch_parallel); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first igfem material UInt igfem_index = 0; while ((dynamic_cast(materials[igfem_index]) == NULL) && igfem_index <= materials.size()) ++igfem_index; AKANTU_DEBUG_ASSERT(igfem_index != materials.size(), "No igfem materials in the material input file"); DefaultMaterialIGFEMSelector * igfem_mat_selector = dynamic_cast(material_selector); if (igfem_mat_selector != NULL) igfem_mat_selector->setIGFEMFallback(igfem_index); SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model, basically pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModelIGFEM::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("IGFEMFEEngine", mesh, spatial_dimension); /// insert the two feengines associated with the model in the map this->fe_engines_per_kind[_ek_regular] = &(this->getFEEngine()); this->fe_engines_per_kind[_ek_igfem] = &(this->getFEEngine("IGFEMFEEngine")); /// add the igfem type connectivities for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType type_ghost = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, type_ghost); Mesh::type_iterator last = mesh.lastType(spatial_dimension, type_ghost); for (; it != last; ++it) { const Array & connectivity = mesh.getConnectivity(*it, type_ghost); if (connectivity.getSize() != 0) { ElementType type = *it; Vector types_igfem = FEEngine::getIGFEMElementTypes(type); for (UInt i = 0; i < types_igfem.size(); ++i) mesh.addConnectivityType(types_igfem(i), type_ghost); } } } getFEEngine("IGFEMFEEngine").initShapeFunctions(_not_ghost); getFEEngine("IGFEMFEEngine").initShapeFunctions(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::onElementsAdded(const Array & elements, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); const NewIGFEMElementsEvent * igfem_event = dynamic_cast(&event); /// insert the new and old elements in the map if (igfem_event != NULL) { this->element_map.zero(); const Array & old_elements = igfem_event->getOldElementsList(); for (UInt e = 0; e < elements.getSize(); ++e) { this->element_map[elements(e)] = old_elements(e); } } /// update shape functions getFEEngine("IGFEMFEEngine").initShapeFunctions(_not_ghost); getFEEngine("IGFEMFEEngine").initShapeFunctions(_ghost); SolidMechanicsModel::onElementsAdded(elements, event); this->reassignMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine("IGFEMFEEngine").initShapeFunctions(_not_ghost); this->getFEEngine("IGFEMFEEngine").initShapeFunctions(_ghost); SolidMechanicsModel::onElementsRemoved(element_list, new_numbering, event); if (synch_parallel) synch_parallel->computeAllBufferSizes(*this); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); const NewIGFEMNodesEvent * igfem_event = dynamic_cast(&event); // update the node type if (igfem_event != NULL) { intersector_sphere.updateNodeType( nodes_list, igfem_event->getNewNodePerElem(), igfem_event->getElementType(), igfem_event->getGhostType()); } UInt nb_nodes = mesh.getNbNodes(); if (real_displacement) real_displacement->resize(nb_nodes); if (real_force) real_force->resize(nb_nodes); if (real_residual) real_residual->resize(nb_nodes); if (mesh.isDistributed()) mesh.getGlobalNodesIds().resize(mesh.getNbNodes()); if (displacement) displacement->resize(nb_nodes); if (mass) mass->resize(nb_nodes); if (velocity) velocity->resize(nb_nodes); if (acceleration) acceleration->resize(nb_nodes); if (force) force->resize(nb_nodes); if (residual) residual->resize(nb_nodes); if (blocked_dofs) blocked_dofs->resize(nb_nodes); if (previous_displacement) previous_displacement->resize(nb_nodes); if (increment_acceleration) increment_acceleration->resize(nb_nodes); if (increment) increment->resize(nb_nodes); if (current_position) current_position->resize(nb_nodes); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) { if (real_displacement) mesh.removeNodesFromArray(*real_displacement, new_numbering); if (real_force) mesh.removeNodesFromArray(*real_force, new_numbering); if (real_residual) mesh.removeNodesFromArray(*real_residual, new_numbering); // communicate global connectivity for slave nodes if (global_ids_updater) global_ids_updater->updateGlobalIDs( mesh.getNbNodes() - intersector_sphere.getNbStandardNodes()); SolidMechanicsModel::onNodesRemoved(nodes_list, new_numbering, event); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); ElementKind _element_kind = element_kind; if (dumper_name == "igfem elements") { _element_kind = _ek_igfem; } SolidMechanicsModel::addDumpGroupFieldToDumper( dumper_name, field_id, group_name, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::onDump() { this->computeValuesOnEnrichedNodes(); this->flattenAllRegisteredInternals(_ek_igfem); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumpers::Field * SolidMechanicsModelIGFEM::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, ElementKind kind) { dumpers::Field * field = NULL; if (kind != _ek_igfem) field = SolidMechanicsModel::createElementalField( field_name, group_name, padding_flag, spatial_dimension, kind); else { if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, spatial_dimension, kind); else if (field_name == "material_index") field = mesh.createElementalField( material_index, group_name, spatial_dimension, kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); if (field_name == "strain" || field_name == "Green strain" || field_name == "principal strain" || field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; bool is_internal = this->isInternal(field_name_copy, kind); if (is_internal) { ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy, kind); field = mesh.createElementalField( internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem); if (field_name == "strain") { dumpers::ComputeStrain * foo = new dumpers::ComputeStrain(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Von Mises stress") { dumpers::ComputeVonMisesStress * foo = new dumpers::ComputeVonMisesStress(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Green strain") { dumpers::ComputeStrain * foo = new dumpers::ComputeStrain(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal strain") { dumpers::ComputePrincipalStrain * foo = new dumpers::ComputePrincipalStrain(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal Green strain") { dumpers::ComputePrincipalStrain * foo = new dumpers::ComputePrincipalStrain(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } /// treat the paddings if (padding_flag) { if (field_name == "stress") { if (spatial_dimension == 2) { dumpers::StressPadder<2> * foo = new dumpers::StressPadder<2>(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } } else if (field_name == "strain" || field_name == "Green strain") { if (spatial_dimension == 2) { dumpers::StrainPadder<2> * foo = new dumpers::StrainPadder<2>(*this); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } } } // homogenize the field dumpers::ComputeFunctorInterface * foo = dumpers::HomogenizerProxy::createHomogenizer(*field); field = dumpers::FieldComputeProxy::createFieldCompute(field, *foo); } } } // } return field; } /* -------------------------------------------------------------------------- */ dumpers::Field * SolidMechanicsModelIGFEM::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["real_displacement"] = real_displacement; dumpers::Field * field = NULL; if (padding_flag) field = mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = mesh.createNodalField(real_nodal_fields[field_name], group_name); if (field == NULL) return SolidMechanicsModel::createNodalFieldReal(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ dumpers::Field * SolidMechanicsModelIGFEM::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, ElementKind kind) { return NULL; } /* -------------------------------------------------------------------------- */ dumpers::Field * SolidMechanicsModelIGFEM::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::computeValuesOnEnrichedNodes() { for (UInt n = 0; n < mesh.getNbNodes(); ++n) { for (UInt s = 0; s < spatial_dimension; ++s) (*real_displacement)(n, s) = (*displacement)(n, s); } Element element; Vector real_coords(spatial_dimension); Vector interpolated(spatial_dimension); Array::const_vector_iterator r_displ_it = this->real_displacement->begin(spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { element.ghost_type = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, *gt, _ek_igfem); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_igfem); for (; it != last; ++it) { element.type = *it; UInt nb_element = mesh.getNbElement(*it, *gt); if (!nb_element) continue; UInt * elem_val = mesh.getConnectivity(*it, *gt).storage(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); Matrix displ_val(spatial_dimension, nb_nodes_per_element); UInt nb_enriched_nodes = IGFEMHelper::getNbEnrichedNodes(*it); UInt nb_parent_nodes = IGFEMHelper::getNbParentNodes(*it); for (UInt el = 0; el < nb_element; ++el) { element.element = el; /// get the node coordinates of the element mesh.extractNodalValuesFromElement( mesh.getNodes(), nodes_coord.storage(), elem_val + el * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); /// get the displacement values at the nodes of the element mesh.extractNodalValuesFromElement( *(this->displacement), displ_val.storage(), elem_val + el * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); for (UInt i = 0; i < nb_enriched_nodes; ++i) { /// coordinates of enriched node real_coords = nodes_coord(nb_parent_nodes + i); /// global index of the enriched node UInt idx = elem_val[el * nb_nodes_per_element + nb_parent_nodes + i]; /// compute the real displacement value this->getFEEngine("IGFEMFEEngine") .interpolate(real_coords, displ_val, interpolated, element); r_displ_it[idx] = interpolated; } } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::transferInternalValues( const ID & internal, std::vector & new_elements, Array & added_quads, Array & internal_values) { /// @todo sort the new elements by their corresponding old element type and /// old material!!! /// get the number of elements for which iternals need to be transfered UInt nb_new_elements = new_elements.size(); UInt nb_new_quads = added_quads.getSize() / nb_new_elements; Array::const_matrix_iterator quad_coords = added_quads.begin_reinterpret(this->spatial_dimension, nb_new_quads, nb_new_elements); UInt nb_internal_component = internal_values.getNbComponent(); Array::matrix_iterator internal_val = internal_values.begin_reinterpret( nb_internal_component, nb_new_quads, nb_new_elements); Vector default_values(nb_internal_component, 0.); for (UInt e = 0; e < nb_new_elements; ++e, ++quad_coords, ++internal_val) { Element new_element = new_elements[e]; Element old_element = this->element_map[new_element]; UInt mat_idx = (this->material_index( old_element.type, old_element.ghost_type))(old_element.element); Material & old_material = *(this->materials[mat_idx]); old_material.extrapolateInternal(internal, old_element, *quad_coords, *internal_val); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelIGFEM::applyEigenGradU( const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type) { AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() == spatial_dimension * spatial_dimension, "The prescribed grad_u is not of the good size"); std::vector::iterator mat_it; for (mat_it = this->materials.begin(); mat_it != this->materials.end(); ++mat_it) { MaterialIGFEM * mat_igfem = dynamic_cast(*mat_it); if (mat_igfem != NULL) mat_igfem->applyEigenGradU(prescribed_eigen_grad_u, material_name, ghost_type); else if ((*mat_it)->getName() == material_name) (*mat_it)->applyEigenGradU(prescribed_eigen_grad_u, ghost_type); } } } // namespace akantu diff --git a/extra_packages/igfem/src/solid_mechanics_model_igfem.hh b/extra_packages/igfem/src/solid_mechanics_model_igfem.hh index 101e40b67..4933932f8 100644 --- a/extra_packages/igfem/src/solid_mechanics_model_igfem.hh +++ b/extra_packages/igfem/src/solid_mechanics_model_igfem.hh @@ -1,197 +1,196 @@ /** * @file solid_mechanics_model_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief solid mechanics model for IGFEM analysis * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_IGFEM_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_IGFEM_HH_ #include "global_ids_updater.hh" #include "igfem_enrichment.hh" #include "solid_mechanics_model.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelIGFEMOptions : public SolidMechanicsModelOptions { SolidMechanicsModelIGFEMOptions(AnalysisMethod analysis_method = _static, bool no_init_materials = false) : SolidMechanicsModelOptions(analysis_method, no_init_materials) {} }; extern const SolidMechanicsModelIGFEMOptions default_solid_mechanics_model_igfem_options; /* -------------------------------------------------------------------------- */ /* Solid Mechanics Model for IGFEM analysis */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModelIGFEM : public SolidMechanicsModel, public SolidMechanicsModelEventHandler, public IGFEMEnrichment { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate MyFEEngineIGFEMType; typedef std::map ElementMap; typedef std::map FEEnginesPerKindMap; SolidMechanicsModelIGFEM(Mesh & mesh, UInt spatial_dimension = _all_dimensions, - const ID & id = "solid_mechanics_model_igfem", - const MemoryID & memory_id = 0); + const ID & id = "solid_mechanics_model_igfem"); virtual ~SolidMechanicsModelIGFEM(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the cohesive model virtual void initFull(const ModelOptions & options = default_solid_mechanics_model_igfem_options); /// initialize the model virtual void initModel(); /// initialize igfem material virtual void initMaterials(); /// register the tags associated with the parallel synchronizer virtual void initParallel(MeshPartition * partition, DataAccessor * data_accessor = NULL); /// allocate all vectors virtual void initArrays(); /// transfer internals from old to new elements void transferInternalValues(const ID & internal, std::vector & new_elements, Array & added_quads, Array & internal_values); /// compute the barycenter for a sub-element inline void getSubElementBarycenter(UInt element, UInt sub_element, ElementType type, Vector & barycenter, GhostType ghost_type) const; /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); private: /// compute the real values of displacement, force, etc. on the enriched nodes void computeValuesOnEnrichedNodes(); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event); virtual void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event); virtual void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag); virtual dumpers::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, ElementKind kind); virtual dumpers::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag); /* -------------------------------------------------------------------------- */ /* Accessors */ /* -------------------------------------------------------------------------- */ public: /// get the fe-engines per kind AKANTU_GET_MACRO(FEEnginesPerKind, fe_engines_per_kind, const FEEnginesPerKindMap &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// real displacements array Array * real_displacement; /// real forces array Array * real_force; /// real residuals array Array * real_residual; /// map between and new elements (needed when the interface is moving) ElementMap element_map; /// global connectivity ids updater GlobalIdsUpdater * global_ids_updater; /// map between element kind and corresponding FEEngine object FEEnginesPerKindMap fe_engines_per_kind; }; /* -------------------------------------------------------------------------- */ /* IGFEMMaterialSelector */ /* -------------------------------------------------------------------------- */ class DefaultMaterialIGFEMSelector : public DefaultMaterialSelector { public: DefaultMaterialIGFEMSelector(const SolidMechanicsModelIGFEM & model) : DefaultMaterialSelector(model.getMaterialByElement()), fallback_value_igfem(0) {} virtual UInt operator()(const Element & element) { if (Mesh::getKind(element.type) == _ek_igfem) return fallback_value_igfem; else return DefaultMaterialSelector::operator()(element); } void setIGFEMFallback(UInt f) { this->fallback_value_igfem = f; } protected: UInt fallback_value_igfem; }; } // namespace akantu #if defined(AKANTU_INCLUDE_INLINE_IMPL) #include "solid_mechanics_model_igfem_inline_impl.hh" #endif #endif /* AKANTU_SOLID_MECHANICS_MODEL_IGFEM_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh index 29ac24991..cdb31ae9b 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh @@ -1,107 +1,106 @@ /** * @file ntn_friclaw_coulomb.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief coulomb friction with \mu_s = \mu_k (constant) * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICLAW_COULOMB_HH_ #define AST_NTN_FRICLAW_COULOMB_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawCoulomb : public Regularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNFricLawCoulomb(NTNBaseContact & contact, const ID & id = "coulomb", - const MemoryID & memory_id = 0); + NTNFricLawCoulomb(NTNBaseContact & contact, const ID & id = "coulomb"); virtual ~NTNFricLawCoulomb(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // friction coefficient SynchronizedArray mu; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const NTNFricLawCoulomb & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_coulomb_tmpl.hh" #endif /* AST_NTN_FRICLAW_COULOMB_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh index 733a6c548..a1b626b67 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh @@ -1,167 +1,166 @@ /** * @file ntn_friclaw_coulomb_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of coulomb friction * * * Copyright (©) 2015-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 "dumper_nodal_field.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawCoulomb::NTNFricLawCoulomb(NTNBaseContact & contact, - const ID & id, - const MemoryID & memory_id) - : Regularisation(contact, id, memory_id), + const ID & id) + : Regularisation(contact, id), mu(0, 1, 0., id + ":mu", 0., "mu") { AKANTU_DEBUG_IN(); Regularisation::registerSynchronizedArray(this->mu); this->registerParam("mu", this->mu, _pat_parsmod, "friction coefficient"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::computeFrictionalStrength() { AKANTU_DEBUG_IN(); // get contact arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); const SynchronizedArray & pressure = this->internalGetContactPressure(); // array to fill SynchronizedArray & strength = this->internalGetFrictionalStrength(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) strength(n) = 0.; // node pair is in contact else { // compute frictional strength strength(n) = this->mu(n) * pressure(n); } } Regularisation::computeFrictionalStrength(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->mu.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->mu.dumpRestartFile(file_name); Regularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->mu.readRestartFile(file_name); Regularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawCoulomb [" << std::endl; Regularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "mu") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->mu.getArray())); } /* else if (field_id == "frictional_contact_pressure") { this->internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::NodalField(this->frictional_contact_pressure.getArray())); } */ else { Regularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh index 5d86d0636..c67547a1e 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh @@ -1,114 +1,113 @@ /** * @file ntn_friclaw_linear_cohesive.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief linear cohesive law * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICLAW_LINEAR_COHESIVE_HH_ #define AST_NTN_FRICLAW_LINEAR_COHESIVE_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawLinearCohesive : public Regularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricLawLinearCohesive(NTNBaseContact & contact, - const ID & id = "linear_cohesive", - const MemoryID & memory_id = 0); + const ID & id = "linear_cohesive"); virtual ~NTNFricLawLinearCohesive(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // fracture energy SynchronizedArray G_c; // peak value of cohesive law SynchronizedArray tau_c; // residual value of cohesive law (for slip > d_c) SynchronizedArray tau_r; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const NTNFricLawLinearCohesive & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_linear_cohesive_tmpl.hh" #endif /* AST_NTN_FRICLAW_LINEAR_COHESIVE_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh index 887bc2e20..ba169ee22 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh @@ -1,188 +1,188 @@ /** * @file ntn_friclaw_linear_cohesive_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of linear cohesive law * * * Copyright (©) 2015-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 "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawLinearCohesive::NTNFricLawLinearCohesive( - NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) - : Regularisation(contact, id, memory_id), + NTNBaseContact & contact, const ID & id) + : Regularisation(contact, id), G_c(0, 1, 0., id + ":G_c", 0., "G_c"), tau_c(0, 1, 0., id + ":tau_c", 0., "tau_c"), tau_r(0, 1, 0., id + ":tau_r", 0., "tau_r") { AKANTU_DEBUG_IN(); Regularisation::registerSynchronizedArray(this->G_c); Regularisation::registerSynchronizedArray(this->tau_c); Regularisation::registerSynchronizedArray(this->tau_r); this->registerParam("G_c", this->G_c, _pat_parsmod, "fracture energy"); this->registerParam("tau_c", this->tau_c, _pat_parsmod, "peak shear strength"); this->registerParam("tau_r", this->tau_r, _pat_parsmod, "residual shear strength"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::computeFrictionalStrength() { AKANTU_DEBUG_IN(); // get arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); // const SynchronizedArray & slip = this->internalGetSlip(); const SynchronizedArray & slip = this->internalGetCumulativeSlip(); // array to fill SynchronizedArray & strength = this->internalGetFrictionalStrength(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) strength(n) = 0.; // node pair is in contact else { if (this->G_c(n) == 0.) { // strength(n) = 0.; strength(n) = this->tau_r(n); } else { Real slope = (this->tau_c(n) - this->tau_r(n)) * (this->tau_c(n) - this->tau_r(n)) / (2 * this->G_c(n)); // no strength < tau_r strength(n) = std::max(this->tau_c(n) - slope * slip(n), this->tau_r(n)); // strength(n) = std::max(this->tau_c(n) - slope * slip(n), 0.); // no // negative strength } } } Regularisation::computeFrictionalStrength(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->G_c.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->G_c.dumpRestartFile(file_name); this->tau_c.dumpRestartFile(file_name); this->tau_r.dumpRestartFile(file_name); Regularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->G_c.readRestartFile(file_name); this->tau_c.readRestartFile(file_name); this->tau_r.readRestartFile(file_name); Regularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawLinearCohesive [" << std::endl; Regularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "G_c") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->G_c.getArray())); } else if (field_id == "tau_c") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->tau_c.getArray())); } else if (field_id == "tau_r") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->tau_r.getArray())); } else { Regularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh index d12963771..60890c82c 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh @@ -1,116 +1,115 @@ /** * @file ntn_friclaw_linear_slip_weakening.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief linear slip weakening * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_HH_ #define AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_friclaw_coulomb.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawLinearSlipWeakening : public NTNFricLawCoulomb { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricLawLinearSlipWeakening(NTNBaseContact & contact, - const ID & id = "linear_slip_weakening", - const MemoryID & memory_id = 0); + const ID & id = "linear_slip_weakening"); virtual ~NTNFricLawLinearSlipWeakening(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /// computes the friction coefficient as a function of slip virtual void computeFrictionCoefficient(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // static coefficient of friction SynchronizedArray mu_s; // kinetic coefficient of friction SynchronizedArray mu_k; // Dc the length over which slip weakening happens SynchronizedArray d_c; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const NTNFricLawLinearSlipWeakening & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_linear_slip_weakening_tmpl.hh" #endif /* AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh index e6e44ba10..cc2c63f4f 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh @@ -1,95 +1,94 @@ /** * @file ntn_friclaw_linear_slip_weakening_no_healing.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief linear slip weakening * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_NO_HEALING_HH_ #define AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_NO_HEALING_HH_ /* -------------------------------------------------------------------------- */ #include "ntn_friclaw_linear_slip_weakening.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawLinearSlipWeakeningNoHealing : public NTNFricLawLinearSlipWeakening { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricLawLinearSlipWeakeningNoHealing( NTNBaseContact & contact, - const ID & id = "linear_slip_weakening_no_healing", - const MemoryID & memory_id = 0); + const ID & id = "linear_slip_weakening_no_healing"); virtual ~NTNFricLawLinearSlipWeakeningNoHealing(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// computes the friction coefficient as a function of slip virtual void computeFrictionCoefficient(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<( std::ostream & stream, const NTNFricLawLinearSlipWeakeningNoHealing & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh" #endif /* AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_NO_HEALING_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh index dee74841e..a8cc4cc41 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh @@ -1,85 +1,84 @@ /** * @file ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of linear slip weakening * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawLinearSlipWeakeningNoHealing:: NTNFricLawLinearSlipWeakeningNoHealing(NTNBaseContact & contact, - const ID & id, - const MemoryID & memory_id) - : NTNFricLawLinearSlipWeakening(contact, id, memory_id) { + const ID & id) + : NTNFricLawLinearSlipWeakening(contact, id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakeningNoHealing< Regularisation>::computeFrictionCoefficient() { AKANTU_DEBUG_IN(); // get arrays const SynchronizedArray & slip = this->internalGetCumulativeSlip(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (slip(n) >= this->d_c(n)) { this->mu(n) = this->mu_k(n); } else { // mu = mu_k + (1 - slip / Dc) * (mu_s - mu_k) this->mu(n) = this->mu_k(n) + (1 - (slip(n) / this->d_c(n))) * (this->mu_s(n) - this->mu_k(n)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakeningNoHealing::printself( std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawLinearSlipWeakeningNoHealing [" << std::endl; NTNFricLawLinearSlipWeakening::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh index 35c3d5d60..114143b36 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh @@ -1,189 +1,189 @@ /** * @file ntn_friclaw_linear_slip_weakening_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of linear slip weakening * * * Copyright (©) 2015-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 "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawLinearSlipWeakening::NTNFricLawLinearSlipWeakening( - NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) - : NTNFricLawCoulomb(contact, id, memory_id), + NTNBaseContact & contact, const ID & id) + : NTNFricLawCoulomb(contact, id), mu_s(0, 1, 0., id + ":mu_s", 0., "mu_s"), mu_k(0, 1, 0., id + ":mu_k", 0., "mu_k"), d_c(0, 1, 0., id + ":d_c", 0., "d_c") { AKANTU_DEBUG_IN(); NTNFricLawCoulomb::registerSynchronizedArray(this->mu_s); NTNFricLawCoulomb::registerSynchronizedArray(this->mu_k); NTNFricLawCoulomb::registerSynchronizedArray(this->d_c); this->registerParam("mu_s", this->mu_s, _pat_parsmod, "static friction coefficient"); this->registerParam("mu_k", this->mu_k, _pat_parsmod, "kinetic friction coefficient"); this->registerParam("d_c", this->d_c, _pat_parsmod, "slip weakening length"); this->setParameterAccessType("mu", _pat_readable); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening< Regularisation>::computeFrictionalStrength() { AKANTU_DEBUG_IN(); computeFrictionCoefficient(); NTNFricLawCoulomb::computeFrictionalStrength(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening< Regularisation>::computeFrictionCoefficient() { AKANTU_DEBUG_IN(); // get arrays const SynchronizedArray & stick = this->internalGetIsSticking(); const SynchronizedArray & slip = this->internalGetSlip(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (stick(n)) { this->mu(n) = this->mu_s(n); } else { if (slip(n) >= this->d_c(n)) { this->mu(n) = this->mu_k(n); } else { // mu = mu_k + (1 - slip / Dc) * (mu_s - mu_k) this->mu(n) = this->mu_k(n) + (1 - (slip(n) / this->d_c(n))) * (this->mu_s(n) - this->mu_k(n)); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->mu_s.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->mu_s.dumpRestartFile(file_name); this->mu_k.dumpRestartFile(file_name); this->d_c.dumpRestartFile(file_name); NTNFricLawCoulomb::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->mu_s.readRestartFile(file_name); this->mu_k.readRestartFile(file_name); this->d_c.readRestartFile(file_name); NTNFricLawCoulomb::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::printself( std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawLinearSlipWeakening [" << std::endl; NTNFricLawCoulomb::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "mu_s") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->mu_s.getArray())); } else if (field_id == "mu_k") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->mu_k.getArray())); } else if (field_id == "d_c") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->d_c.getArray())); } else { NTNFricLawCoulomb::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc index 85cab8044..76080511e 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc @@ -1,167 +1,167 @@ /** * @file ntn_fricreg_no_regularisation.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of no regularisation * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTNFricRegNoRegularisation::NTNFricRegNoRegularisation( - NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) - : NTNBaseFriction(contact, id, memory_id), + NTNBaseContact & contact, const ID & id) + : NTNBaseFriction(contact, id), frictional_contact_pressure(0, 1, 0., id + ":frictional_contact_pressure", 0., "frictional_contact_pressure") { AKANTU_DEBUG_IN(); NTNBaseFriction::registerSynchronizedArray(this->frictional_contact_pressure); this->registerParam("frictional_contact_pressure", this->frictional_contact_pressure, _pat_internal, "contact pressure used for friction law"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const SynchronizedArray & NTNFricRegNoRegularisation::internalGetContactPressure() { AKANTU_DEBUG_IN(); this->computeFrictionalContactPressure(); AKANTU_DEBUG_OUT(); return this->frictional_contact_pressure; } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::computeFrictionalContactPressure() { AKANTU_DEBUG_IN(); SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); // get contact arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); const Array & pressure = this->contact.getContactPressure().getArray(); Array::const_iterator> it = pressure.begin(dim); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) this->frictional_contact_pressure(n) = 0.; // node pair is in contact else { // compute frictional contact pressure const Vector & pres = it[n]; this->frictional_contact_pressure(n) = pres.norm(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->frictional_contact_pressure.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->frictional_contact_pressure.dumpRestartFile(file_name); NTNBaseFriction::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->frictional_contact_pressure.readRestartFile(file_name); NTNBaseFriction::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricRegNoRegularisation [" << std::endl; NTNBaseFriction::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "frictional_contact_pressure") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>( this->frictional_contact_pressure.getArray())); } else { NTNBaseFriction::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh index f7df9f108..75f399819 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh @@ -1,133 +1,132 @@ /** * @file ntn_fricreg_no_regularisation.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief regularisation that does nothing * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICREG_NO_REGULARISATION_HH_ #define AST_NTN_FRICREG_NO_REGULARISATION_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_friction.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNFricRegNoRegularisation : public NTNBaseFriction { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricRegNoRegularisation(NTNBaseContact & contact, - const ID & id = "no_regularisation", - const MemoryID & memory_id = 0); + const ID & id = "no_regularisation"); virtual ~NTNFricRegNoRegularisation(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set to steady state for no regularisation -> do nothing virtual void setToSteadyState(){}; virtual void registerSynchronizedArray(SynchronizedArrayBase & array); virtual void dumpRestart(const std::string & file_name) const; virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: virtual void computeFrictionalContactPressure(); /// compute frictional strength according to friction law virtual void computeFrictionalStrength(){}; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// get the is_in_contact array virtual const SynchronizedArray & internalGetIsInContact() { return this->contact.getIsInContact(); }; /// get the contact pressure (the norm: scalar value) virtual const SynchronizedArray & internalGetContactPressure(); /// get the frictional strength array virtual SynchronizedArray & internalGetFrictionalStrength() { return this->frictional_strength; }; /// get the is_sticking array virtual SynchronizedArray & internalGetIsSticking() { return this->is_sticking; } /// get the slip array virtual SynchronizedArray & internalGetSlip() { return this->slip; } /// get the slip array virtual SynchronizedArray & internalGetCumulativeSlip() { return this->cumulative_slip; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // contact pressure (absolut value) for computation of friction SynchronizedArray frictional_contact_pressure; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_fricreg_no_regularisation_inline_impl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNFricRegNoRegularisation & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_NTN_FRICREG_NO_REGULARISATION_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc index dfecefbb4..6877a5bf0 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc @@ -1,175 +1,174 @@ /** * @file ntn_fricreg_rubin_ampuero.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of no regularisation * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_rubin_ampuero.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTNFricRegRubinAmpuero::NTNFricRegRubinAmpuero(NTNBaseContact & contact, - const ID & id, - const MemoryID & memory_id) - : NTNFricRegNoRegularisation(contact, id, memory_id), + const ID & id) + : NTNFricRegNoRegularisation(contact, id), t_star(0, 1, 0., id + ":t_star", 0., "t_star") { AKANTU_DEBUG_IN(); NTNFricRegNoRegularisation::registerSynchronizedArray(this->t_star); this->registerParam("t_star", this->t_star, _pat_parsmod, "time scale of regularization"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const SynchronizedArray & NTNFricRegRubinAmpuero::internalGetContactPressure() { AKANTU_DEBUG_IN(); SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); Real delta_t = model.getTimeStep(); // get contact arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); const Array & pressure = this->contact.getContactPressure().getArray(); Array::const_iterator> it = pressure.begin(dim); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) this->frictional_contact_pressure(n) = 0.; // if t_star is too small compute like Coulomb friction (without // regularization) else if (Math::are_float_equal(this->t_star(n), 0.)) { const Vector & pres = it[n]; this->frictional_contact_pressure(n) = pres.norm(); } else { // compute frictional contact pressure // backward euler method: first order implicit numerical integration // method // \reg_pres_n+1 = (\reg_pres_n + \delta_t / \t_star * \cur_pres) // / (1 + \delta_t / \t_star) Real alpha = delta_t / this->t_star(n); const Vector & pres = it[n]; this->frictional_contact_pressure(n) += alpha * pres.norm(); this->frictional_contact_pressure(n) /= 1 + alpha; } } AKANTU_DEBUG_OUT(); return this->frictional_contact_pressure; } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::setToSteadyState() { AKANTU_DEBUG_IN(); NTNFricRegNoRegularisation::computeFrictionalContactPressure(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->t_star.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->t_star.dumpRestartFile(file_name); NTNFricRegNoRegularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->t_star.readRestartFile(file_name); NTNFricRegNoRegularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricRegRubinAmpuero [" << std::endl; NTNFricRegNoRegularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "t_star") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->t_star.getArray())); } else { NTNFricRegNoRegularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh index 71593c4b8..2a6b85580 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh @@ -1,101 +1,100 @@ /** * @file ntn_fricreg_rubin_ampuero.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief regularisation that regularizes the contact pressure * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICREG_RUBIN_AMPUERO_HH_ #define AST_NTN_FRICREG_RUBIN_AMPUERO_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNFricRegRubinAmpuero : public NTNFricRegNoRegularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricRegRubinAmpuero(NTNBaseContact & contact, - const ID & id = "rubin_ampuero", - const MemoryID & memory_id = 0); + const ID & id = "rubin_ampuero"); virtual ~NTNFricRegRubinAmpuero(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void registerSynchronizedArray(SynchronizedArrayBase & array); virtual void dumpRestart(const std::string & file_name) const; virtual void readRestart(const std::string & file_name); virtual void setToSteadyState(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// get the contact pressure (the norm: scalar value) virtual const SynchronizedArray & internalGetContactPressure(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: SynchronizedArray t_star; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_fricreg_rubin_ampuero_inline_impl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNFricRegRubinAmpuero & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_NTN_FRICREG_RUBIN_AMPUERO_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc index 5812f7779..dfacc1585 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc @@ -1,162 +1,162 @@ /** * @file ntn_fricreg_simplified_prakash_clifton.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of simplified prakash clifton with one parameter * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_simplified_prakash_clifton.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTNFricRegSimplifiedPrakashClifton::NTNFricRegSimplifiedPrakashClifton( - NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) - : NTNFricRegNoRegularisation(contact, id, memory_id), + NTNBaseContact & contact, const ID & id) + : NTNFricRegNoRegularisation(contact, id), t_star(0, 1, 0., id + ":t_star", 0., "t_star"), spc_internal(0, 1, 0., id + ":spc_internal", 0., "spc_internal") { AKANTU_DEBUG_IN(); NTNFricRegNoRegularisation::registerSynchronizedArray(this->t_star); NTNFricRegNoRegularisation::registerSynchronizedArray(this->spc_internal); this->registerParam("t_star", this->t_star, _pat_parsmod, "time scale of regularisation"); this->registerParam("spc_internal", this->spc_internal, _pat_internal, ""); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::computeFrictionalStrength() { AKANTU_DEBUG_IN(); SolidMechanicsModel & model = this->contact.getModel(); Real delta_t = model.getTimeStep(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { Real alpha = delta_t / this->t_star(n); this->frictional_strength(n) += alpha * this->spc_internal(n); this->frictional_strength(n) /= 1 + alpha; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::setToSteadyState() { AKANTU_DEBUG_IN(); /// fill the spc_internal array computeFrictionalStrength(); /// set strength without regularisation UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { this->frictional_strength(n) = this->spc_internal(n); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->t_star.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->t_star.dumpRestartFile(file_name); this->spc_internal.dumpRestartFile(file_name); NTNFricRegNoRegularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->t_star.readRestartFile(file_name); this->spc_internal.readRestartFile(file_name); NTNFricRegNoRegularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricRegSimplifiedPrakashClifton [" << std::endl; NTNFricRegNoRegularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "t_star") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->t_star.getArray())); } else { NTNFricRegNoRegularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh index c513fdeca..e8d285288 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh @@ -1,112 +1,111 @@ /** * @file ntn_fricreg_simplified_prakash_clifton.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief regularisation that regularizes the frictional strength with one * parameter * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICREG_SIMPLIFIED_PRAKASH_CLIFTON_HH_ #define AST_NTN_FRICREG_SIMPLIFIED_PRAKASH_CLIFTON_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNFricRegSimplifiedPrakashClifton : public NTNFricRegNoRegularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricRegSimplifiedPrakashClifton( - NTNBaseContact & contact, const ID & id = "simplified_prakash_clifton", - const MemoryID & memory_id = 0); + NTNBaseContact & contact, const ID & id = "simplified_prakash_clifton"); virtual ~NTNFricRegSimplifiedPrakashClifton(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void registerSynchronizedArray(SynchronizedArrayBase & array); virtual void dumpRestart(const std::string & file_name) const; virtual void readRestart(const std::string & file_name); virtual void setToSteadyState(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// get the frictional strength array virtual SynchronizedArray & internalGetFrictionalStrength() { return this->spc_internal; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: SynchronizedArray t_star; // to get the incremental frictional strength SynchronizedArray spc_internal; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_fricreg_simplified_prakash_clifton_inline_impl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNFricRegSimplifiedPrakashClifton & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_NTN_FRICREG_SIMPLIFIED_PRAKASH_CLIFTON_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc index 5e221e535..014b71886 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc @@ -1,119 +1,118 @@ /** * @file mIIasym_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "mIIasym_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -MIIASYMContact::MIIASYMContact(SolidMechanicsModel & model, const ID & id, - const MemoryID & memory_id) - : NTRFContact(model, id, memory_id) { +MIIASYMContact::MIIASYMContact(SolidMechanicsModel & model, const ID & id) + : NTRFContact(model, id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::updateImpedance() { AKANTU_DEBUG_IN(); NTRFContact::updateImpedance(); for (UInt i = 0; i < this->impedance.size(); ++i) { this->impedance(i) *= 0.5; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// WARNING: this is only valid for the acceleration in equilibrium void MIIASYMContact::computeRelativeNormalField( const Array & field, Array & rel_normal_field) const { AKANTU_DEBUG_IN(); NTRFContact::computeRelativeNormalField(field, rel_normal_field); for (auto it_rtfield = rel_normal_field.begin(); it_rtfield != rel_normal_field.end(); ++it_rtfield) { // in the anti-symmetric case *it_rtfield *= 2.; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::computeRelativeTangentialField( const Array & field, Array & rel_tang_field) const { AKANTU_DEBUG_IN(); NTRFContact::computeRelativeTangentialField(field, rel_tang_field); UInt dim = this->model.getSpatialDimension(); for (Array::iterator> it_rtfield = rel_tang_field.begin(dim); it_rtfield != rel_tang_field.end(dim); ++it_rtfield) { // in the anti-symmetric case, the tangential fields become twice as large *it_rtfield *= 2.; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::computeContactPressureInEquilibrium() { AKANTU_DEBUG_IN(); NTRFContact::computeContactPressure(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "MIIASYMContact [" << std::endl; NTRFContact::printself(stream, indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh index 3d2855cc6..5de74ea4d 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh @@ -1,90 +1,90 @@ /** * @file mIIasym_contact.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief contact for mode II anti-symmetric simulations * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_MIIASYM_CONTACT_HH_ #define AST_MIIASYM_CONTACT_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntrf_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class MIIASYMContact : public NTRFContact { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - MIIASYMContact(SolidMechanicsModel & model, const ID & id = "contact", - const MemoryID & memory_id = 0); - virtual ~MIIASYMContact() = default; + MIIASYMContact(SolidMechanicsModel & model, const ID & id = "contact"); + ~MIIASYMContact() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// update the impedance matrix virtual void updateImpedance(); /// compute contact pressure -> do nothing because can only compute it in /// equilibrium - virtual void computeContactPressure(){}; + void computeContactPressure() override{}; /// compute relative normal field (only value that has to be multiplied with /// the normal) /// WARNING: this is only valid for the acceleration in equilibrium - virtual void computeRelativeNormalField(const Array & field, - Array & rel_normal_field) const; + void + computeRelativeNormalField(const Array & field, + Array & rel_normal_field) const override; /// compute relative tangential field (complet array) /// relative to master nodes - virtual void + void computeRelativeTangentialField(const Array & field, - Array & rel_tang_field) const; + Array & rel_tang_field) const override; /// compute contact pressure that is used over the entire time virtual void computeContactPressureInEquilibrium(); /// function to print the contain of the class - virtual void printself(std::ostream & stream, int indent = 0) const; + void printself(std::ostream & stream, int indent = 0) const override; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const MIIASYMContact & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_MIIASYM_CONTACT_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc index 4f9c385a6..b2a7c43c2 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc @@ -1,567 +1,565 @@ /** * @file ntn_base_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn base contact * * * Copyright (©) 2015-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 "ntn_base_contact.hh" #include "dof_manager_default.hh" #include "dumpable_inline_impl.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" #include "element_synchronizer.hh" #include "mesh_utils.hh" #include "non_linear_solver_lumped.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ // NTNContactSynchElementFilter::NTNContactSynchElementFilter( // NTNBaseContact & contact) // : contact(contact), // connectivity(contact.getModel().getMesh().getConnectivities()) { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // bool NTNContactSynchElementFilter::operator()(const Element & e) { // AKANTU_DEBUG_IN(); // ElementType type = e.type; // UInt element = e.element; // GhostType ghost_type = e.ghost_type; // // loop over all nodes of this element // bool need_element = false; // UInt nb_nodes = Mesh::getNbNodesPerElement(type); // for (UInt n = 0; n < nb_nodes; ++n) { // UInt nn = this->connectivity(type, ghost_type)(element, n); // // if one nodes is in this contact, we need this element // if (this->contact.getNodeIndex(nn) >= 0) { // need_element = true; // break; // } // } // AKANTU_DEBUG_OUT(); // return need_element; // } /* -------------------------------------------------------------------------- */ -NTNBaseContact::NTNBaseContact(SolidMechanicsModel & model, const ID & id, - const MemoryID & memory_id) - : Memory(id, memory_id), Dumpable(), model(model), +NTNBaseContact::NTNBaseContact(SolidMechanicsModel & model, const ID & id) + : id(id), model(model), slaves(0, 1, 0, id + ":slaves", std::numeric_limits::quiet_NaN(), "slaves"), normals(0, model.getSpatialDimension(), 0, id + ":normals", std::numeric_limits::quiet_NaN(), "normals"), contact_pressure( 0, model.getSpatialDimension(), 0, id + ":contact_pressure", std::numeric_limits::quiet_NaN(), "contact_pressure"), is_in_contact(0, 1, false, id + ":is_in_contact", false, "is_in_contact"), lumped_boundary_slaves(0, 1, 0, id + ":lumped_boundary_slaves", std::numeric_limits::quiet_NaN(), "lumped_boundary_slaves"), impedance(0, 1, 0, id + ":impedance", std::numeric_limits::quiet_NaN(), "impedance"), - contact_surfaces(), slave_elements("slave_elements", id, memory_id), - node_to_elements() { + slave_elements("slave_elements", id) { AKANTU_DEBUG_IN(); auto & boundary_fem = this->model.getFEEngineBoundary(); for (auto && ghost_type : ghost_types) { boundary_fem.initShapeFunctions(ghost_type); } auto & mesh = this->model.getMesh(); auto spatial_dimension = this->model.getSpatialDimension(); this->slave_elements.initialize(mesh, _spatial_dimension = spatial_dimension - 1); MeshUtils::buildNode2Elements(mesh, this->node_to_elements, spatial_dimension - 1); this->registerDumper("text_all", id, true); this->addDumpFilteredMesh(mesh, slave_elements, slaves.getArray(), spatial_dimension - 1, _not_ghost, _ek_regular); // parallelisation this->synch_registry = std::make_unique(); this->synch_registry->registerDataAccessor(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NTNBaseContact::~NTNBaseContact() = default; /* -------------------------------------------------------------------------- */ void NTNBaseContact::initParallel() { AKANTU_DEBUG_IN(); this->synchronizer = std::make_unique( this->model.getMesh().getElementSynchronizer()); this->synchronizer->filterScheme([&](auto && element) { // loop over all nodes of this element Vector conn = const_cast(this->model.getMesh()) .getConnectivity(element); for (auto & node : conn) { // if one nodes is in this contact, we need this element if (this->getNodeIndex(node) >= 0) { return true; } } return false; }); this->synch_registry->registerSynchronizer(*(this->synchronizer), SynchronizationTag::_cf_nodal); this->synch_registry->registerSynchronizer(*(this->synchronizer), SynchronizationTag::_cf_incr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::findBoundaryElements( const Array & interface_nodes, ElementTypeMapArray & elements) { AKANTU_DEBUG_IN(); // add connected boundary elements that have all nodes on this contact for (const auto & node : interface_nodes) { for (const auto & element : this->node_to_elements.getRow(node)) { Vector conn = const_cast(this->model.getMesh()) .getConnectivity(element); auto nb_nodes = conn.size(); decltype(nb_nodes) nb_found_nodes = 0; for (auto & nn : conn) { if (interface_nodes.find(nn) != UInt(-1)) { nb_found_nodes++; } else { break; } } // this is an element between all contact nodes // and is not already in the elements if ((nb_found_nodes == nb_nodes) && (elements(element.type, element.ghost_type).find(element.element) == UInt(-1))) { elements(element.type, element.ghost_type).push_back(element.element); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::addSplitNode(UInt node, UInt) { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); // add to node arrays this->slaves.push_back(node); // set contact as false this->is_in_contact.push_back(false); // before initializing // set contact pressure, normal, lumped_boundary to Nan this->contact_pressure.push_back(std::numeric_limits::quiet_NaN()); this->impedance.push_back(std::numeric_limits::quiet_NaN()); this->lumped_boundary_slaves.push_back( std::numeric_limits::quiet_NaN()); Vector nan_normal(dim, std::numeric_limits::quiet_NaN()); this->normals.push_back(nan_normal); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::registerSynchronizedArray(SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->slaves.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->slaves.dumpRestartFile(file_name); this->normals.dumpRestartFile(file_name); this->is_in_contact.dumpRestartFile(file_name); this->contact_pressure.dumpRestartFile(file_name); this->lumped_boundary_slaves.dumpRestartFile(file_name); this->impedance.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->slaves.readRestartFile(file_name); this->normals.readRestartFile(file_name); this->is_in_contact.readRestartFile(file_name); this->contact_pressure.readRestartFile(file_name); this->lumped_boundary_slaves.readRestartFile(file_name); this->impedance.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt NTNBaseContact::getNbNodesInContact() const { AKANTU_DEBUG_IN(); UInt nb_contact = 0; UInt nb_nodes = this->getNbContactNodes(); const Mesh & mesh = this->model.getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(this->slaves(n)); bool is_pbc_slave_node = mesh.isPeriodicSlave(this->slaves(n)); if (is_local_node && !is_pbc_slave_node && this->is_in_contact(n)) { nb_contact++; } } mesh.getCommunicator().allReduce(nb_contact, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_contact; } /* -------------------------------------------------------------------------- */ void NTNBaseContact::updateInternalData() { AKANTU_DEBUG_IN(); updateNormals(); updateLumpedBoundary(); updateImpedance(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::updateLumpedBoundary() { AKANTU_DEBUG_IN(); this->internalUpdateLumpedBoundary(this->slaves.getArray(), this->slave_elements, this->lumped_boundary_slaves); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::internalUpdateLumpedBoundary( const Array & nodes, const ElementTypeMapArray & elements, SynchronizedArray & boundary) { AKANTU_DEBUG_IN(); // set all values in lumped_boundary to zero boundary.zero(); UInt dim = this->model.getSpatialDimension(); // UInt nb_contact_nodes = getNbContactNodes(); const FEEngine & boundary_fem = this->model.getFEEngineBoundary(); const Mesh & mesh = this->model.getMesh(); for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(dim - 1, ghost_type)) { UInt nb_elements = mesh.getNbElement(type, ghost_type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, ghost_type); // get shapes and compute integral const Array & shapes = boundary_fem.getShapes(type, ghost_type); Array area(nb_elements, nb_nodes_per_element); boundary_fem.integrate(shapes, area, nb_nodes_per_element, type, ghost_type); if (this->contact_surfaces.size() == 0) { AKANTU_DEBUG_WARNING( "No surfaces in ntn base contact." << " You have to define the lumped boundary by yourself."); } Array::const_iterator elem_it = (elements)(type, ghost_type).begin(); Array::const_iterator elem_it_end = (elements)(type, ghost_type).end(); // loop over contact nodes for (; elem_it != elem_it_end; ++elem_it) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { UInt node = connectivity(*elem_it, q); UInt node_index = nodes.find(node); AKANTU_DEBUG_ASSERT(node_index != UInt(-1), "Could not find node " << node << " in the array!"); Real area_to_add = area(*elem_it, q); boundary(node_index) += area_to_add; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::computeAcceleration(Array & acceleration) const { auto && dof_manager = dynamic_cast(model.getDOFManager()); const auto & b = dof_manager.getResidual(); acceleration.resize(b.size()); const auto & blocked_dofs = dof_manager.getGlobalBlockedDOFs(); const auto & A = dof_manager.getLumpedMatrix("M"); Array blocked_dofs_bool(blocked_dofs.size()); for (auto && data : zip(blocked_dofs, blocked_dofs_bool)) { std::get<1>(data) = std::get<0>(data); } // pre-compute the acceleration // (not increment acceleration, because residual is still Kf) NonLinearSolverLumped::solveLumped(A, acceleration, b, this->model.getF_M2A(), blocked_dofs_bool); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::computeContactPressure() { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); Real delta_t = this->model.getTimeStep(); UInt nb_contact_nodes = getNbContactNodes(); AKANTU_DEBUG_ASSERT(delta_t > 0., "Cannot compute contact pressure if no time step is set"); // synchronize data this->synch_registry->synchronize(SynchronizationTag::_cf_nodal); Array acceleration(0, dim); this->computeAcceleration(acceleration); // compute relative normal fields of displacement, velocity and acceleration Array r_disp(0, 1); Array r_velo(0, 1); Array r_acce(0, 1); Array r_old_acce(0, 1); computeNormalGap(r_disp); // computeRelativeNormalField(this->model.getCurrentPosition(), r_disp); computeRelativeNormalField(this->model.getVelocity(), r_velo); computeRelativeNormalField(acceleration, r_acce); computeRelativeNormalField(this->model.getAcceleration(), r_old_acce); AKANTU_DEBUG_ASSERT(r_disp.size() == nb_contact_nodes, "computeRelativeNormalField does not give back arrays " << "size == nb_contact_nodes. nb_contact_nodes = " << nb_contact_nodes << " | array size = " << r_disp.size()); // compute gap array for all nodes Array gap(nb_contact_nodes, 1); Real * gap_p = gap.storage(); Real * r_disp_p = r_disp.storage(); Real * r_velo_p = r_velo.storage(); Real * r_acce_p = r_acce.storage(); Real * r_old_acce_p = r_old_acce.storage(); for (UInt i = 0; i < nb_contact_nodes; ++i) { *gap_p = *r_disp_p + delta_t * *r_velo_p + delta_t * delta_t * *r_acce_p - 0.5 * delta_t * delta_t * *r_old_acce_p; // increment pointers gap_p++; r_disp_p++; r_velo_p++; r_acce_p++; r_old_acce_p++; } // check if gap is negative -> is in contact for (UInt n = 0; n < nb_contact_nodes; ++n) { if (gap(n) <= 0.) { for (UInt d = 0; d < dim; ++d) { this->contact_pressure(n, d) = this->impedance(n) * gap(n) / (2 * delta_t) * this->normals(n, d); } this->is_in_contact(n) = true; } else { for (UInt d = 0; d < dim; ++d) this->contact_pressure(n, d) = 0.; this->is_in_contact(n) = false; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::applyContactPressure() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); UInt dim = this->model.getSpatialDimension(); Array & residual = this->model.getInternalForce(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt slave = this->slaves(n); for (UInt d = 0; d < dim; ++d) { // residual(master,d) += this->lumped_boundary(n,0) * // this->contact_pressure(n,d); residual(slave, d) -= this->lumped_boundary_slaves(n) * this->contact_pressure(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int NTNBaseContact::getNodeIndex(UInt node) const { return this->slaves.find(node); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNBaseContact [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + slaves : " << std::endl; this->slaves.printself(stream, indent + 2); stream << space << " + normals : " << std::endl; this->normals.printself(stream, indent + 2); stream << space << " + contact_pressure : " << std::endl; this->contact_pressure.printself(stream, indent + 2); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::syncArrays(SyncChoice sync_choice) { AKANTU_DEBUG_IN(); this->slaves.syncElements(sync_choice); this->normals.syncElements(sync_choice); this->is_in_contact.syncElements(sync_choice); this->contact_pressure.syncElements(sync_choice); this->lumped_boundary_slaves.syncElements(sync_choice); this->impedance.syncElements(sync_choice); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER const Array & nodal_filter = this->slaves.getArray(); #define ADD_FIELD(field_id, field, type) \ internalAddDumpFieldToDumper( \ dumper_name, field_id, \ std::make_unique< \ - dumpers::NodalField, Array>>( \ + dumpers::NodalField, Array>>( \ field, 0, 0, &nodal_filter)) if (field_id == "displacement") { ADD_FIELD(field_id, this->model.getDisplacement(), Real); } else if (field_id == "mass") { ADD_FIELD(field_id, this->model.getMass(), Real); } else if (field_id == "velocity") { ADD_FIELD(field_id, this->model.getVelocity(), Real); } else if (field_id == "acceleration") { ADD_FIELD(field_id, this->model.getAcceleration(), Real); } else if (field_id == "external_force") { ADD_FIELD(field_id, this->model.getExternalForce(), Real); } else if (field_id == "internal_force") { ADD_FIELD(field_id, this->model.getInternalForce(), Real); } else if (field_id == "blocked_dofs") { ADD_FIELD(field_id, this->model.getBlockedDOFs(), bool); } else if (field_id == "increment") { ADD_FIELD(field_id, this->model.getIncrement(), Real); } else if (field_id == "normal") { internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->normals.getArray())); } else if (field_id == "contact_pressure") { internalAddDumpFieldToDumper(dumper_name, field_id, std::make_unique>( this->contact_pressure.getArray())); } else if (field_id == "is_in_contact") { internalAddDumpFieldToDumper(dumper_name, field_id, std::make_unique>( this->is_in_contact.getArray())); } else if (field_id == "lumped_boundary_slave") { internalAddDumpFieldToDumper(dumper_name, field_id, std::make_unique>( this->lumped_boundary_slaves.getArray())); } else if (field_id == "impedance") { - internalAddDumpFieldToDumper( - dumper_name, field_id, - std::make_unique>(this->impedance.getArray())); + internalAddDumpFieldToDumper(dumper_name, field_id, + std::make_unique>( + this->impedance.getArray())); } else { std::cerr << "Could not add field '" << field_id << "' to the dumper. Just ignored it." << std::endl; } #undef ADD_FIELD #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh index 8c9c3cf9d..d786ccab5 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh @@ -1,251 +1,250 @@ /** * @file ntn_base_contact.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief base contact for ntn and ntrf contact * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_BASE_CONTACT_HH_ #define AST_NTN_BASE_CONTACT_HH_ /* -------------------------------------------------------------------------- */ // akantu #include "aka_csr.hh" #include "solid_mechanics_model.hh" // simtools #include "synchronized_array.hh" namespace akantu { class NTNBaseContact; /* -------------------------------------------------------------------------- */ // class NTNContactSynchElementFilter : public SynchElementFilter { // public: // // constructor // NTNContactSynchElementFilter(NTNBaseContact & contact); // // answer to: do we need this element ? // virtual bool operator()(const Element & e); // private: // const NTNBaseContact & contact; // const ElementTypeMapArray & connectivity; // }; /* -------------------------------------------------------------------------- */ -class NTNBaseContact : protected Memory, - public DataAccessor, - public Dumpable { +class NTNBaseContact : public DataAccessor, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNBaseContact(SolidMechanicsModel & model, const ID & id = "contact", - const MemoryID & memory_id = 0); - virtual ~NTNBaseContact(); + NTNBaseContact(SolidMechanicsModel & model, const ID & id = "contact"); + ~NTNBaseContact() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initializes ntn contact parallel virtual void initParallel(); /// add split node virtual void addSplitNode(UInt node, UInt = 0); /// update normals, lumped boundary, and impedance virtual void updateInternalData(); /// update (compute the normals) virtual void updateNormals() = 0; /// update the lumped boundary B matrix virtual void updateLumpedBoundary(); /// update the impedance matrix virtual void updateImpedance() = 0; /// compute the normal contact force virtual void computeContactPressure(); /// impose the normal contact force - virtual void applyContactPressure() ; + virtual void applyContactPressure(); /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// compute the normal gap virtual void computeNormalGap(Array & gap) const = 0; /// compute relative normal field (only value that has to be multiplied with /// the normal) /// relative to master nodes virtual void computeRelativeNormalField(const Array & field, Array & rel_normal_field) const = 0; /// compute relative tangential field (complet array) /// relative to master nodes virtual void computeRelativeTangentialField(const Array & field, Array & rel_tang_field) const = 0; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// computes the acceleration void computeAcceleration(Array & acceleration) const; protected: /// updateLumpedBoundary virtual void internalUpdateLumpedBoundary(const Array & nodes, const ElementTypeMapArray & elements, SynchronizedArray & boundary); // to find the slave_elements or master_elements virtual void findBoundaryElements(const Array & interface_nodes, ElementTypeMapArray & elements); /// synchronize arrays virtual void syncArrays(SyncChoice sync_choice); /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: void addDumpFieldToDumper(const std::string & dumper_name, - const std::string & field_id) override; + const std::string & field_id) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, model, SolidMechanicsModel &) AKANTU_GET_MACRO(Slaves, slaves, const SynchronizedArray &) AKANTU_GET_MACRO(Normals, normals, const SynchronizedArray &) AKANTU_GET_MACRO(ContactPressure, contact_pressure, const SynchronizedArray &) AKANTU_GET_MACRO(LumpedBoundarySlaves, lumped_boundary_slaves, const SynchronizedArray &) AKANTU_GET_MACRO(Impedance, impedance, const SynchronizedArray &) AKANTU_GET_MACRO(IsInContact, is_in_contact, const SynchronizedArray &) AKANTU_GET_MACRO(SlaveElements, slave_elements, const ElementTypeMapArray &) AKANTU_GET_MACRO(SynchronizerRegistry, *synch_registry, SynchronizerRegistry &) /// get number of nodes that are in contact (globally, on all procs together) /// is_in_contact = true virtual UInt getNbNodesInContact() const; /// get index of node in either slaves or masters array /// if node is in neither of them, return -1 virtual Int getNodeIndex(UInt node) const; /// get number of contact nodes: nodes in the system locally (on this proc) /// is_in_contact = true and false, because just in the system virtual UInt getNbContactNodes() const { return this->slaves.size(); } bool isNTNContact() const { return this->is_ntn_contact; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - typedef std::set SurfacePtrSet; + using SurfacePtrSet = std::set; + + ID id; SolidMechanicsModel & model; /// array of slave nodes SynchronizedArray slaves; /// array of normals SynchronizedArray normals; /// array indicating if nodes are in contact SynchronizedArray contact_pressure; /// array indicating if nodes are in contact SynchronizedArray is_in_contact; /// boundary matrix for slave nodes SynchronizedArray lumped_boundary_slaves; /// impedance matrix SynchronizedArray impedance; /// contact surface SurfacePtrSet contact_surfaces; /// element list for dump and lumped_boundary ElementTypeMapArray slave_elements; CSR node_to_elements; /// parallelisation std::unique_ptr synch_registry; std::unique_ptr synchronizer; bool is_ntn_contact{true}; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNBaseContact & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "ntn_base_contact_inline_impl.hh" #endif /* AST_NTN_BASE_CONTACT_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc index efa895782..e28e034fe 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc @@ -1,380 +1,379 @@ /** * @file ntn_base_friction.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn base friction * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_friction.hh" #include "dof_manager_default.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" #include "non_linear_solver_lumped.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -NTNBaseFriction::NTNBaseFriction(NTNBaseContact & contact, const ID & id, - const MemoryID & memory_id) - : Memory(id, memory_id), Parsable(ParserType::_friction, id), Dumpable(), +NTNBaseFriction::NTNBaseFriction(NTNBaseContact & contact, const ID & id) + : Parsable(ParserType::_friction, id), Dumpable(), contact(contact), is_sticking(0, 1, true, id + ":is_sticking", true, "is_sticking"), frictional_strength(0, 1, 0., id + ":frictional_strength", 0., "frictional_strength"), friction_traction(0, contact.getModel().getSpatialDimension(), 0., id + ":friction_traction", 0., "friction_traction"), slip(0, 1, 0., id + ":slip", 0., "slip"), cumulative_slip(0, 1, 0., id + ":cumulative_slip", 0., "cumulative_slip"), slip_velocity(0, contact.getModel().getSpatialDimension(), 0., id + ":slip_velocity", 0., "slip_velocity") { AKANTU_DEBUG_IN(); this->contact.registerSynchronizedArray(this->is_sticking); this->contact.registerSynchronizedArray(this->frictional_strength); this->contact.registerSynchronizedArray(this->friction_traction); this->contact.registerSynchronizedArray(this->slip); this->contact.registerSynchronizedArray(this->cumulative_slip); this->contact.registerSynchronizedArray(this->slip_velocity); this->registerExternalDumper(contact.getDumper().shared_from_this(), contact.getDefaultDumperName(), true); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::updateSlip() { AKANTU_DEBUG_IN(); SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); // synchronize increment this->contact.getSynchronizerRegistry().synchronize( SynchronizationTag::_cf_incr); Array rel_tan_incr(0, dim); this->contact.computeRelativeTangentialField(model.getIncrement(), rel_tan_incr); Array::const_iterator> it = rel_tan_incr.begin(dim); UInt nb_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->is_sticking(n)) { this->slip(n) = 0.; } else { const Vector & rti = it[n]; this->slip(n) += rti.norm(); this->cumulative_slip(n) += rti.norm(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::computeFrictionTraction() { AKANTU_DEBUG_IN(); this->computeStickTraction(); this->computeFrictionalStrength(); SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); // get contact arrays const SynchronizedArray & is_in_contact = this->contact.getIsInContact(); Array & traction = const_cast &>(this->friction_traction.getArray()); Array::iterator> it_fric_trac = traction.begin(dim); this->is_sticking.zero(); // set to not sticking UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is in contact if (is_in_contact(n)) { Vector fric_trac = it_fric_trac[n]; // check if it is larger than frictional strength Real abs_fric = fric_trac.norm(); if (abs_fric != 0.) { Real alpha = this->frictional_strength(n) / abs_fric; // larger -> sliding if (alpha < 1.) { fric_trac *= alpha; } else this->is_sticking(n) = true; } else { // frictional traction is already zero this->is_sticking(n) = true; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::computeStickTraction() { AKANTU_DEBUG_IN(); SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); Real delta_t = model.getTimeStep(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); // get contact arrays const SynchronizedArray & impedance = this->contact.getImpedance(); const SynchronizedArray & is_in_contact = this->contact.getIsInContact(); Array acceleration(0, dim); this->contact.computeAcceleration(acceleration); // compute relative normal fields of velocity and acceleration Array r_velo(0, dim); Array r_acce(0, dim); Array r_old_acce(0, dim); this->contact.computeRelativeTangentialField(model.getVelocity(), r_velo); this->contact.computeRelativeTangentialField(acceleration, r_acce); this->contact.computeRelativeTangentialField(model.getAcceleration(), r_old_acce); AKANTU_DEBUG_ASSERT(r_velo.size() == nb_contact_nodes, "computeRelativeNormalField does not give back arrays " << "size == nb_contact_nodes. nb_contact_nodes = " << nb_contact_nodes << " | array size = " << r_velo.size()); // compute tangential gap_dot array for all nodes Array gap_dot(nb_contact_nodes, dim); for (auto && data : zip(make_view(gap_dot), make_view(r_velo), make_view(r_acce), make_view(r_old_acce))) { auto & gap_dot = std::get<0>(data); auto & r_velo = std::get<1>(data); auto & r_acce = std::get<2>(data); auto & r_old_acce = std::get<3>(data); gap_dot = r_velo + delta_t * r_acce - 1. / 2. * delta_t * r_old_acce; } // compute friction traction to stop sliding Array & traction = const_cast &>(this->friction_traction.getArray()); auto it_fric_trac = traction.begin(dim); for (UInt n = 0; n < nb_contact_nodes; ++n) { Vector fric_trac = it_fric_trac[n]; // node pair is NOT in contact if (!is_in_contact(n)) { fric_trac.zero(); // set to zero } // node pair is in contact else { // compute friction traction for (UInt d = 0; d < dim; ++d) fric_trac(d) = impedance(n) * gap_dot(n, d) / 2.; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::applyFrictionTraction() { AKANTU_DEBUG_IN(); SolidMechanicsModel & model = this->contact.getModel(); Array & residual = model.getInternalForce(); UInt dim = model.getSpatialDimension(); const SynchronizedArray & slaves = this->contact.getSlaves(); const SynchronizedArray & lumped_boundary_slaves = this->contact.getLumpedBoundarySlaves(); UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt slave = slaves(n); for (UInt d = 0; d < dim; ++d) { residual(slave, d) -= lumped_boundary_slaves(n) * this->friction_traction(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::registerSynchronizedArray(SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->frictional_strength.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->is_sticking.dumpRestartFile(file_name); this->frictional_strength.dumpRestartFile(file_name); this->friction_traction.dumpRestartFile(file_name); this->slip.dumpRestartFile(file_name); this->cumulative_slip.dumpRestartFile(file_name); this->slip_velocity.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->is_sticking.readRestartFile(file_name); this->frictional_strength.readRestartFile(file_name); this->friction_traction.readRestartFile(file_name); this->cumulative_slip.readRestartFile(file_name); this->slip_velocity.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::setParam(const std::string & name, UInt node, Real value) { AKANTU_DEBUG_IN(); SynchronizedArray & array = this->get(name).get>(); Int index = this->contact.getNodeIndex(node); if (index < 0) { AKANTU_DEBUG_WARNING("Node " << node << " is not a contact node. " << "Therefore, cannot set interface parameter!!"); } else { array(index) = value; // put value } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt NTNBaseFriction::getNbStickingNodes() const { AKANTU_DEBUG_IN(); UInt nb_stick = 0; UInt nb_nodes = this->contact.getNbContactNodes(); const SynchronizedArray & nodes = this->contact.getSlaves(); const SynchronizedArray & is_in_contact = this->contact.getIsInContact(); const Mesh & mesh = this->contact.getModel().getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(nodes(n)); bool is_pbc_slave_node = mesh.isPeriodicSlave(nodes(n)); if (is_local_node && !is_pbc_slave_node && is_in_contact(n) && this->is_sticking(n)) { nb_stick++; } } mesh.getCommunicator().allReduce(nb_stick, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_stick; } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNBaseFriction [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = // &(this->contact.getSlaves()); if (field_id == "is_sticking") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>( this->is_sticking.getArray())); } else if (field_id == "frictional_strength") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>( this->frictional_strength.getArray())); } else if (field_id == "friction_traction") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>( this->friction_traction.getArray())); } else if (field_id == "slip") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>(this->slip.getArray())); } else if (field_id == "cumulative_slip") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>( this->cumulative_slip.getArray())); } else if (field_id == "slip_velocity") { this->internalAddDumpFieldToDumper( dumper_name, field_id, std::make_unique>( this->slip_velocity.getArray())); } else { this->contact.addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh index 538d02f65..7b6ec0b96 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh @@ -1,177 +1,176 @@ /** * @file ntn_base_friction.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief base class for ntn and ntrf friction * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_BASE_FRICTION_HH_ #define AST_NTN_BASE_FRICTION_HH_ /* -------------------------------------------------------------------------- */ // akantu #include "parsable.hh" // simtools #include "ntn_base_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto( const ParserParameter & in_param) { Parameter::setAuto(in_param); Real r = in_param; param.setAndChangeDefault(r); } /* -------------------------------------------------------------------------- */ template <> template <> inline void ParameterTyped>::setTyped( const Real & value) { param.setAndChangeDefault(value); } /* -------------------------------------------------------------------------- */ -class NTNBaseFriction : protected Memory, public Parsable, public Dumpable { +class NTNBaseFriction : public Parsable, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNBaseFriction(NTNBaseContact & contact, const ID & id = "friction", - const MemoryID & memory_id = 0); + NTNBaseFriction(NTNBaseContact & contact, const ID & id = "friction"); virtual ~NTNBaseFriction() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// compute friction traction virtual void computeFrictionTraction(); /// compute stick traction (friction traction needed to stick the nodes) virtual void computeStickTraction(); /// apply the friction force virtual void applyFrictionTraction(); /// compute slip virtual void updateSlip(); /// register Syncronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// set to steady state virtual void setToSteadyState() { AKANTU_TO_IMPLEMENT(); }; /// get the number of sticking nodes (in parallel) /// a node that is not in contact does not count as sticking virtual UInt getNbStickingNodes() const; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength() = 0; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Contact, contact, const NTNBaseContact &) AKANTU_GET_MACRO(IsSticking, is_sticking, const SynchronizedArray &) AKANTU_GET_MACRO(FrictionalStrength, frictional_strength, const SynchronizedArray &) AKANTU_GET_MACRO(FrictionTraction, friction_traction, const SynchronizedArray &) AKANTU_GET_MACRO(Slip, slip, const SynchronizedArray &) AKANTU_GET_MACRO(CumulativeSlip, cumulative_slip, const SynchronizedArray &) AKANTU_GET_MACRO(SlipVelocity, slip_velocity, const SynchronizedArray &) /// set parameter of a given node /// (if you need to set to all: used the setMixed function of the Parsable). virtual void setParam(const std::string & name, UInt node, Real value); // replaced by the setMixed of the Parsable // virtual void setParam(const std::string & param, Real value) { // AKANTU_ERROR("Friction does not know the following parameter: " << // param); // }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: NTNBaseContact & contact; // if node is sticking SynchronizedArray is_sticking; // frictional strength SynchronizedArray frictional_strength; // friction force SynchronizedArray friction_traction; // slip SynchronizedArray slip; SynchronizedArray cumulative_slip; // slip velocity (tangential vector) SynchronizedArray slip_velocity; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_base_friction_inline_impl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNBaseFriction & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_NTN_BASE_FRICTION_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc index 17a44cac8..e5ad24f75 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc @@ -1,555 +1,554 @@ /** * @file ntn_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn_contact * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_contact.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -NTNContact::NTNContact(SolidMechanicsModel & model, const ID & id, - const MemoryID & memory_id) - : NTNBaseContact(model, id, memory_id), +NTNContact::NTNContact(SolidMechanicsModel & model, const ID & id) + : NTNBaseContact(model, id), masters(0, 1, 0, id + ":masters", std::numeric_limits::quiet_NaN(), "masters"), lumped_boundary_masters(0, 1, 0, id + ":lumped_boundary_masters", std::numeric_limits::quiet_NaN(), "lumped_boundary_masters"), - master_elements("master_elements", id, memory_id) { + master_elements("master_elements", id) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = this->model.getSpatialDimension(); this->master_elements.initialize(mesh, _nb_component = 1, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::pairInterfaceNodes(const ElementGroup & slave_boundary, const ElementGroup & master_boundary, UInt surface_normal_dir, const Mesh & mesh, Array & pairs) { AKANTU_DEBUG_IN(); pairs.resize(0); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(surface_normal_dir < dim, "Mesh is of " << dim << " dimensions" << " and cannot have direction " << surface_normal_dir << " for surface normal"); // offset for projection computation Vector offset(dim - 1); for (UInt i = 0, j = 0; i < dim; ++i) { if (surface_normal_dir != i) { offset(j) = i; ++j; } } // find projected node coordinates const Array & coordinates = mesh.getNodes(); // find slave nodes Array proj_slave_coord(slave_boundary.getNbNodes(), dim - 1, 0.); Array slave_nodes(slave_boundary.getNbNodes()); UInt n(0); for (auto && slave_node : slave_boundary.getNodeGroup().getNodes()) { for (UInt d = 0; d < dim - 1; ++d) { proj_slave_coord(n, d) = coordinates(slave_node, offset[d]); slave_nodes(n) = slave_node; } ++n; } // find master nodes Array proj_master_coord(master_boundary.getNbNodes(), dim - 1, 0.); Array master_nodes(master_boundary.getNbNodes()); n = 0; for (auto && master_node : master_boundary.getNodeGroup().getNodes()) { for (UInt d = 0; d < dim - 1; ++d) { proj_master_coord(n, d) = coordinates(master_node, offset[d]); master_nodes(n) = master_node; } ++n; } // find minimum distance between slave nodes to define tolerance Real min_dist = std::numeric_limits::max(); for (UInt i = 0; i < proj_slave_coord.size(); ++i) { for (UInt j = i + 1; j < proj_slave_coord.size(); ++j) { Real dist = 0.; for (UInt d = 0; d < dim - 1; ++d) { dist += (proj_slave_coord(i, d) - proj_slave_coord(j, d)) * (proj_slave_coord(i, d) - proj_slave_coord(j, d)); } if (dist < min_dist) { min_dist = dist; } } } min_dist = std::sqrt(min_dist); Real local_tol = 0.1 * min_dist; // find master slave node pairs for (UInt i = 0; i < proj_slave_coord.size(); ++i) { for (UInt j = 0; j < proj_master_coord.size(); ++j) { Real dist = 0.; for (UInt d = 0; d < dim - 1; ++d) { dist += (proj_slave_coord(i, d) - proj_master_coord(j, d)) * (proj_slave_coord(i, d) - proj_master_coord(j, d)); } dist = std::sqrt(dist); if (dist < local_tol) { // it is a pair Vector pair(2); pair[0] = slave_nodes(i); pair[1] = master_nodes(j); pairs.push_back(pair); continue; // found master do not need to search further for this slave } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addSurfacePair(const ID & slave, const ID & master, UInt surface_normal_dir) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model.getMesh(); const ElementGroup & slave_boundary = mesh.getElementGroup(slave); const ElementGroup & master_boundary = mesh.getElementGroup(master); this->contact_surfaces.insert(&slave_boundary); this->contact_surfaces.insert(&master_boundary); Array pairs(0, 2); NTNContact::pairInterfaceNodes(slave_boundary, master_boundary, surface_normal_dir, this->model.getMesh(), pairs); // eliminate pairs which contain a pbc slave node Array pairs_no_PBC_slaves(0, 2); Array::const_vector_iterator it = pairs.begin(2); Array::const_vector_iterator end = pairs.end(2); for (; it != end; ++it) { const Vector & pair = *it; if (not mesh.isPeriodicSlave(pair(0)) and not mesh.isPeriodicSlave(pair(1))) { pairs_no_PBC_slaves.push_back(pair); } } this->addNodePairs(pairs_no_PBC_slaves); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addNodePairs(const Array & pairs) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt nb_pairs = pairs.size(); for (UInt n = 0; n < nb_pairs; ++n) { this->addSplitNode(pairs(n, 0), pairs(n, 1)); } // synchronize with depending nodes findBoundaryElements(this->slaves.getArray(), this->slave_elements); findBoundaryElements(this->masters.getArray(), this->master_elements); updateInternalData(); syncArrays(_added); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::getNodePairs(Array & pairs) const { AKANTU_DEBUG_IN(); pairs.resize(0); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt nb_pairs = this->getNbContactNodes(); for (UInt n = 0; n < nb_pairs; ++n) { Vector pair{this->slaves(n), this->masters(n)}; pairs.push_back(pair); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addSplitNode(UInt slave, UInt master) { AKANTU_DEBUG_IN(); NTNBaseContact::addSplitNode(slave); this->masters.push_back(master); this->lumped_boundary_masters.push_back( std::numeric_limits::quiet_NaN()); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* This function only works for surface elements with one quad point. For surface elements with more quad points, it computes still, but the result might not be what you are looking for. */ void NTNContact::updateNormals() { AKANTU_DEBUG_IN(); // set normals to zero this->normals.zero(); // contact information UInt dim = this->model.getSpatialDimension(); UInt nb_contact_nodes = this->getNbContactNodes(); this->synch_registry->synchronize( SynchronizationTag::_cf_nodal); // synchronize current pos const Array & cur_pos = this->model.getCurrentPosition(); FEEngine & boundary_fem = this->model.getFEEngineBoundary(); const Mesh & mesh = this->model.getMesh(); for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(dim - 1, ghost_type)) { // compute the normals Array quad_normals(0, dim); boundary_fem.computeNormalsOnIntegrationPoints(cur_pos, quad_normals, type, ghost_type); UInt nb_quad_points = boundary_fem.getNbIntegrationPoints(type, ghost_type); // new version: compute normals only based on master elements (and not all // boundary elements) // ------------------------------------------------------------------------------------- UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, ghost_type); // loop over contact nodes for (auto & element : (this->master_elements)(type, ghost_type)) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { UInt node = connectivity(element, q); UInt node_index = this->masters.find(node); AKANTU_DEBUG_ASSERT(node_index != UInt(-1), "Could not find node " << node << " in the array!"); for (UInt q = 0; q < nb_quad_points; ++q) { // add quad normal to master normal for (UInt d = 0; d < dim; ++d) { this->normals(node_index, d) += quad_normals(element * nb_quad_points + q, d); } } } } } } Real * master_normals = this->normals.storage(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (dim == 2) Math::normalize2(&(master_normals[n * dim])); else if (dim == 3) Math::normalize3(&(master_normals[n * dim])); } // // normalize normals // auto nit = this->normals.begin(); // auto nend = this->normals.end(); // for (; nit != nend; ++nit) { // nit->normalize(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); NTNBaseContact::dumpRestart(file_name); this->masters.dumpRestartFile(file_name); this->lumped_boundary_masters.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); NTNBaseContact::readRestart(file_name); this->masters.readRestartFile(file_name); this->lumped_boundary_masters.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::updateImpedance() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); Real delta_t = this->model.getTimeStep(); AKANTU_DEBUG_ASSERT(delta_t != NAN, "Time step is NAN. Have you set it already?"); const Array & mass = this->model.getMass(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt master = this->masters(n); UInt slave = this->slaves(n); Real imp = (this->lumped_boundary_masters(n) / mass(master)) + (this->lumped_boundary_slaves(n) / mass(slave)); imp = 2 / delta_t / imp; this->impedance(n) = imp; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::updateLumpedBoundary() { AKANTU_DEBUG_IN(); internalUpdateLumpedBoundary(this->slaves.getArray(), this->slave_elements, this->lumped_boundary_slaves); internalUpdateLumpedBoundary(this->masters.getArray(), this->master_elements, this->lumped_boundary_masters); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::applyContactPressure() { AKANTU_DEBUG_IN(); UInt nb_ntn_pairs = getNbContactNodes(); UInt dim = this->model.getSpatialDimension(); Array & residual = this->model.getInternalForce(); for (UInt n = 0; n < nb_ntn_pairs; ++n) { UInt master = this->masters(n); UInt slave = this->slaves(n); for (UInt d = 0; d < dim; ++d) { residual(master, d) += this->lumped_boundary_masters(n) * this->contact_pressure(n, d); residual(slave, d) -= this->lumped_boundary_slaves(n) * this->contact_pressure(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::computeRelativeTangentialField( const Array & field, Array & rel_tang_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_tang_field.resize(0); UInt dim = this->model.getSpatialDimension(); auto it_field = field.begin(dim); auto it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); Vector np_rfv(dim); UInt nb_contact_nodes = this->slaves.size(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt slave = this->slaves(n); UInt master = this->masters(n); // relative field vector (slave - master) rfv = Vector(it_field[slave]); rfv -= Vector(it_field[master]); // normal projection of relative field const Vector normal_v = it_normal[n]; np_rfv = normal_v; np_rfv *= rfv.dot(normal_v); // subract normal projection from relative field to get the tangential // projection rfv -= np_rfv; rel_tang_field.push_back(rfv); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::computeRelativeNormalField( const Array & field, Array & rel_normal_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_normal_field.resize(0); UInt dim = this->model.getSpatialDimension(); // Real * field_p = field.storage(); // Real * normals_p = this->normals.storage(); Array::const_iterator> it_field = field.begin(dim); Array::const_iterator> it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); UInt nb_contact_nodes = this->getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt slave = this->slaves(n); UInt master = this->masters(n); // relative field vector (slave - master) rfv = Vector(it_field[slave]); rfv -= Vector(it_field[master]); // length of normal projection of relative field const Vector normal_v = it_normal[n]; rel_normal_field.push_back(rfv.dot(normal_v)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int NTNContact::getNodeIndex(UInt node) const { AKANTU_DEBUG_IN(); Int slave_i = NTNBaseContact::getNodeIndex(node); Int master_i = this->masters.find(node); AKANTU_DEBUG_OUT(); return std::max(slave_i, master_i); } /* -------------------------------------------------------------------------- */ void NTNContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNContact [" << std::endl; NTNBaseContact::printself(stream, indent); stream << space << " + masters : " << std::endl; this->masters.printself(stream, indent + 2); stream << space << " + lumped_boundary_mastres : " << std::endl; this->lumped_boundary_masters.printself(stream, indent + 2); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::syncArrays(SyncChoice sync_choice) { AKANTU_DEBUG_IN(); NTNBaseContact::syncArrays(sync_choice); this->masters.syncElements(sync_choice); this->lumped_boundary_masters.syncElements(sync_choice); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); /* #ifdef AKANTU_USE_IOHELPER const Array & nodal_filter = this->slaves.getArray(); #define ADD_FIELD(field_id, field, type) \ internalAddDumpFieldToDumper(dumper_name, \ field_id, \ new DumperIOHelper::NodalField< type, true, \ Array, \ Array >(field, 0, 0, &nodal_filter)) */ if (field_id == "lumped_boundary_master") { internalAddDumpFieldToDumper(dumper_name, field_id, std::make_unique>( this->lumped_boundary_masters.getArray())); } else { NTNBaseContact::addDumpFieldToDumper(dumper_name, field_id); } /* #undef ADD_FIELD #endif */ AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh index eb445dc20..dcd3beb69 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh @@ -1,166 +1,165 @@ /** * @file ntn_contact.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief contact for node to node discretization * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_CONTACT_HH_ #define AST_NTN_CONTACT_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNContact : public NTNBaseContact { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNContact(SolidMechanicsModel & model, const ID & id = "contact", - const MemoryID & memory_id = 0); - virtual ~NTNContact(){}; + NTNContact(SolidMechanicsModel & model, const ID & id = "contact"); + ~NTNContact() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// add surface pair and pair nodes according to the surface normal void addSurfacePair(const ID & slave, const ID & master, UInt surface_normal_dir); /// fills the pairs vector with interface node pairs (*,0)=slaves, /// (*,1)=masters static void pairInterfaceNodes(const ElementGroup & slave_boundary, const ElementGroup & master_boundary, UInt surface_normal_dir, const Mesh & mesh, Array & pairs); // add node pairs from a list with pairs(*,0)=slaves and pairs(*,1)=masters void addNodePairs(const Array & pairs); /// add node pair void addSplitNode(UInt slave, UInt master) override; /// update (compute the normals on the master nodes) void updateNormals() override; /// update the lumped boundary B matrix void updateLumpedBoundary() override; /// update the impedance matrix void updateImpedance() override; /// impose the normal contact force void applyContactPressure() override; /// dump restart file void dumpRestart(const std::string & file_name) const override; /// read restart file void readRestart(const std::string & file_name) override; /// compute the normal gap void computeNormalGap(Array & gap) const override { this->computeRelativeNormalField(this->model.getCurrentPosition(), gap); }; /// compute relative normal field (only value that has to be multiplied with /// the normal) /// relative to master nodes void computeRelativeNormalField(const Array & field, Array & rel_normal_field) const override; /// compute relative tangential field (complet array) /// relative to master nodes void computeRelativeTangentialField(const Array & field, Array & rel_tang_field) const override; /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; protected: /// synchronize arrays void syncArrays(SyncChoice sync_choice) override; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) override; // virtual void addDumpFieldVector(const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Masters, masters, const SynchronizedArray &) AKANTU_GET_MACRO(LumpedBoundaryMasters, lumped_boundary_masters, const SynchronizedArray &) /// get interface node pairs (*,0) are slaves, (*,1) are masters void getNodePairs(Array & pairs) const; /// get index of node in either slaves or masters array /// if node is in neither of them, return -1 Int getNodeIndex(UInt node) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of master nodes SynchronizedArray masters; /// lumped boundary of master nodes SynchronizedArray lumped_boundary_masters; // element list for dump and lumped_boundary ElementTypeMapArray master_elements; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_contact_inline_impl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNContact & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_NTN_CONTACT_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh index dec3a550a..a5f43b29a 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh @@ -1,99 +1,98 @@ /** * @file ntn_friction.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief implementation of friction for node to node contact * * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_NTN_FRICTION_HH_ #define AST_NTN_FRICTION_HH_ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_friction.hh" #include "ntn_friclaw_coulomb.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template