diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 65747616a..ab607bb50 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,350 +1,345 @@ # yaml-language-server: $schema=gitlab-ci # yaml-language-server: $format.enable=false stages: - configure - build_libs - build_tests - test - code_quality - deploy include: local: '.gitlab-ci.d/*.yaml' #------------------------------------------------------------------------------- # Rebuilding the docker images if needed #------------------------------------------------------------------------------- docker build:debian-bullseye: variables: IMAGE_NAME: debian:bullseye extends: .docker_build docker build:ubuntu-lts: variables: IMAGE_NAME: ubuntu:lts extends: .docker_build docker build:manylinux: variables: IMAGE_NAME: manylinux:2010_x86_64 extends: .docker_build # ------------------------------------------------------------------------------ # Debian bullseye compiled with GCC # ------------------------------------------------------------------------------ configure:debian_bullseye_gcc: extends: - .debian_bullseye_gcc - .build_coverage - .configure build:debian_bullseye_gcc: extends: - .debian_bullseye_gcc - .build_coverage - .build_all needs: - job: configure:debian_bullseye_gcc test:debian_bullseye_gcc: extends: - .debian_bullseye_gcc - .build_coverage - .tests coverage: '/^lines: (\d+\.\d+\%)/' needs: - job: build:debian_bullseye_gcc # ------------------------------------------------------------------------------ # Debian bullseye compiled with Clang # ------------------------------------------------------------------------------ configure:debian_bullseye_clang: extends: - .debian_bullseye_clang - .build_coverage - .configure build:debian_bullseye_clang: extends: - .debian_bullseye_clang - .build_coverage - .build_all needs: - job: configure:debian_bullseye_clang test:debian_bullseye_clang: extends: - .debian_bullseye_clang - .build_coverage - .tests coverage: '/^lines: (\d+\.\d+\%)/' needs: - job: build:debian_bullseye_clang # ------------------------------------------------------------------------------ # Ubuntu LTS compiled with GCC # ------------------------------------------------------------------------------ configure:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build_release - .configure build:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build_release - .build_all needs: - job: configure:ubuntu_lts_gcc test:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build_release - .tests needs: - job: build:ubuntu_lts_gcc # ------------------------------------------------------------------------------ # Debian bullseye compiled with GCC tested with valgrind # ------------------------------------------------------------------------------ configure:ubuntu_lts_gcc_valgrind: extends: - .ubuntu_lts_gcc - .build_valgrind - .configure build:ubuntu_lts_gcc_valgrind: extends: - .ubuntu_lts_gcc - .build_valgrind - .build_all needs: - job: configure:ubuntu_lts_gcc_valgrind test:ubuntu_lts_gcc_valgrind: extends: - .ubuntu_lts_gcc - .build_valgrind - .tests needs: - job: build:ubuntu_lts_gcc_valgrind # ------------------------------------------------------------------------------ # Manylinux to build python packages # ------------------------------------------------------------------------------ configure:python_package: stage: configure extends: - .manylinux_2010_x64_gcc - .build_release script: # create the build folder - cmake -E make_directory build - cd build # Variables for cmake - export CMAKE_PREFIX_PATH=/softs/view - export BOOST_ROOT=/softs/view # Configure in sequential and without tests or examples - cmake -DAKANTU_COHESIVE_ELEMENT:BOOL=TRUE -DAKANTU_IMPLICIT:BOOL=TRUE -DAKANTU_PARALLEL:BOOL=FALSE -DAKANTU_STRUCTURAL_MECHANICS:BOOL=TRUE -DAKANTU_HEAT_TRANSFER:BOOL=TRUE -DAKANTU_DAMAGE_NON_LOCAL:BOOL=TRUE -DAKANTU_PHASE_FIELD:BOOL=TRUE -DAKANTU_PYTHON_INTERFACE:BOOL=FALSE -DAKANTU_CONTACT_MECHANICS:BOOL=TRUE -DAKANTU_EXAMPLES:BOOL=FALSE -DAKANTU_TESTS:BOOL=FALSE -DMUMPS_DETECT_DEBUG:BOOL=TRUE -DCMAKE_INSTALL_PREFIX:PATH=${CI_PROJECT_DIR}/install -DCMAKE_BUILD_TYPE:STRING=${BUILD_TYPE} .. artifacts: when: on_success paths: - build/ expire_in: 10h build_akantu:python_package: extends: - .build_libs - .build_release - .manylinux_2010_x64_gcc - script: - stage: build_libs script: - cmake --build build --target akantu -j1 - cmake --install build artifacts: when: on_success paths: - install/ expire_in: 10h needs: - job: configure:python_package build_pip:python_package: stage: build_tests extends: - .build_release - .manylinux_2010_x64_gcc script: - export CI_AKANTU_INSTALL_PREFIX=${CI_PROJECT_DIR}/install - export CMAKE_PREFIX_PATH=/softs/view:${CI_AKANTU_INSTALL_PREFIX} - test/ci/make-wheels.sh needs: - job: build_akantu:python_package artifacts: when: on_success paths: - wheelhouse expire_in: 10h test:python_package: stage: test image: python:3.8 needs: - job: build_pip:python_package script: - pip install numpy scipy - pip install akantu --no-index --find-links=${PWD}/wheelhouse - python -c "import akantu" - cd examples/python/dynamics/ - apt update && apt install -y gmsh - gmsh -2 bar.geo - python ./dynamics.py package:python_gitlab: stage: deploy image: python:latest script: - pip install twine - TWINE_PASSWORD=${CI_JOB_TOKEN} TWINE_USERNAME=gitlab-ci-token python3 -m twine upload --repository-url https://gitlab.com/api/v4/projects/${CI_PROJECT_ID}/packages/pypi wheelhouse/* needs: - job: build_pip:python_package - job: test:python_package only: - master package:python_pypi: stage: deploy image: python:latest script: - pip install twine - TWINE_PASSWORD=${PYPI_TOKEN} TWINE_USERNAME=__token__ python3 -m twine upload --verbose wheelhouse/* needs: - job: build_pip:python_package - job: test:python_package only: - tags # ------------------------------------------------------------------------------ # Code Quality # ------------------------------------------------------------------------------ cq:code_quality: extends: - .code_quality_gitlab_template needs: - job: build:debian_bullseye_clang artifacts: paths: - gl-code-quality-report.json cq:clang_tidy: extends: - .clang_tools script: - test/ci/scripts/cq -x third-party -x extra-packages -x pybind11 -x test ${FILE_LIST_ARG} clang-tidy -p ${CI_PROJECT_DIR}/build > gl-clang-tidy-report.json needs: - job: build:debian_bullseye_clang artifacts: paths: - gl-clang-tidy-report.json cq:clang_format: extends: - .clang_tools script: - test/ci/scripts/cq -x third-party -x extra-packages clang-format -p ${CI_PROJECT_DIR}/build > gl-clang-format-report.json needs: - job: build:debian_bullseye_clang artifacts: paths: - gl-clang-format-report.json cq:compilation_warnings: stage: code_quality image: python:latest script: - pip install warning-parser termcolor Click - ls build-*-err.log - test/ci/scripts/cq -x third-party -x extra-packages warnings build-*-err.log > gl-warnings-report.json needs: - job: build:debian_bullseye_clang - job: build:debian_bullseye_gcc - job: build:ubuntu_lts_gcc artifacts: paths: - gl-warnings-report.json cq:merge_code_quality: stage: deploy extends: - .debian_bullseye_clang script: - jq -Ms '[.[][]]' gl-*-report.json | tee gl-codequality.json | jq -C needs: - job: cq:code_quality - job: cq:clang_tidy - job: cq:clang_format - job: cq:compilation_warnings - artifacts: - paths: - - gl-codequality.json artifacts: reports: codequality: [gl-codequality.json] # ------------------------------------------------------------------------------ # Deploy pages # ------------------------------------------------------------------------------ pages: stage: deploy extends: - .debian_bullseye_gcc script: - cd build - cmake -DAKANTU_DOCUMENTATION=ON .. - cmake --build . -t sphinx-doc - mv doc/dev-doc/html ../public needs: - job: build:debian_bullseye_gcc artifacts: paths: - public only: - master diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 4c76a6e89..5cd6a4456 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,144 +1,145 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Fri Dec 12 2014 # @date last modification: Fri May 07 2021 # # @brief CMake file for the python wrapping of akantu # # # @section LICENSE # # Copyright (©) 2015-2021 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 . # #=============================================================================== set(PYBIND11_PYTHON_VERSION ${AKANTU_PREFERRED_PYTHON_VERSION} CACHE INTERNAL "") if(NOT SKBUILD) package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) package_get_all_external_informations( PRIVATE_INCLUDE AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR INTERFACE_INCLUDE AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR LIBRARIES AKANTU_EXTERNAL_LIBRARIES ) endif() set(PYAKANTU_SRCS py_aka_common.cc py_aka_error.cc py_akantu.cc py_boundary_conditions.cc py_dof_manager.cc py_fe_engine.cc py_group_manager.cc + py_integration_scheme.cc py_mesh.cc py_model.cc py_parser.cc py_solver.cc py_dumpable.cc ) package_is_activated(solid_mechanics _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model.cc py_material.cc py_material_selector.cc ) endif() package_is_activated(cohesive_element _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model_cohesive.cc py_fragment_manager.cc ) endif() package_is_activated(heat_transfer _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_heat_transfer_model.cc ) endif() package_is_activated(contact_mechanics _is_activated) if(_is_activated) list(APPEND PYAKANTU_SRCS py_contact_mechanics_model.cc py_model_couplers.cc ) endif() package_is_activated(phase_field _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_phase_field_model.cc ) endif() package_is_activated(structural_mechanics _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_structural_mechanics_model.cc ) endif() pybind11_add_module(py11_akantu ${PYAKANTU_SRCS}) # to avoid compilation warnings from pybind11 target_include_directories(py11_akantu SYSTEM BEFORE PRIVATE ${PYBIND11_INCLUDE_DIR} PRIVATE ${pybind11_INCLUDE_DIR} PRIVATE ${Python_INCLUDE_DIRS}) target_link_libraries(py11_akantu PUBLIC akantu) set_target_properties(py11_akantu PROPERTIES DEBUG_POSTFIX "" LIBRARY_OUTPUT_DIRECTORY akantu) file(COPY akantu DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) if(NOT Python_MAJOR) set(Python_VERSION_MAJOR ${PYTHON_VERSION_MAJOR}) set(Python_VERSION_MINOR ${PYTHON_VERSION_MINOR}) endif() if(NOT SKBUILD) set(_python_install_dir ${CMAKE_INSTALL_LIBDIR}/python${Python_VERSION_MAJOR}.${Python_VERSION_MINOR}/site-packages/akantu) else() set(_python_install_dir python/akantu) endif() install(TARGETS py11_akantu LIBRARY DESTINATION ${_python_install_dir}) if(NOT SKBUILD) install(DIRECTORY akantu DESTINATION ${_python_install_dir} FILES_MATCHING PATTERN "*.py") endif() diff --git a/python/py_integration_scheme.cc b/python/py_integration_scheme.cc index 62d31731a..9485d3bb6 100644 --- a/python/py_integration_scheme.cc +++ b/python/py_integration_scheme.cc @@ -1,56 +1,63 @@ /** * Copyright (©) 2022 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 "py_integration_scheme.hh" /* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; namespace akantu { /* -------------------------------------------------------------------------- */ void register_integration_schemes(py::module & mod) { - py::class_(mod, "IntegrationScheme") - .def(py::init()); + auto integration_scheme_class = + py::class_(mod, "IntegrationScheme"); - py::enum_() - .def("_not_defined", _not_defined) - .def("_displacement", _displacement) - .def("_temperature", _temperature) - .def("_damage", _damage) - .def("_velocity", _velocity) - .def("_temperature_rate", _temperature_rate) - .def("_acceleration", _acceleration); + py::enum_(integration_scheme_class, + "SolutionType") + .value("_not_defined", IntegrationScheme::SolutionType::_not_defined) + .value("_displacement", IntegrationScheme::SolutionType::_displacement) + .value("_temperature", IntegrationScheme::SolutionType::_temperature) + .value("_damage", IntegrationScheme::SolutionType::_damage) + .value("_velocity", IntegrationScheme::SolutionType::_velocity) + .value("_temperature_rate", + IntegrationScheme::SolutionType::_temperature_rate) + .value("_acceleration", IntegrationScheme::SolutionType::_acceleration) + .export_values(); py::class_( - mod, "IntegrationScheme2ndOrder") - .def(py::init()); + mod, "IntegrationScheme2ndOrder"); py::class_(mod, "NewmarkBeta") - .def(py::init(), + .def(py::init(), py::arg("dof_manager"), py::arg("id"), py::arg("alpha"), py::arg("beta")); py::class_(mod, "CentralDifference"); py::class_(mod, "TrapezoidalRule2"); py::class_(mod, "FoxGoodwin"); py::class_(mod, "LinearAceleration"); } } // namespace akantu diff --git a/python/py_mesh.cc b/python/py_mesh.cc index c371d9c6a..a73250d8f 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,237 +1,238 @@ /** * @file py_mesh.cc * * @author Guillaume Anciaux * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Mon Mar 15 2021 * * @brief pybind11 interface to Mesh * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { /* ------------------------------------------------------------------------ */ template decltype(auto) register_element_type_map_array(py::module & mod, const std::string & name) { return py::class_, std::shared_ptr>>( mod, ("ElementTypeMapArray" + name).c_str()) .def(py::init(), py::arg("id") = "by_element_type_array", py::arg("parent_id") = "no_parent") .def( "__call__", [](ElementTypeMapArray & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference, py::keep_alive<0, 1>()) .def( "elementTypes", [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, ElementKind _kind) -> std::vector { auto types = self.elementTypes(_dim, _ghost_type, _kind); std::vector _types; for (auto && t : types) { _types.push_back(t); } return _types; }, py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular) .def( "initialize", [](ElementTypeMapArray & self, const Mesh & mesh, GhostType ghost_type = _casper, UInt nb_component = 1, UInt spatial_dimension = UInt(-2), ElementKind element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false, T default_value = T(), bool do_not_default = false) { self.initialize( mesh, _ghost_type = ghost_type, _nb_component = nb_component, _spatial_dimension = (spatial_dimension == UInt(-2) ? mesh.getSpatialDimension() : spatial_dimension), _element_kind = element_kind, _with_nb_element = with_nb_element, _with_nb_nodes_per_element = with_nb_nodes_per_element, _default_value = default_value, _do_not_default = do_not_default); }, py::arg("mesh"), py::arg("ghost_type") = _casper, py::arg("nb_component") = 1, py::arg("spatial_dimension") = UInt(-2), py::arg("element_kind") = _ek_not_defined, py::arg("with_nb_element") = false, py::arg("with_nb_nodes_per_element") = false, py::arg("default_value") = T(), py::arg("do_not_default") = false); } } // namespace /* -------------------------------------------------------------------------- */ void register_mesh(py::module & mod) { py::class_(mod, "PeriodicSlaves") .def( "__iter__", [](Mesh::PeriodicSlaves & _this) { py::make_iterator(_this.begin(), _this.end()); }, py::keep_alive<0, 1>()); py::class_(mod, "MeshData") .def( "getElementalDataUInt", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference) .def( "getElementalDataReal", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference); py::class_(mod, "Mesh", py::multiple_inheritance()) .def(py::init(), py::arg("spatial_dimension"), py::arg("id") = "mesh") .def("read", &Mesh::read, py::arg("filename"), py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file") .def( "getNodes", [](Mesh & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("getNbNodes", &Mesh::getNbNodes) .def( "getConnectivity", [](Mesh & self, ElementType type) -> decltype(auto) { return self.getConnectivity(type); }, py::return_value_policy::reference) .def( "addConnectivityType", [](Mesh & self, ElementType type, GhostType ghost_type) -> void { self.addConnectivityType(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("distribute", [](Mesh & self) { self.distribute(); }) + .def("isDistributed", [](const Mesh& self) { return self.isDistributed(); }) .def("fillNodesToElements", &Mesh::fillNodesToElements, py::arg("dimension") = _all_dimensions) .def("getAssociatedElements", [](Mesh & self, const UInt & node, py::list list) { Array elements; self.getAssociatedElements(node, elements); for (auto && element : elements) { list.append(element); } }) .def("makePeriodic", [](Mesh & self, const SpatialDirection & direction) { self.makePeriodic(direction); }) .def( "getNbElement", [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type, ElementKind kind) { return self.getNbElement(spatial_dimension, ghost_type, kind); }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined) .def( "getNbElement", [](Mesh & self, ElementType type, GhostType ghost_type) { return self.getNbElement(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def_static( "getSpatialDimension", [](ElementType & type) { return Mesh::getSpatialDimension(type); }) .def( "getDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> decltype(auto) { return _this.getData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "hasDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> bool { return _this.hasData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); }) .def("getPeriodicMaster", &Mesh::getPeriodicMaster) .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves) .def("isPeriodicSlave", &Mesh::isPeriodicSlave) .def("isPeriodicMaster", &Mesh::isPeriodicMaster); /* ------------------------------------------------------------------------ */ py::class_(mod, "MeshUtils") .def_static("buildFacets", &MeshUtils::buildFacets); py::class_(mod, "MeshAccessor") .def(py::init(), py::arg("mesh")) .def( "resizeConnectivity", [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt) -> void { self.resizeConnectivity(new_size, type, gt); }, py::arg("new_size"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "resizeNodes", [](MeshAccessor & self, UInt new_size) -> void { self.resizeNodes(new_size); }, py::arg("new_size")) .def("makeReady", &MeshAccessor::makeReady); register_element_type_map_array(mod, "Real"); register_element_type_map_array(mod, "UInt"); // register_element_type_map_array(mod, "String"); } } // namespace akantu diff --git a/python/py_model.cc b/python/py_model.cc index 4328a40c0..4f0a9312c 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,145 +1,149 @@ /** * @file py_model.cc * * @author Guillaume Anciaux * @author Emil Gallyamov * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Sat Mar 13 2021 * * @brief pybind11 interface to Model and parent classes * * * @section LICENSE * * Copyright (©) 2018-2021 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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_model(py::module & mod) { py::class_(mod, "ModelSolver", py::multiple_inheritance()) .def( "getNonLinearSolver", [](ModelSolver & self, const ID & solver_id) -> NonLinearSolver & { return self.getNonLinearSolver(solver_id); }, py::arg("solver_id") = "", py::return_value_policy::reference) .def( "getTimeStepSolver", [](ModelSolver & self, const ID & solver_id) -> TimeStepSolver & { return self.getTimeStepSolver(solver_id); }, py::arg("solver_id") = "", py::return_value_policy::reference) .def( "solveStep", [](ModelSolver & self, const ID & solver_id) { self.solveStep(solver_id); }, py::arg("solver_id") = "") .def( "solveStep", [](ModelSolver & self, SolverCallback & callback, const ID & solver_id) { self.solveStep(callback, solver_id); }, py::arg("callback"), py::arg("solver_id") = ""); py::class_(mod, "Model", py::multiple_inheritance()) .def("setBaseName", &Model::setBaseName) .def("setDirectory", &Model::setDirectory) .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("addDumpFieldVector", &Model::addDumpFieldVector) .def("addDumpField", &Model::addDumpField) .def("setBaseNameToDumper", &Model::setBaseNameToDumper) .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper) .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper) .def("dump", [](Model & self) { self.dump(); }) .def( "dump", [](Model & self, UInt step) { self.dump(step); }, py::arg("step")) .def( "dump", [](Model & self, Real time, UInt step) { self.dump(time, step); }, py::arg("time"), py::arg("step")) .def( "dump", [](Model & self, const std::string & dumper) { self.dump(dumper); }, py::arg("dumper_name")) .def( "dump", [](Model & self, const std::string & dumper, UInt step) { self.dump(dumper, step); }, py::arg("dumper_name"), py::arg("step")) .def( "dump", [](Model & self, const std::string & dumper, Real time, UInt step) { self.dump(dumper, time, step); }, py::arg("dumper_name"), py::arg("time"), py::arg("step")) .def("initNewSolver", &Model::initNewSolver) .def( "getNewSolver", [](Model & self, const std::string id, const TimeStepSolverType & time, const NonLinearSolverType & type) { self.getNewSolver(id, time, type); }, py::return_value_policy::reference) - .def("setIntegrationScheme", - [](Model & self, const std::string id, const std::string primal, - const IntegrationSchemeType & scheme_type, - IntegrationScheme::SolutionType solution_type) { - self.setIntegrationScheme(id, primal, scheme_type, solution_type); - }) + .def( + "setIntegrationScheme", + [](Model & self, const std::string id, const std::string primal, + const IntegrationSchemeType & scheme_type, + IntegrationScheme::SolutionType solution_type) { + self.setIntegrationScheme(id, primal, scheme_type, solution_type); + }, + py::arg("id"), py::arg("primal"), py::arg("scheme_type"), + py::arg("solution_type") = + IntegrationScheme::SolutionType::_not_defined) // .def("setIntegrationScheme", // [](Model & self, const std::string id, const std::string primal, // std::unique_ptr & scheme, // IntegrationScheme::SolutionType solution_type) { // self.setIntegrationScheme(id, primal, scheme, solution_type); // }) .def("getDOFManager", &Model::getDOFManager, py::return_value_policy::reference) .def("assembleMatrix", &Model::assembleMatrix); } } // namespace akantu diff --git a/python/py_solver.cc b/python/py_solver.cc index 1f16c3acb..c1239854a 100644 --- a/python/py_solver.cc +++ b/python/py_solver.cc @@ -1,111 +1,119 @@ /** * @file py_solver.cc * * @author Nicolas Richart * * @date creation: Tue Sep 29 2020 * @date last modification: Sat Mar 06 2021 * * @brief pybind11 interface to Solver and SparseMatrix * * * @section LICENSE * * Copyright (©) 2018-2021 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 "py_solver.hh" #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_solvers(py::module & mod) { py::class_(mod, "SparseMatrix") .def("getMatrixType", &SparseMatrix::getMatrixType) .def("size", &SparseMatrix::size) .def("zero", &SparseMatrix::zero) .def("saveProfile", &SparseMatrix::saveProfile) .def("saveMatrix", &SparseMatrix::saveMatrix) .def( "add", [](SparseMatrix & self, UInt i, UInt j) { self.add(i, j); }, "Add entry in the profile") .def( "add", [](SparseMatrix & self, UInt i, UInt j, Real value) { self.add(i, j, value); }, "Add the value to the matrix") .def( "add", [](SparseMatrix & self, SparseMatrix & A, Real alpha) { self.add(A, alpha); }, "Add a matrix to the matrix", py::arg("A"), py::arg("alpha") = 1.) .def("isFinite", &SparseMatrix::isFinite) .def("getRelease", [](const SparseMatrix& self) -> UInt { return self.getRelease(); }) .def("__call__", [](const SparseMatrix & self, UInt i, UInt j) { return self(i, j); }) .def("getRelease", &SparseMatrix::getRelease); py::class_(mod, "SparseMatrixAIJ") .def("getIRN", &SparseMatrixAIJ::getIRN) .def("getJCN", &SparseMatrixAIJ::getJCN) .def("getA", &SparseMatrixAIJ::getA); - py::class_(mod, "SolverVector"); + py::class_(mod, "SolverVector") + .def("getValues", + [](SolverVector& self) -> decltype(auto) { + return static_cast& >(self); + }, + py::return_value_policy::reference_internal, + "Transform this into a vector, Is not copied.") + .def("isDistributed", [](const SolverVector& self) { return self.isDistributed(); }) + ; py::class_(mod, "TermToAssemble") .def(py::init()) .def(py::self += Real()) .def_property_readonly("i", &TermsToAssemble::TermToAssemble::i) .def_property_readonly("j", &TermsToAssemble::TermToAssemble::j); py::class_(mod, "TermsToAssemble") .def(py::init()) .def("getDOFIdM", &TermsToAssemble::getDOFIdM) .def("getDOFIdN", &TermsToAssemble::getDOFIdN) .def( "__call__", [](TermsToAssemble & self, UInt i, UInt j, Real val) { auto & term = self(i, j); term = val; return term; }, py::arg("i"), py::arg("j"), py::arg("val") = 0., py::return_value_policy::reference); } } // namespace akantu diff --git a/python/py_structural_mechanics_model.cc b/python/py_structural_mechanics_model.cc index 9d160bc79..79801516f 100644 --- a/python/py_structural_mechanics_model.cc +++ b/python/py_structural_mechanics_model.cc @@ -1,172 +1,177 @@ /** * @file py_structural_mechanics_model.cc * * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Wed Feb 03 2021 * @date last modification: Thu Apr 01 2021 * * @brief pybind11 interface to StructuralMechanicsModel * * * @section LICENSE * * Copyright (©) 2018-2021 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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](StructuralMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function_(func_name) \ def(#func_name, [](StructuralMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) #define def_plainmember(M) def_readwrite(#M, &StructuralMaterial::M) /* -------------------------------------------------------------------------- */ void register_structural_mechanics_model(pybind11::module & mod) { /* First we have to register the material class * The wrapper aims to mimic the behaviour of the real material. */ py::class_(mod, "StructuralMaterial") .def(py::init<>()) .def(py::init()) .def_plainmember(E) .def_plainmember(A) .def_plainmember(I) .def_plainmember(Iz) .def_plainmember(Iy) .def_plainmember(GJ) .def_plainmember(rho) .def_plainmember(t) .def_plainmember(nu); /* Now we create the structural model wrapper * Note that this is basically a port from the solid mechanic part. */ py::class_(mod, "StructuralMechanicsModel") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "structural_mechanics_model") .def( "initFull", [](StructuralMechanicsModel & self, const AnalysisMethod & analysis_method) -> void { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def("initFull", [](StructuralMechanicsModel & self) -> void { self.initFull(); }) .def_function_nocopy(getExternalForce) .def_function_nocopy(getDisplacement) .def_function_nocopy(getInternalForce) .def_function_nocopy(getVelocity) .def_function_nocopy(getAcceleration) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getMesh) .def("setTimeStep", &StructuralMechanicsModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def( "getElementMaterial", [](StructuralMechanicsModel & self, const ElementType & type, GhostType ghost_type) -> decltype(auto) { return self.getElementMaterial(type, ghost_type); }, "This function returns the map that maps elements to materials.", py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getMaterialByElement", [](StructuralMechanicsModel & self, Element element) -> decltype(auto) { return self.getMaterialByElement(element); }, "This function returns the `StructuralMaterial` instance that is " "associated with element `element`.", py::arg("element"), py::return_value_policy::reference) .def( "addMaterial", [](StructuralMechanicsModel & self, StructuralMaterial & mat, const ID & name) -> UInt { return self.addMaterial(mat, name); }, "This function adds the `StructuralMaterial` `mat` to `self`." " The function returns the ID of the new material.", py::arg("mat"), py::arg("name") = "") .def( "getMaterial", - [](StructuralMechanicsModel & self, UInt material_index) - -> decltype(auto) { return self.getMaterial(material_index); }, + [](const StructuralMechanicsModel & self, UInt material_index) + -> StructuralMaterial { return self.getMaterial(material_index); }, "This function returns the `i`th material of `self`." - " Note a reference is returned which allows to modify the material " - "inside `self`.", - py::arg("material_index"), py::return_value_policy::reference) + " The function returns a copy of the material.", + py::arg("material_index"), py::return_value_policy::copy) .def( "getMaterial", - [](StructuralMechanicsModel & self, const ID & name) - -> decltype(auto) { return self.getMaterial(name); }, - "This function returns the material with name `i` of `self`." - " Note a reference is returned which allows to modify the material " - "inside `self`.", - py::arg("material_index"), py::return_value_policy::reference) + [](const StructuralMechanicsModel & self, const ID & name) + -> StructuralMaterial { return self.getMaterial(name); }, + "This function returns the `i`th material of `self`." + " The function returns a copy of the material.", + py::arg("material_index"), py::return_value_policy::copy) .def( "getNbMaterials", [](StructuralMechanicsModel & self) { return self.getNbMaterials(); }, "Returns the number of different materials inside `self`.") + .def("getKineticEnergy", &StructuralMechanicsModel::getKineticEnergy, "Compute kinetic energy") .def("getPotentialEnergy", &StructuralMechanicsModel::getPotentialEnergy, "Compute potential energy") .def("getEnergy", &StructuralMechanicsModel::getEnergy, "Compute the specified energy") + .def( "getLumpedMass", - [](StructuralMechanicsModel & self) -> decltype(auto) { + [](const StructuralMechanicsModel & self) -> decltype(auto) { return self.getLumpedMass(); }, py::return_value_policy::reference_internal) + .def( + "getMass", + [](const StructuralMechanicsModel & self) -> decltype(auto) { + return self.getMass(); + }, + py::return_value_policy::reference_internal) .def("assembleLumpedMassMatrix", &StructuralMechanicsModel::assembleLumpedMassMatrix, "Assembles the lumped mass matrix") - .def("hasLumpedMass", &StructuralMechanicsModel::hasLumpedMass); } } // namespace akantu diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index 267a8944b..61e4d3df0 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,656 +1,658 @@ /** * @file mesh.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Mohit Pundir * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 09 2021 * * @brief class handling meshes * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_global_data_updater.hh" #include "mesh_io.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" #include "periodic_node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "dumper_internal_material_field.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, Communicator & communicator) : GroupManager(*this, id + ":group_manager"), MeshData("mesh_data", id), id(id), connectivities("connectivities", id), ghosts_counters("ghosts_counters", id), normals("normals", id), spatial_dimension(spatial_dimension), size(spatial_dimension, 0.), bbox(spatial_dimension), bbox_local(spatial_dimension), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id) : Mesh(spatial_dimension, id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); this->nodes_flags = std::make_shared>(0, 1, NodeFlag::_normal, id + ":nodes_flags"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id) : Mesh(spatial_dimension, id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, ElementType type, GhostType ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } class FacetGlobalConnectivityAccessor : public DataAccessor { public: FacetGlobalConnectivityAccessor(Mesh & mesh) : global_connectivity("global_connectivity", "facet_connectivity_synchronizer") { global_connectivity.initialize( mesh, _spatial_dimension = _all_dimensions, _with_nb_element = true, _with_nb_nodes_per_element = true, _element_kind = _ek_regular); mesh.getGlobalConnectivity(global_connectivity); } UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override { UInt size = 0; if (tag == SynchronizationTag::_smmc_facets_conn) { UInt nb_nodes = Mesh::getNbNodesPerElementList(elements); size += nb_nodes * sizeof(UInt); } return size; } void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_smmc_facets_conn) { for (const auto & element : elements) { const auto & conns = global_connectivity(element.type, element.ghost_type); for (auto n : arange(conns.getNbComponent())) { buffer << conns(element.element, n); } } } } void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_smmc_facets_conn) { for (const auto & element : elements) { auto & conns = global_connectivity(element.type, element.ghost_type); for (auto n : arange(conns.getNbComponent())) { buffer >> conns(element.element, n); } } } } AKANTU_GET_MACRO(GlobalConnectivity, (global_connectivity), decltype(auto)); protected: ElementTypeMapArray global_connectivity; }; /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (mesh_facets) { AKANTU_DEBUG_OUT(); return *mesh_facets; } mesh_facets = std::make_unique(spatial_dimension, this->nodes, getID() + ":" + id); mesh_facets->mesh_parent = this; mesh_facets->is_mesh_facets = true; mesh_facets->nodes_flags = this->nodes_flags; mesh_facets->nodes_global_ids = this->nodes_global_ids; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); FacetGlobalConnectivityAccessor data_accessor(*mesh_facets); /// communicate mesh_facets->element_synchronizer->synchronizeOnce( data_accessor, SynchronizationTag::_smmc_facets_conn); /// flip facets MeshUtils::flipFacets(*mesh_facets, data_accessor.getGlobalConnectivity(), _ghost); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } auto & mesh_phys_data = this->getData("physical_names"); auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes(spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) { tolerance *= norm_barycenter; } Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view( barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name auto && facet_element = Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}; phys_data(facet_element) = mesh_phys_data(element); }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { AKANTU_DEBUG_ASSERT(not is_distributed, "You cannot read a mesh that is already distributed"); MeshIO::read(filename, *this, mesh_io_type); auto types = this->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); auto it = types.begin(); auto last = types.end(); if (it == last) { AKANTU_DEBUG_WARNING( "The mesh contained in the file " << filename << " does not seem to be of the good dimension." << " No element of dimension " << spatial_dimension << " were read."); } this->makeReady(); } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO::write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::makeReady() { this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_flags->resize(nodes->size(), NodeFlag::_normal); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); bbox_local.reset(); for (auto & pos : make_view(*nodes, spatial_dimension)) { // if(!isPureGhostNode(i)) bbox_local += pos; } if (this->is_distributed) { bbox = bbox_local.allSum(*communicator); } else { bbox = bbox_local; } size = bbox.size(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_not_defined, _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) { continue; } auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [&](UInt l) -> UInt { return this->getNodeGlobalId(l); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name, const std::string & group_name) { if (group_name == "all") { return this->getDumper(dumper_name); } return element_groups[group_name]->getDumper(dumper_name); } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & arrays) { ElementTypeMap nb_data_per_elem; for (auto type : arrays.elementTypes(_element_kind = _ek_not_defined)) { UInt nb_elements = this->getNbElement(type); auto & array = arrays(type); nb_data_per_elem(type) = array.getNbComponent() * array.size(); nb_data_per_elem(type) /= nb_elements; } return nb_data_per_elem; } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array); /* -------------------------------------------------------------------------- */ template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind) { std::shared_ptr field; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind); template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind); /* -------------------------------------------------------------------------- */ void Mesh::distributeImpl( Communicator & communicator, const std::function & edge_weight_function [[gnu::unused]], const std::function & vertex_weight_function [[gnu::unused]]) { AKANTU_DEBUG_ASSERT(is_distributed == false, "This mesh is already distribute"); this->communicator = &communicator; this->element_synchronizer = std::make_unique( *this, this->getID() + ":element_synchronizer", true); this->node_synchronizer = std::make_unique( *this, this->getID() + ":node_synchronizer", true); Int psize = this->communicator->getNbProc(); if (psize > 1) { #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize, edge_weight_function, vertex_weight_function); MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition); } else { MeshUtilsDistribution::distributeMeshCentralized(*this, 0); } + #else - if (psize > 1) { - AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool"); - } + AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool"); #endif } // if (psize > 1) this->is_distributed = true; this->computeBoundingBox(); + + MeshIsDistributedEvent event(*this, AKANTU_CURRENT_FUNCTION); + this->sendEvent(event); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) const { for (const auto & node : node_list) { for (const auto & element : *nodes_to_elements[node]) { elements.push_back(element); } } } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const UInt & node, Array & elements) const { for (const auto & element : *nodes_to_elements[node]) { elements.push_back(element); } } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements(UInt dimension) { Element e; UInt nb_nodes = nodes->size(); this->nodes_to_elements.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) { this->nodes_to_elements[n]->clear(); } else { this->nodes_to_elements[n] = std::make_unique>(); } } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); auto connectivity = connectivities(type, ghost_type); auto conn_it = connectivity.begin(connectivity.getNbComponent()); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (auto node : conn) { nodes_to_elements[node]->insert(e); } } } } } /* -------------------------------------------------------------------------- */ std::tuple Mesh::updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { if (global_data_updater) { return this->global_data_updater->updateData(nodes_event, elements_event); } return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } /* -------------------------------------------------------------------------- */ void Mesh::registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->global_data_updater = std::move(global_data_updater); } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this, "new_numbering", AKANTU_CURRENT_FUNCTION); auto & remove_list = event.getList(); auto & new_numbering = event.getNewNumbering(); for (auto && el : elements) { if (el.ghost_type != _not_ghost) { auto & count = ghosts_counters(el); --count; if (count > 0) { continue; } } remove_list.push_back(el); if (not new_numbering.exists(el.type, el.ghost_type)) { auto nb_element = mesh.getNbElement(el.type, el.ghost_type); auto & numbering = new_numbering.alloc(nb_element, 1, el.type, el.ghost_type); for (auto && pair : enumerate(numbering)) { std::get<1>(pair) = std::get<0>(pair); } } new_numbering(el) = UInt(-1); } auto find_last_not_deleted = [](auto && array, Int start) -> Int { do { --start; } while (start >= 0 and array[start] == UInt(-1)); return start; }; auto find_first_deleted = [](auto && array, Int start) -> Int { auto begin = array.begin(); auto it = std::find_if(begin + start, array.end(), [](auto & el) { return el == UInt(-1); }); return Int(it - begin); }; for (auto ghost_type : ghost_types) { for (auto type : new_numbering.elementTypes(_ghost_type = ghost_type)) { auto & numbering = new_numbering(type, ghost_type); auto last_not_delete = find_last_not_deleted(numbering, numbering.size()); if (last_not_delete < 0) { continue; } auto pos = find_first_deleted(numbering, 0); while (pos < last_not_delete) { std::swap(numbering[pos], numbering[last_not_delete]); last_not_delete = find_last_not_deleted(numbering, last_not_delete); pos = find_first_deleted(numbering, pos + 1); } } } this->sendEvent(event); } } // namespace akantu diff --git a/src/mesh/mesh_events.hh b/src/mesh/mesh_events.hh index bceb620f7..38422f7fc 100644 --- a/src/mesh/mesh_events.hh +++ b/src/mesh/mesh_events.hh @@ -1,206 +1,228 @@ /** * @file mesh_events.hh * * @author Nicolas Richart * * @date creation: Fri Feb 20 2015 * @date last modification: Thu Feb 20 2020 * * @brief Classes corresponding to mesh events type * * * @section LICENSE * * Copyright (©) 2015-2021 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 #include "aka_array.hh" #include "element.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MESH_EVENTS_HH_ #define AKANTU_MESH_EVENTS_HH_ namespace akantu { /// akantu::MeshEvent is the base event for meshes template class MeshEvent { public: MeshEvent(const std::string & origin = "") : origin_(origin) {} virtual ~MeshEvent() = default; /// Get the list of entity modified by the event nodes or elements const Array & getList() const { return list; } /// Get the list of entity modified by the event nodes or elements Array & getList() { return list; } std::string origin() const { return origin_; } protected: Array list; private: std::string origin_; }; class Mesh; /// akantu::MeshEvent related to new nodes in the mesh class NewNodesEvent : public MeshEvent { public: NewNodesEvent(const std::string & origin = "") : MeshEvent(origin) {} ~NewNodesEvent() override = default; }; /// akantu::MeshEvent related to nodes removed from the mesh class RemovedNodesEvent : public MeshEvent { public: inline RemovedNodesEvent(const Mesh & mesh, const std::string & origin = ""); ~RemovedNodesEvent() override = default; /// Get the new numbering following suppression of nodes from nodes arrays AKANTU_GET_MACRO_NOT_CONST(NewNumbering, new_numbering, Array &); /// Get the new numbering following suppression of nodes from nodes arrays AKANTU_GET_MACRO(NewNumbering, new_numbering, const Array &); private: Array new_numbering; }; /// akantu::MeshEvent related to new elements in the mesh class NewElementsEvent : public MeshEvent { public: NewElementsEvent(const std::string & origin = "") : MeshEvent(origin) {} ~NewElementsEvent() override = default; }; +/// akantu::MeshEvent related to the case the mesh is made distributed. +/// Note that the `list` has no meaning for this event. +class MeshIsDistributedEvent : public MeshEvent { +public: + MeshIsDistributedEvent(const Mesh & mesh, const std::string & origin = "") + : MeshEvent(origin), mesh(mesh) {} + ~MeshIsDistributedEvent() override = default; + + const Mesh & getMesh() const noexcept { return this->mesh; } + +private: + const Mesh & mesh; +}; + /// akantu::MeshEvent related to elements removed from the mesh class RemovedElementsEvent : public MeshEvent { public: inline RemovedElementsEvent(const Mesh & mesh, const ID & new_numbering_id = "new_numbering", const std::string & origin = ""); ~RemovedElementsEvent() override = default; /// Get the new numbering following suppression of elements from elements /// arrays AKANTU_GET_MACRO(NewNumbering, new_numbering, const ElementTypeMapArray &); /// Get the new numbering following suppression of elements from elements /// arrays AKANTU_GET_MACRO_NOT_CONST(NewNumbering, new_numbering, ElementTypeMapArray &); /// Get the new numbering following suppression of elements from elements /// arrays AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NewNumbering, new_numbering, UInt); /// Get the new numbering following suppression of elements from elements /// arrays AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NewNumbering, new_numbering, UInt); protected: ElementTypeMapArray new_numbering; }; /// akantu::MeshEvent for element that changed in some sort, can be seen as a /// combination of removed and added elements class ChangedElementsEvent : public RemovedElementsEvent { public: inline ChangedElementsEvent( const Mesh & mesh, const ID & new_numbering_id = "changed_event:new_numbering", const std::string & origin = "") : RemovedElementsEvent(mesh, new_numbering_id, origin) {} ~ChangedElementsEvent() override = default; AKANTU_GET_MACRO(ListOld, list, const Array &); AKANTU_GET_MACRO_NOT_CONST(ListOld, list, Array &); AKANTU_GET_MACRO(ListNew, new_list, const Array &); AKANTU_GET_MACRO_NOT_CONST(ListNew, new_list, Array &); protected: Array new_list; }; /* -------------------------------------------------------------------------- */ class MeshEventHandler { public: virtual ~MeshEventHandler() = default; /* ------------------------------------------------------------------------ */ /* Internal code */ /* ------------------------------------------------------------------------ */ private: /// send a akantu::NewNodesEvent inline void sendEvent(const NewNodesEvent & event) { onNodesAdded(event.getList(), event); } /// send a akantu::RemovedNodesEvent inline void sendEvent(const RemovedNodesEvent & event) { onNodesRemoved(event.getList(), event.getNewNumbering(), event); } /// send a akantu::NewElementsEvent inline void sendEvent(const NewElementsEvent & event) { onElementsAdded(event.getList(), event); } /// send a akantu::RemovedElementsEvent inline void sendEvent(const RemovedElementsEvent & event) { onElementsRemoved(event.getList(), event.getNewNumbering(), event); } /// send a akantu::ChangedElementsEvent inline void sendEvent(const ChangedElementsEvent & event) { onElementsChanged(event.getListOld(), event.getListNew(), event.getNewNumbering(), event); } + /// send a akantu::MeshIsDistributedEvent + inline void sendEvent(const MeshIsDistributedEvent & event) { + onMeshIsDistributed(event.getMesh(), event); + } template friend class EventHandlerManager; /* ------------------------------------------------------------------------ */ /* Interface */ /* ------------------------------------------------------------------------ */ public: /// function to implement to react on akantu::NewNodesEvent virtual void onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent & /*event*/) {} /// function to implement to react on akantu::RemovedNodesEvent virtual void onNodesRemoved(const Array & /*nodes_list*/, const Array & /*new_numbering*/, const RemovedNodesEvent & /*event*/) {} /// function to implement to react on akantu::NewElementsEvent virtual void onElementsAdded(const Array & /*elements_list*/, const NewElementsEvent & /*event*/) {} /// function to implement to react on akantu::RemovedElementsEvent virtual void onElementsRemoved(const Array & /*elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const RemovedElementsEvent & /*event*/) {} /// function to implement to react on akantu::ChangedElementsEvent virtual void onElementsChanged(const Array & /*old_elements_list*/, const Array & /*new_elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const ChangedElementsEvent & /*event*/) {} + + /// function to implement to react on akantu::MeshIsDistributedEvent + virtual void onMeshIsDistributed(const Mesh & /*mesh*/, + const MeshIsDistributedEvent & /*event*/) {} }; } // namespace akantu #endif /* AKANTU_MESH_EVENTS_HH_ */ diff --git a/src/model/common/dof_manager/dof_manager.cc b/src/model/common/dof_manager/dof_manager.cc index 137192f11..e8e351d26 100644 --- a/src/model/common/dof_manager/dof_manager.cc +++ b/src/model/common/dof_manager/dof_manager.cc @@ -1,1017 +1,1041 @@ /** * @file dof_manager.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Sat Mar 06 2021 * * @brief Implementation of the common parts of the DOFManagers * * * @section LICENSE * * Copyright (©) 2015-2021 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 "dof_manager.hh" #include "communicator.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver.hh" #include "periodic_node_synchronizer.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id) : id(id), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(Communicator::getStaticCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id) : id(id), mesh(&mesh), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ std::vector DOFManager::getDOFIDs() const { std::vector keys; for (const auto & dof_data : this->dofs) { keys.push_back(dof_data.first); } return keys; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_vect.size() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); Array::const_matrix_iterator elem_it = elementary_vect.begin(nb_degree_of_freedom, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++elem_it) { UInt element = el; if (filter_it != nullptr) { // conn_it = conn_begin + *filter_it; element = *filter_it; } // const Vector & conn = *conn_it; const Matrix & elemental_val = *elem_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_node = connectivity(element, n) * nb_degree_of_freedom; Vector assemble(array_assembeled.storage() + offset_node, nb_degree_of_freedom); Vector elem_val = elemental_val(n); assemble.aXplusY(elem_val, scale_factor); } if (filter_it != nullptr) { ++filter_it; } // else // ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToResidual(dof_id, array_localy_assembeled, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToLumpedMatrix(dof_id, array_localy_assembeled, lumped_mtx, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = *pair.second; this->assembleMatMulVectToResidual(dof_id, A_id, *dof_data.dof, scale_factor); } } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { for (auto && data : this->dofs) { auto & dof_data = *data.second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(data.first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ void DOFManager::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getSolution(), solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); this->assembleToGlobalArray(dof_id, array_to_assemble, this->getResidual(), scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); auto & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number"), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ template auto DOFManager::countDOFsForNodes(const DOFData & dof_data, UInt nb_nodes, Func && getNode) { auto nb_local_dofs = nb_nodes; decltype(nb_local_dofs) nb_pure_local = 0; for (auto n : arange(nb_nodes)) { UInt node = getNode(n); // http://www.open-std.org/jtc1/sc22/open/n2356/conv.html // bool are by convention casted to 0 and 1 when promoted to int nb_pure_local += this->mesh->isLocalOrMasterNode(node); nb_local_dofs -= this->mesh->isPeriodicSlave(node); } const auto & dofs_array = *dof_data.dof; nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); return std::make_pair(nb_local_dofs, nb_pure_local); } /* -------------------------------------------------------------------------- */ auto DOFManager::getNewDOFDataInternal(const ID & dof_id) -> DOFData & { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dof_data_ptr = this->getNewDOFData(dof_id); DOFData & dof_data = *dof_data_ptr; this->dofs[dof_id] = std::move(dof_data_ptr); return dof_data; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ std::tuple DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dof_data = this->getDOFData(dof_id); dof_data.dof = &dofs_array; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; const auto & support_type = dof_data.support_type; switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; std::function getNode; if (group == "__mesh__") { AKANTU_DEBUG_ASSERT( dofs_array.size() == this->mesh->getNbNodes(), "The array of dof is too short to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes( dof_data, this->mesh->getNbNodes(), [](auto && n) { return n; }); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); AKANTU_DEBUG_ASSERT( dofs_array.size() == node_group.size(), "The array of dof is too shot to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes(dof_data, node_group.size(), [&node_group](auto && n) { return node_group(n); }); } break; } case _dst_generic: { nb_local_dofs = nb_pure_local = dofs_array.size() * dofs_array.getNbComponent(); break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } dof_data.local_nb_dofs = nb_local_dofs; dof_data.pure_local_nb_dofs = nb_pure_local; dof_data.ghosts_nb_dofs = nb_local_dofs - nb_pure_local; this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; auto nb_total_pure_local = nb_pure_local; communicator.allReduce(nb_total_pure_local, SynchronizerOperation::_sum); this->system_size += nb_total_pure_local; // updating the dofs data after counting is finished switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; if (group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_local_dofs, nb_pure_local, support_nodes.size(), [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local, mesh->getNbNodes(), [](UInt node) -> UInt { return node; }); } break; } case _dst_generic: { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local); break; } } return std::make_tuple(nb_local_dofs, nb_pure_local, nb_total_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsPrevious(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.previous != nullptr) { AKANTU_EXCEPTION("The previous dofs array for " << dof_id << " has already been registered"); } dof.previous = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsIncrement(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.increment != nullptr) { AKANTU_EXCEPTION("The dofs increment array for " << dof_id << " has already been registered"); } dof.increment = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative) { DOFData & dof = this->getDOFData(dof_id); std::vector *> & derivatives = dof.dof_derivatives; if (derivatives.size() < order) { derivatives.resize(order, nullptr); } else { if (derivatives[order - 1] != nullptr) { AKANTU_EXCEPTION("The dof derivatives of order " << order << " already been registered for this dof (" << dof_id << ")"); } } derivatives[order - 1] = &dofs_derivative; } /* -------------------------------------------------------------------------- */ void DOFManager::registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs) { DOFData & dof = this->getDOFData(dof_id); if (dof.blocked_dofs != nullptr) { AKANTU_EXCEPTION("The blocked dofs array for " << dof_id << " has already been registered"); } dof.blocked_dofs = &blocked_dofs; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix SolverVector & DOFManager::registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->lumped_matrices.find(matrix_id); if (it != this->lumped_matrices.end()) { AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->lumped_matrices[matrix_id] = std::move(matrix); ret.resize(); return ret; } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::registerNonLinearSolver( const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver) { NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it != this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " already exists in " << this->id); } NonLinearSolver & ret = *non_linear_solver; this->non_linear_solvers[non_linear_solver_id] = std::move(non_linear_solver); return ret; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::registerTimeStepSolver( const ID & time_step_solver_id, std::unique_ptr & time_step_solver) { TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it != this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " already exists in " << this->id); } TimeStepSolver & ret = *time_step_solver; this->time_step_solvers[time_step_solver_id] = std::move(time_step_solver); return ret; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::getMatrix(const ID & id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it == this->matrices.end()) { AKANTU_SILENT_EXCEPTION("The matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasMatrix(const ID & id) const { ID mtx_id = this->id + ":mtx:" + id; auto it = this->matrices.find(mtx_id); return it != this->matrices.end(); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManager::getLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const SolverVector & DOFManager::getLumpedMatrix(const ID & id) const { ID matrix_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasLumpedMatrix(const ID & id) const { ID mtx_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(mtx_id); return it != this->lumped_matrices.end(); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::getNonLinearSolver(const ID & id) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it == this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasNonLinearSolver(const ID & id) const { ID solver_id = this->id + ":nls:" + id; auto it = this->non_linear_solvers.find(solver_id); return it != this->non_linear_solvers.end(); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::getTimeStepSolver(const ID & id) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it == this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasTimeStepSolver(const ID & solver_id) const { ID time_step_solver_id = this->id + ":tss:" + solver_id; auto it = this->time_step_solvers.find(time_step_solver_id); return it != this->time_step_solvers.end(); } /* -------------------------------------------------------------------------- */ void DOFManager::savePreviousDOFs(const ID & dofs_id) { this->getPreviousDOFs(dofs_id).copy(this->getDOFs(dofs_id)); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroResidual() { this->residual->zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroMatrix(const ID & mtx) { this->getMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ std::pair DOFManager::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { auto & dof_data = this->getDOFData(dof_id); UInt nb_new_local_dofs; UInt nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = countDOFsForNodes(dof_data, nodes_list.size(), [&nodes_list](auto && n) { return nodes_list(n); }); this->pure_local_system_size += nb_new_pure_local; this->local_system_size += nb_new_local_dofs; UInt nb_new_global = nb_new_pure_local; communicator.allReduce(nb_new_global, SynchronizerOperation::_sum); this->system_size += nb_new_global; dof_data.solution.resize(local_system_size); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, nodes_list.size(), [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::resizeGlobalArrays() { // resize all relevant arrays this->residual->resize(); this->solution->resize(); this->data_cache->resize(); for (auto & lumped_matrix : lumped_matrices) { lumped_matrix.second->resize(); } for (auto & matrix : matrices) { matrix.second->clearProfile(); } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, const NewNodesEvent & /*unused*/) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) { continue; } const auto & group = dof_data.group_support; if (group == "__mesh__") { this->updateNodalDOFs(dof_id, nodes_list); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup(); Array new_nodes_list; for (const auto & node : nodes_list) { if (node_group.find(node) != UInt(-1)) { new_nodes_list.push_back(node); } } this->updateNodalDOFs(dof_id, new_nodes_list); } } this->resizeGlobalArrays(); } +/* -------------------------------------------------------------------------- */ +void DOFManager::onMeshIsDistributed(const Mesh & mesh_, + const MeshIsDistributedEvent & event) { + AKANTU_DEBUG_ASSERT(this->mesh != nullptr, "The `Mesh` pointer is not set."); + + // check if the distributed state of the residual and the mesh are the same. + if (this->mesh->isDistributed() != this->residual->isDistributed()) { + // TODO: Allow to reallocate the internals, in that case one could actually + // react on that event. + auto is_or_is_not = [](bool q) { + return ((q) ? std::string("is") : std::string("is not")); + }; + AKANTU_EXCEPTION("There is an inconsistency about the distribution state " + "of the `DOFManager`." + " It seams that the `Mesh` " + << is_or_is_not(this->mesh->isDistributed()) + << " distributed, but the `DOFManager`'s residual " + << is_or_is_not(this->residual->isDistributed()) + << ", which is of type " + << debug::demangle(typeid(*this->residual).name()) << "."); + + } +} + /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor(DOFManager::DOFData & dof_data, DOFManager & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { for (auto && pair : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { UInt node; Int dof; std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(Int); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node.at(node); for (const auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node[node]; for (const auto & dof : dofs) { Int global_dof; buffer >> global_dof; AKANTU_DEBUG_ASSERT( (dof_manager.global_equation_number(dof) == -1 or dof_manager.global_equation_number(dof) == global_dof), "This dof already had a global_dof_id which is different from " "the received one. " << dof_manager.global_equation_number(dof) << " != " << global_dof); dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManager::DOFData & dof_data; DOFManager & dof_manager; }; /* -------------------------------------------------------------------------- */ auto DOFManager::computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local) { // determine the first local/global dof id to use UInt offset = 0; this->communicator.exclusiveScan(nb_new_pure_local, offset); auto first_global_dof_id = this->first_global_dof_id + offset; auto first_local_dof_id = this->local_system_size - nb_new_local_dofs; offset = nb_new_pure_local; this->communicator.allReduce(offset); this->first_global_dof_id += offset; return std::make_pair(first_local_dof_id, first_global_dof_id); } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, UInt nb_node, const std::function & getNode) { auto nb_local_dofs_added = nb_node * dof_data.dof->getNbComponent(); auto first_dof_pos = dof_data.local_equation_number.size(); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_local_dofs_added); dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_local_dofs_added); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); std::unordered_map, UInt> masters_dofs; // update per dof info UInt local_eq_num; UInt first_global_dof_id; std::tie(local_eq_num, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); for (auto d : arange(nb_local_dofs_added)) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof_flag = this->mesh->getNodeFlag(node); dof_data.associated_nodes.push_back(node); auto is_local_dof = this->mesh->isLocalOrMasterNode(node); auto is_periodic_slave = this->mesh->isPeriodicSlave(node); auto is_periodic_master = this->mesh->isPeriodicMaster(node); if (is_periodic_slave) { dof_data.local_equation_number.push_back(-1); continue; } // update equation numbers this->dofs_flag(local_eq_num) = dof_flag; dof_data.local_equation_number.push_back(local_eq_num); if (is_local_dof) { this->global_equation_number(local_eq_num) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = local_eq_num; ++first_global_dof_id; } else { this->global_equation_number(local_eq_num) = -1; } if (is_periodic_master) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof = d % dof_data.dof->getNbComponent(); masters_dofs.insert( std::make_pair(std::make_pair(node, dof), local_eq_num)); } ++local_eq_num; } // correct periodic slave equation numbers if (this->mesh->isPeriodic()) { auto assoc_begin = dof_data.associated_nodes.begin(); for (auto d : arange(nb_local_dofs_added)) { auto node = dof_data.associated_nodes(first_dof_pos + d); if (not this->mesh->isPeriodicSlave(node)) { continue; } auto master_node = this->mesh->getPeriodicMaster(node); auto dof = d % dof_data.dof->getNbComponent(); dof_data.local_equation_number(first_dof_pos + d) = masters_dofs[std::make_pair(master_node, dof)]; } } // synchronize the global numbering for slaves nodes if (this->mesh->isDistributed()) { GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); if (this->mesh->isPeriodic()) { mesh->getPeriodicNodeSynchronizer().synchronizeOnce( data_accessor, SynchronizationTag::_giu_global_conn); } auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, SynchronizationTag::_ask_nodes); } } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local) { dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); UInt first_local_dof_id; UInt first_global_dof_id; std::tie(first_local_dof_id, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); // update per dof info for (auto _ [[gnu::unused]] : arange(nb_new_local_dofs)) { // update equation numbers this->dofs_flag(first_local_dof_id) = NodeFlag::_normal; dof_data.local_equation_number.push_back(first_local_dof_id); this->global_equation_number(first_local_dof_id) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = first_local_dof_id; ++first_global_dof_id; ++first_local_dof_id; } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array & /*unused*/, const Array & /*unused*/, const RemovedNodesEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const RemovedElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::updateGlobalBlockedDofs() { this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.reserve(this->local_system_size, 0); this->previous_global_blocked_dofs_release = this->global_blocked_dofs_release; for (auto & pair : dofs) { if (not this->hasBlockedDOFs(pair.first)) { continue; } DOFData & dof_data = *pair.second; for (auto && data : zip(dof_data.getLocalEquationsNumbers(), make_view(*dof_data.blocked_dofs))) { const auto & dof = std::get<0>(data); const auto & is_blocked = std::get<1>(data); if (is_blocked) { this->global_blocked_dofs.push_back(dof); } } } std::sort(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); auto last = std::unique(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); this->global_blocked_dofs.resize(last - this->global_blocked_dofs.begin()); auto are_equal = global_blocked_dofs.size() == previous_global_blocked_dofs.size() and std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) { ++this->global_blocked_dofs_release; } } /* -------------------------------------------------------------------------- */ void DOFManager::applyBoundary(const ID & matrix_id) { auto & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getRelease()) { if (this->hasBlockedDOFsChanged()) { J.applyBoundary(); } previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getRelease(); this->previous_global_blocked_dofs_release = this->global_blocked_dofs_release; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor) { auto & A = this->getMatrix(A_id); data_cache->resize(); data_cache->zero(); this->assembleToGlobalArray(dof_id, x, *data_cache, 1.); A.matVecMul(*data_cache, array, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { assembleMatMulVectToGlobalArray(dof_id, A_id, x, *residual, scale_factor); } } // namespace akantu diff --git a/src/model/common/dof_manager/dof_manager.hh b/src/model/common/dof_manager/dof_manager.hh index b9ccc7b99..73c9bbf98 100644 --- a/src/model/common/dof_manager/dof_manager.hh +++ b/src/model/common/dof_manager/dof_manager.hh @@ -1,715 +1,720 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Fri Jul 24 2020 * * @brief Class handling the different types of dofs * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_DOF_MANAGER_HH_ #define AKANTU_DOF_MANAGER_HH_ namespace akantu { class TermsToAssemble; class NonLinearSolver; class TimeStepSolver; class SparseMatrix; class SolverVector; class SolverCallback; } // namespace akantu namespace akantu { class DOFManager : protected MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ protected: struct DOFData; public: DOFManager(const ID & id = "dof_manager"); DOFManager(Mesh & mesh, const ID & id = "dof_manager"); ~DOFManager() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register an array of degree of freedom virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// the dof as an implied type of _dst_nodal and is defined only on a subset /// of nodes virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group); /// register an array of previous values of the degree of freedom virtual void registerDOFsPrevious(const ID & dof_id, Array & dofs_array); /// register an array of increment of degree of freedom virtual void registerDOFsIncrement(const ID & dof_id, Array & dofs_array); /// register an array of derivatives for a particular dof array virtual void registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative); /// register array representing the blocked degree of freedoms virtual void registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor = 1.); /// Assemble an array to the global lumped matrix array virtual void assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.); /** * Assemble elementary values to a local array of the size nb_nodes * * nb_dof_per_node. The dof number is implicitly considered as * conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to a global array corresponding to a lumped * matrix */ virtual void assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. With 0 < * n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type = _not_ghost, const MatrixType & elemental_matrix_type = _symmetric, const Array & filter_elements = empty_filter) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToArray(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor = 1) = 0; /// multiply a vector by a lumped matrix and assemble the result to the /// residual virtual void assembleLumpedMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) = 0; /// assemble coupling terms between to dofs virtual void assemblePreassembledMatrix(const ID & matrix_id, const TermsToAssemble & terms) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1); /// multiply the dofs by a matrix and assemble the result to the residual virtual void assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor = 1); /// updates the global blocked_dofs array virtual void updateGlobalBlockedDofs(); /// sets the residual to 0 virtual void zeroResidual(); /// sets the matrix to 0 virtual void zeroMatrix(const ID & mtx); /// sets the lumped matrix to 0 virtual void zeroLumpedMatrix(const ID & mtx); virtual void applyBoundary(const ID & matrix_id = "J"); // virtual void applyBoundaryLumped(const ID & matrix_id = "J"); /// extract a lumped matrix part corresponding to a given dof virtual void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped); /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); private: /// dispatch the creation of the dof data and register it DOFData & getNewDOFDataInternal(const ID & dof_id); protected: /// common function to help registering dofs the return values are the add new /// numbers of local dofs, pure local dofs, and system size virtual std::tuple registerDOFsInternal(const ID & dof_id, Array & dofs_array); /// minimum functionality to implement per derived version of the DOFManager /// to allow the splitSolutionPerDOFs function to work virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array); /// fill a Vector with the equation numbers corresponding to the given /// connectivity static inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number); /// Assemble a array to a global one void assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor = 1.); /// common function that can be called by derived class with proper matrice /// types template void assemblePreassembledMatrix_(Mat & A, const TermsToAssemble & terms); template void assembleElementalMatricesToMatrix_(Mat & A, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements); template void assembleMatMulVectToArray_(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isSlaveDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isPureGhostDOF(UInt local_dof_num); /// tells if the dof manager knows about a global dof bool hasGlobalEquationNumber(Int global) const; /// return the local index of the global equation number inline Int globalToLocalEquationNumber(Int global) const; /// converts local equation numbers to global equation numbers; inline Int localToGlobalEquationNumber(Int local) const; /// get the array of dof types (use only if you know what you do...) inline NodeFlag getDOFFlag(Int local_id) const; /// defines if the boundary changed bool hasBlockedDOFsChanged() const { return this->global_blocked_dofs_release != this->previous_global_blocked_dofs_release; } /// Global number of dofs AKANTU_GET_MACRO(SystemSize, this->system_size, UInt); /// Local number of dofs AKANTU_GET_MACRO(LocalSystemSize, this->local_system_size, UInt); /// Pure local number of dofs AKANTU_GET_MACRO(PureLocalSystemSize, this->pure_local_system_size, UInt); /// Retrieve all the registered DOFs std::vector getDOFIDs() const; /* ------------------------------------------------------------------------ */ /* DOFs and derivatives accessors */ /* ------------------------------------------------------------------------ */ /// Get a reference to the registered dof array for a given id inline Array & getDOFs(const ID & dofs_id); /// Get the support type of a given dof inline DOFSupportType getSupportType(const ID & dofs_id) const; /// are the dofs registered inline bool hasDOFs(const ID & dof_id) const; /// Get a reference to the registered dof derivatives array for a given id inline Array & getDOFsDerivatives(const ID & dofs_id, UInt order); /// Does the dof has derivatives inline bool hasDOFsDerivatives(const ID & dofs_id, UInt order) const; /// Get a reference to the blocked dofs array registered for the given id inline const Array & getBlockedDOFs(const ID & dofs_id) const; /// Does the dof has a blocked array inline bool hasBlockedDOFs(const ID & dofs_id) const; /// Get a reference to the registered dof increment array for a given id inline Array & getDOFsIncrement(const ID & dofs_id); /// Does the dof has a increment array inline bool hasDOFsIncrement(const ID & dofs_id) const; /// Does the dof has a previous array inline Array & getPreviousDOFs(const ID & dofs_id); /// Get a reference to the registered dof array for previous step values a /// given id inline bool hasPreviousDOFs(const ID & dofs_id) const; /// saves the values from dofs to previous dofs virtual void savePreviousDOFs(const ID & dofs_id); /// Get a reference to the solution array registered for the given id inline const Array & getSolution(const ID & dofs_id) const; /// Get a reference to the solution array registered for the given id inline Array & getSolution(const ID & dofs_id); /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(PreviousGlobalBlockedDOFs, previous_global_blocked_dofs, const Array &); /* ------------------------------------------------------------------------ */ /* Matrices accessors */ /* ------------------------------------------------------------------------ */ /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) = 0; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) = 0; /// Get the equation numbers corresponding to a dof_id. This might be used to /// access the matrix. inline const Array & getLocalEquationsNumbers(const ID & dof_id) const; protected: /// get the array of dof types (use only if you know what you do...) inline const Array & getDOFsAssociatedNodes(const ID & dof_id) const; protected: /* ------------------------------------------------------------------------ */ /// register a matrix SparseMatrix & registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a lumped matrix (aka a Vector) SolverVector & registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a non linear solver instantiated by a derived class NonLinearSolver & registerNonLinearSolver(const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver); /// register a time step solver instantiated by a derived class TimeStepSolver & registerTimeStepSolver(const ID & time_step_solver_id, std::unique_ptr & time_step_solver); template NonLinearSolver & registerNonLinearSolver(DMType & dm, const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls = std::make_unique(dm, type, non_linear_solver_id); return this->registerNonLinearSolver(non_linear_solver_id, nls); } template TimeStepSolver & registerTimeStepSolver(DMType & dm, const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique( dm, type, non_linear_solver, solver_callback, time_step_solver_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } template SparseMatrix & registerSparseMatrix(DMType & dm, const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(dm, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } template SparseMatrix & registerSparseMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; auto & sm_to_copy = aka::as_type(this->getMatrix(matrix_to_copy_id)); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } template SolverVector & registerLumpedMatrix(DMType & dm, const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; std::unique_ptr sm = std::make_unique(dm, matrix_id); return this->registerLumpedMatrix(matrix_id, sm); } protected: virtual void makeConsistentForPeriodicity(const ID & dof_id, SolverVector & array) = 0; virtual void assembleToGlobalArray(const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array, Real scale_factor) = 0; public: /// extract degrees of freedom (identified by ID) from a global solver array virtual void getArrayPerDOFs(const ID & dof_id, const SolverVector & global, Array & local) = 0; /// Get the reference of an existing matrix SparseMatrix & getMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasMatrix(const ID & matrix_id) const; /// Get an instance of a new lumped matrix virtual SolverVector & getNewLumpedMatrix(const ID & matrix_id) = 0; /// Get the lumped version of a given matrix const SolverVector & getLumpedMatrix(const ID & matrix_id) const; /// Get the lumped version of a given matrix SolverVector & getLumpedMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasLumpedMatrix(const ID & matrix_id) const; /* ------------------------------------------------------------------------ */ /* Non linear system solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) = 0; /// get instance of a non linear solver virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id); /// check if the given solver exists bool hasNonLinearSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver virtual TimeStepSolver & getNewTimeStepSolver(const ID & time_step_solver_id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) = 0; /// get instance of a time step solver virtual TimeStepSolver & getTimeStepSolver(const ID & time_step_solver_id); /// check if the given solver exists bool hasTimeStepSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ const Mesh & getMesh() { if (mesh != nullptr) { return *mesh; } AKANTU_EXCEPTION("No mesh registered in this dof manager"); } /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Communicator, communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, communicator, auto &); /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Solution, *(solution.get()), const auto &); AKANTU_GET_MACRO_NOT_CONST(Solution, *(solution.get()), auto &); AKANTU_GET_MACRO(Residual, *(residual.get()), const auto &); AKANTU_GET_MACRO_NOT_CONST(Residual, *(residual.get()), auto &); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: friend class GlobalDOFInfoDataAccessor; /// helper function for the DOFManager::onNodesAdded method virtual std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list); template auto countDOFsForNodes(const DOFData & dof_data, UInt nb_nodes, Func && getNode); void updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, UInt nb_nodes, const std::function & getNode); void updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local); auto computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local); /// resize all the global information and takes the needed measure like /// cleaning matrices profiles virtual void resizeGlobalArrays(); public: /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) override; /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array & elements_list, const NewElementsEvent & event) override; /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; + /// function to implement to react on akantu::MeshIsDistributedEvent + void + onMeshIsDistributed(const Mesh& mesh, + const MeshIsDistributedEvent & event) override; + protected: inline DOFData & getDOFData(const ID & dof_id); inline const DOFData & getDOFData(const ID & dof_id) const; template inline DOFData_ & getDOFDataTyped(const ID & dof_id); template inline const DOFData_ & getDOFDataTyped(const ID & dof_id) const; virtual std::unique_ptr getNewDOFData(const ID & dof_id) = 0; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// dof representations in the dof manager struct DOFData { DOFData() = delete; explicit DOFData(const ID & dof_id); virtual ~DOFData(); /// DOF support type (nodal, general) this is needed to determine how the /// dof are shared among processors DOFSupportType support_type; ID group_support; /// Degree of freedom array Array * dof{nullptr}; /// Blocked degree of freedoms array Array * blocked_dofs{nullptr}; /// Degree of freedoms increment Array * increment{nullptr}; /// Degree of freedoms at previous step Array * previous{nullptr}; /// Solution associated to the dof Array solution; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; /* ---------------------------------------------------------------------- */ /// number of dofs to consider locally for this dof id UInt local_nb_dofs{0}; /// Number of purely local dofs UInt pure_local_nb_dofs{0}; /// number of ghost dofs UInt ghosts_nb_dofs{0}; /// local numbering equation numbers Array local_equation_number; /// associated node for _dst_nodal dofs only Array associated_nodes; virtual Array & getLocalEquationsNumbers() { return local_equation_number; } }; /// type to store dofs information using DOFStorage = std::map>; /// type to store all the matrices using SparseMatricesMap = std::map>; /// type to store all the lumped matrices using LumpedMatricesMap = std::map>; /// type to store all the non linear solver using NonLinearSolversMap = std::map>; /// type to store all the time step solver using TimeStepSolversMap = std::map>; ID id; /// store a reference to the dof arrays DOFStorage dofs; /// list of sparse matrices that where created SparseMatricesMap matrices; /// list of lumped matrices LumpedMatricesMap lumped_matrices; /// non linear solvers storage NonLinearSolversMap non_linear_solvers; /// time step solvers storage TimeStepSolversMap time_step_solvers; /// reference to the underlying mesh Mesh * mesh{nullptr}; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size{0}; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size{0}; /// Total number of degrees of freedom UInt system_size{0}; /// rhs to the system of equation corresponding to the residual linked to the /// different dofs std::unique_ptr residual; /// solution of the system of equation corresponding to the different dofs std::unique_ptr solution; /// a vector that helps internally to perform some tasks std::unique_ptr data_cache; /// define the dofs type, local, shared, ghost Array dofs_flag; /// equation number in global numbering Array global_equation_number; using equation_numbers_map = std::unordered_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// Communicator used for this manager, should be the same as in the mesh if a /// mesh is registered Communicator & communicator; /// accumulator to know what would be the next global id to use UInt first_global_dof_id{0}; /// Release at last apply boundary on jacobian UInt jacobian_release{0}; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; UInt global_blocked_dofs_release{0}; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array previous_global_blocked_dofs; UInt previous_global_blocked_dofs_release{0}; private: /// This is for unit testing friend class DOFManagerTester; }; using DefaultDOFManagerFactory = Factory; using DOFManagerFactory = Factory; } // namespace akantu #include "dof_manager_inline_impl.hh" #endif /* AKANTU_DOF_MANAGER_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 21dc28321..f33020eb4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1250 +1,1254 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Mauro Corrado * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Implementation of the SolidMechanicsModel class * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model * @param model_type this is an internal parameter for inheritance purposes */ SolidMechanicsModel::SolidMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, std::move(dof_manager), dim, id), material_index("material index", id), material_local_numbering("material local numbering", id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); this->mesh.registerDumper("solid_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); material_selector = std::make_shared(material_index); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper().setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param options * \parblock * contains the different options to initialize the model * \li \c analysis_method specify the type of solver to use * \endparblock */ void SolidMechanicsModel::initFullImpl(const ModelOptions & options) { material_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize the materials if (not this->parser.getLastParsedFile().empty()) { this->instantiateMaterials(); this->initMaterials(); } this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ std::tuple SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { auto & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->previous_displacement, spatial_dimension, "previous_displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForces(); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) const { // \TODO check the materials to know what is the correct answer if (matrix_id == "C") { return _mt_not_defined; } if (matrix_id == "K") { auto matrix_type = _unsymmetric; for (auto & material : materials) { matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id)); } } return _symmetric; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { this->assembleMass(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::beforeSolveStep() { for (auto & material : materials) { material->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::afterSolveStep(bool converged) { for (auto & material : materials) { material->afterSolveStep(converged); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as \f$F_{int} = \int_{\Omega} N * \sigma d\Omega@\f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->zero(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if (this->non_local_manager) { this->non_local_manager->computeAllNonLocalStresses(); } // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix(bool need_to_reassemble) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix for (auto & material : materials) { need_to_reassemble |= material->hasMatrixChanged("K"); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").zero(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) { return; } this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal(field_name, _ek_regular)) { continue; } for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { + this->assembleLumpedMatrix("M"); + auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const auto & v = *v_it; const auto & m = *m_it; Real mv2 = 0.; auto is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); auto count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) { mv2 += v(i) * v(i) * m(i); } } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { + this->assembleMatrix("M"); + Array Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension), make_view(*this->velocity, spatial_dimension))) { ekin += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } } else { AKANTU_ERROR("No function called to assemble the mass matrix."); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); auto vit = vel_on_quad.begin(Model::spatial_dimension); auto vend = vel_on_quad.end(Model::spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); decltype(ext_force_it) incr_or_velo_it; if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } else { incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); } Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) { work -= int_force(i) * incr_or_velo(i); } else { work += ext_force(i) * incr_or_velo(i); } } } } mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum); if (this->method != _static) { work *= this->getTimeStep(); } AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; for (auto & material : materials) { energy += material->getEnergy(energy_id); } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const ID & energy_id, const ID & group_id) { auto && group = mesh.getElementGroup(group_id); auto energy = 0.; for (auto && type : group.elementTypes()) { for (auto el : group.getElementsIterable(type)) { energy += getEnergy(energy_id, el); } } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->material_index.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); this->material_local_numbering.initialize( mesh, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); ElementTypeMapArray filter("new_element_filter", this->getID()); for (const auto & elem : element_list) { if (mesh.getSpatialDimension(elem.type) != spatial_dimension) { continue; } if (!filter.exists(elem.type, elem.ghost_type)) { filter.alloc(0, 1, elem.type, elem.ghost_type); } filter(elem.type, elem.ghost_type).push_back(elem.element); } // this fails in parallel if the event is sent on facet between constructor // and initFull \todo: to debug... this->assignMaterialToElements(&filter); for (auto & material : materials) { material->onElementsAdded(element_list, event); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) { mass->resize(nb_nodes, 0.); } if (velocity) { velocity->resize(nb_nodes, 0.); } if (acceleration) { acceleration->resize(nb_nodes, 0.); } if (external_force) { external_force->resize(nb_nodes, 0.); } if (internal_force) { internal_force->resize(nb_nodes, 0.); } if (blocked_dofs) { blocked_dofs->resize(nb_nodes, false); } if (current_position) { current_position->resize(nb_nodes, 0.); } if (previous_displacement) { previous_displacement->resize(nb_nodes, 0.); } if (displacement_increment) { displacement_increment->resize(nb_nodes, 0.); } for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } need_to_reassemble_lumped_mass = true; need_to_reassemble_mass = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Array & /*element_list*/, const Array & new_numbering, const RemovedNodesEvent & /*event*/) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) { mesh.removeNodesFromArray(*mass, new_numbering); } if (velocity) { mesh.removeNodesFromArray(*velocity, new_numbering); } if (acceleration) { mesh.removeNodesFromArray(*acceleration, new_numbering); } if (internal_force) { mesh.removeNodesFromArray(*internal_force, new_numbering); } if (external_force) { mesh.removeNodesFromArray(*external_force, new_numbering); } if (blocked_dofs) { mesh.removeNodesFromArray(*blocked_dofs, new_numbering); } // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) { mesh.removeNodesFromArray(*displacement_increment, new_numbering); } if (previous_displacement) { mesh.removeNodesFromArray(*previous_displacement, new_numbering); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); if (velocity) { velocity->printself(stream, indent + 2); } if (acceleration) { acceleration->printself(stream, indent + 2); } if (mass) { mass->printself(stream, indent + 2); } external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + materials [" << std::endl; for (const auto & material : materials) { material->printself(stream, indent + 2); } stream << space << " ]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeNonLocal() { this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( GhostType ghost_type) { for (auto & mat : materials) { MaterialNonLocalInterface * mat_non_local; if ((mat_non_local = dynamic_cast(mat.get())) == nullptr) { continue; } ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); for (const auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } mat_non_local->initMaterialNonLocal(); mat_non_local->insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) { for (auto & mat : materials) { if (not aka::is_of_type(*mat)) { continue; } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal(field_name, kind)) { material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { if (not aka::is_of_type(*mat)) { continue; } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return getFEEngineClassBoundary(name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::splitElementByMaterial( const Array & elements, std::vector> & elements_per_mat) const { for (const auto & el : elements) { Element mat_el = el; mat_el.element = this->material_local_numbering(el); elements_per_mat[this->material_index(el)].push_back(mat_el); } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_material_id: { size += elements.size() * sizeof(UInt); break; } case SynchronizationTag::_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case SynchronizationTag::_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case SynchronizationTag::_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case SynchronizationTag::_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { packElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case SynchronizationTag::_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { for (auto && element : elements) { UInt recv_mat_index; buffer >> recv_mat_index; UInt & mat_index = material_index(element); if (mat_index != UInt(-1)) { continue; } // add ghosts element to the correct material mat_index = recv_mat_index; UInt index = materials[mat_index]->addElement(element); material_local_numbering(element) = index; } break; } case SynchronizationTag::_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case SynchronizationTag::_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case SynchronizationTag::_smm_res: /* FALLTHRU */ case SynchronizationTag::_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case SynchronizationTag::_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 54fd44268..6c01e0cb1 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,638 +1,777 @@ /** * @file structural_mechanics_model.cc * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Mon Mar 15 2021 * * @brief Model implementation for Structural Mechanics elements * * * @section LICENSE * * Copyright (©) 2010-2021 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 "dof_manager.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "shape_structural.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #include "group_manager_inline_impl.hh" /* -------------------------------------------------------------------------- */ #include "structural_element_bernoulli_beam_2.hh" #include "structural_element_bernoulli_beam_3.hh" #include "structural_element_kirchhoff_shell.hh" /* -------------------------------------------------------------------------- */ //#include "structural_mechanics_model_inline_impl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt StructuralMechanicsModel::getNbDegreeOfFreedom(ElementType type) { UInt ndof = 0; #define GET_(type) ndof = ElementClass::getNbDegreeOfFreedom() AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural); #undef GET_ return ndof; } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id), f_m2a(1.0), stress("stress", id), element_material("element_material", id), set_ID("beam sets", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) { nb_degree_of_freedom = 3; } else if (spatial_dimension == 3) { nb_degree_of_freedom = 6; } else { AKANTU_TO_IMPLEMENT(); } this->mesh.registerDumper("structural_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); this->dumper_default_element_kind = _ek_structural; mesh.getElementalData("extra_normal") .initialize(mesh, _element_kind = _ek_structural, _nb_component = spatial_dimension, _with_nb_element = true, _default_value = 0.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // Initializing stresses ElementTypeMap stress_components; /// TODO this is ugly af, maybe add a function to FEEngine for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { UInt nb_components = 0; // Getting number of components for each element type #define GET_(type) nb_components = ElementClass::getNbStressComponents() AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_); #undef GET_ stress_components(nb_components, type); } stress.initialize( getFEEngine(), _spatial_dimension = _all_dimensions, _element_kind = _ek_structural, _nb_component = [&stress_components](ElementType type, GhostType /*unused*/) -> UInt { return stress_components(type); }); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFEEngineBoundary() { /// TODO: this function should not be reimplemented /// we're just avoiding a call to Model::initFEEngineBoundary() } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper().setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (not dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(velocity, nb_degree_of_freedom, "velocity"); this->allocNodalField(acceleration, nb_degree_of_freedom, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } /* Only allocate the mass if the "lumped" mode is ennabled. * Also it is not a 1D array, but has an element for every * DOF, which are most of the time equal, but makes handling * some operations a bit simpler. */ if (time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { - this->allocNodalField(this->mass, this->nb_degree_of_freedom, - "lumpedmass"); + this->allocateLumpedMassArray(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { element_material.initialize(mesh, _element_kind = _ek_structural, _default_value = 0, _with_nb_element = true); getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); if (not need_to_reassemble_stiffness) { return; } if (not getDOFManager().hasMatrix("K")) { getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().zeroMatrix("K"); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } need_to_reassemble_stiffness = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +bool StructuralMechanicsModel::allocateLumpedMassArray() { + if (this->mass != nullptr) // Already allocated, so nothing to do. + { + return true; + }; + + // now allocate it + this->allocNodalField(this->mass, this->nb_degree_of_freedom, "lumped_mass"); + + return true; +}; + /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs.get(); return mesh.createNodalField(uint_nodal_fields[field_name], group_name); } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } UInt padding_size = 0; if (padding_flag) { padding_size = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation.get(), group_name, n, 0, padding_size); } if (field_name == "velocity") { return mesh.createStridedNodalField(velocity.get(), group_name, n, 0, padding_size); } if (field_name == "acceleration") { return mesh.createStridedNodalField(acceleration.get(), group_name, n, 0, padding_size); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "force") { return mesh.createStridedNodalField(external_force.get(), group_name, n, 0, padding_size); } if (field_name == "external_force") { return mesh.createStridedNodalField(external_force.get(), group_name, n, 0, padding_size); } if (field_name == "momentum") { return mesh.createStridedNodalField(external_force.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force.get(), group_name, n, 0, padding_size); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField(internal_force.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "mass") { + AKANTU_DEBUG_ASSERT(this->mass.get() != nullptr, + "The lumped mass matrix was not allocated."); return mesh.createStridedNodalField(this->mass.get(), group_name, n, 0, padding_size); } return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*unused*/, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; if (field_name == "element_index_by_material") { field = mesh.createElementalField( field_name, group_name, spatial_dimension, kind); } if (field_name == "stress") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(stress); field = mesh.createElementalField( stress, group_name, this->spatial_dimension, kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) const { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") { assembleStiffnessMatrix(); } else if (id == "M") { assembleMassMatrix(); } } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & id) { - if (id == "M") { - assembleLumpedMassMatrix(); + if ("M" == id) { + this->assembleLumpedMassMatrix(); } + return; } /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); auto & dof_manager = getDOFManager(); assembleInternalForce(); + // Ensures that the matrix are assembled. + if (dof_manager.hasMatrix("K")) { + this->assembleMatrix("K"); + } + if (dof_manager.hasMatrix("M")) { + this->assembleMatrix("M"); + } + if (dof_manager.hasLumpedMatrix("M")) { + this->assembleLumpedMassMatrix(); + } + + /* This is essentially a summing up of forces + * first the external forces are counted for and then stored inside the + * residual. + */ dof_manager.assembleToResidual("displacement", *external_force, 1); dof_manager.assembleToResidual("displacement", *internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); + auto & dof_manager = this->getDOFManager(); + if ("external" == residual_part) { - this->getDOFManager().assembleToResidual("displacement", + dof_manager.assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForce(); - this->getDOFManager().assembleToResidual("displacement", + dof_manager.assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes std::tuple StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } + case _explicit_lumped_mass: { // Taken from the solid mechanics part + return std::make_tuple("explicit_lumped", + TimeStepSolverType::_dynamic_lumped); + } + case _explicit_consistent_mass: { // Taken from the solid mechanics part + return std::make_tuple("explicit", TimeStepSolverType::_dynamic); + } default: + std::cout << "UNKOWN." << std::endl; return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { + case TimeStepSolverType::_dynamic_lumped: { // Taken from the solid mechanic + // part + options.non_linear_solver_type = NonLinearSolverType::_lumped; + options.integration_scheme_type["displacement"] = + IntegrationSchemeType::_central_difference; + options.solution_type["displacement"] = IntegrationScheme::_acceleration; + break; + } case TimeStepSolverType::_static: { - options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; + options.non_linear_solver_type = + NonLinearSolverType::_newton_raphson; // _linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } - case TimeStepSolverType::_dynamic: { +#if 1 + case TimeStepSolverType::_dynamic: { // Copied from solid + if (this->method == _explicit_consistent_mass) { + options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; + options.integration_scheme_type["displacement"] = + IntegrationSchemeType::_central_difference; + options.solution_type["displacement"] = IntegrationScheme::_acceleration; + } else { + options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; + options.integration_scheme_type["displacement"] = + IntegrationSchemeType::_trapezoidal_rule_2; + options.solution_type["displacement"] = IntegrationScheme::_displacement; + } + break; + } +#else + case TimeStepSolverType::_dynamic: { // Original options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; break; } +#endif default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce() { + internal_force->zero(); computeStresses(); for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { assembleInternalForce(type, _not_ghost); // assembleInternalForce(type, _ghost); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce(ElementType type, GhostType gt) { auto & fem = getFEEngine(); auto & sigma = stress(type, gt); auto ndof = getNbDegreeOfFreedom(type); auto nb_nodes = mesh.getNbNodesPerElement(type); auto ndof_per_elem = ndof * nb_nodes; Array BtSigma(fem.getNbIntegrationPoints(type) * mesh.getNbElement(type), ndof_per_elem, "BtSigma"); fem.computeBtD(sigma, BtSigma, type, gt); Array intBtSigma(0, ndof_per_elem, "intBtSigma"); fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt); getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force, type, gt, -1.); } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getKineticEnergy() { - if (not this->getDOFManager().hasMatrix("M")) { - return 0.; - } - Real ekin = 0.; - UInt nb_nodes = mesh.getNbNodes(); + const UInt nb_nodes = mesh.getNbNodes(); + const UInt nb_degree_of_freedom = this->nb_degree_of_freedom; + Real ekin = 0.; // used to sum up energy (is divided by two at the very end) - Array Mv(nb_nodes, nb_degree_of_freedom); - this->getDOFManager().assembleMatMulVectToArray("displacement", "M", - *this->velocity, Mv); + //if mass matrix was not assembled, assemble it now + this->assembleMassMatrix(); - for (auto && data : zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom), - make_view(*this->velocity, nb_degree_of_freedom))) { - ekin += std::get<2>(data).dot(std::get<1>(data)) * - static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); + if (this->getDOFManager().hasLumpedMatrix("M")) { + /* This code computes the kinetic energy for the case when the mass is + * lumped. It is based on the solid mechanic equivalent. + */ + AKANTU_DEBUG_ASSERT(this->mass != nullptr, + "The lumped mass is not allocated."); + + if (this->need_to_reassemble_lumpedMass) { + this->assembleLumpedMatrix("M"); + } + + /* Iterating over all nodes. + * Important the velocity and mass also contains the rotational parts. + * However, they can be handled in an uniform way. */ + for (auto && data : + zip(arange(nb_nodes), make_view(*this->velocity, nb_degree_of_freedom), + make_view(*this->mass, nb_degree_of_freedom))) { + const UInt n = std::get<0>(data); // This is the ID of the current node + + if (not mesh.isLocalOrMasterNode( + n)) // Only handle the node if it belongs to us. + { + continue; + } + + const auto & v = + std::get<1>(data); // Get the velocity and mass of that node. + const auto & m = std::get<2>(data); + Real mv2 = 0.; // Contribution of this node. + + for (UInt i = 0; i < nb_degree_of_freedom; ++i) { + /* In the solid mechanics part, only masses that are above a certain + * value are considered. + * However, the structural part, does not do this. */ + const Real v_ = v(i); + const Real m_ = m(i); + mv2 += v_ * v_ * m_; + } // end for(i): going through the components + + ekin += mv2; // add continution + } // end for(n): iterating through all nodes + } else if (this->getDOFManager().hasMatrix("M")) { + /* Handle the case where no lumped mass is there. + * This is basically the original code. + */ + if (this->need_to_reassemble_mass) { + this->assembleMassMatrix(); + } + + Array Mv(nb_nodes, nb_degree_of_freedom); + this->getDOFManager().assembleMatMulVectToArray("displacement", "M", + *this->velocity, Mv); + + for (auto && data : + zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom), + make_view(*this->velocity, nb_degree_of_freedom))) { + if (mesh.isLocalOrMasterNode(std::get<0>( + data))) // only consider the node if we are blonging to it + { + ekin += std::get<2>(data).dot(std::get<1>(data)); + } + } + } else { + /* This is the case where no mass is present, for whatever reason, such as + * the static case. We handle it specially be returning directly zero. + * However, by doing that there will not be a syncronizing event as in the + * other cases. Which is faster, but could be a problem in case the user + * expects this. + * + * Another not is, that the solid mechanics part, would generate an error in + * this clause. But, since the original implementation of the structural + * part, did not do that, I, Philip, decided to refrain from that. However, + * it is an option that should be considered. + */ + return 0.; } + // Sum up across the comunicator mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); - return ekin / 2.; + return ekin / 2.; // finally divide the energy by two } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getPotentialEnergy() { Real epot = 0.; UInt nb_nodes = mesh.getNbNodes(); + //if stiffness matrix is not assembled, do it + this->assembleStiffnessMatrix(); + Array Ku(nb_nodes, nb_degree_of_freedom); this->getDOFManager().assembleMatMulVectToArray( "displacement", "K", *this->displacement_rotation, Ku); for (auto && data : zip(arange(nb_nodes), make_view(Ku, nb_degree_of_freedom), make_view(*this->displacement_rotation, nb_degree_of_freedom))) { epot += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } mesh.getCommunicator().allReduce(epot, SynchronizerOperation::_sum); return epot / 2.; } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getEnergy(const ID & energy) { if (energy == "kinetic") { return getKineticEnergy(); } if (energy == "potential") { return getPotentialEnergy(); } return 0; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeForcesByLocalTractionArray( const Array & tractions, ElementType type) { AKANTU_DEBUG_IN(); auto nb_element = getFEEngine().getMesh().getNbElement(type); auto nb_nodes_per_element = getFEEngine().getMesh().getNbNodesPerElement(type); auto nb_quad = getFEEngine().getNbIntegrationPoints(type); // check dimension match AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(), "element type dimension does not match the dimension of boundaries : " << getFEEngine().getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT( tractions.size() == nb_quad * nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom, "the number of components should be the spatial " "dimension of the problem"); Array Ntbs(nb_element * nb_quad, nb_degree_of_freedom * nb_nodes_per_element); auto & fem = getFEEngine(); fem.computeNtb(tractions, Ntbs, type); // allocate the vector that will contain the integrated values auto name = id + std::to_string(type) + ":integral_boundary"; Array int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element, name); // do the integration getFEEngine().integrate(Ntbs, int_funct, nb_degree_of_freedom * nb_nodes_per_element, type); // assemble the result into force vector getDOFManager().assembleElementalArrayLocalArray(int_funct, *external_force, type, _not_ghost, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeForcesByGlobalTractionArray( const Array & traction_global, ElementType type) { AKANTU_DEBUG_IN(); UInt nb_element = mesh.getNbElement(type); UInt nb_quad = getFEEngine().getNbIntegrationPoints(type); Array traction_local(nb_element * nb_quad, nb_degree_of_freedom, id + ":structuralmechanics:imposed_linear_load"); auto R_it = getFEEngineClass() .getShapeFunctions() .getRotations(type) .begin(nb_degree_of_freedom, nb_degree_of_freedom); auto Te_it = traction_global.begin(nb_degree_of_freedom); auto te_it = traction_local.begin(nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it) { for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) { // turn the traction in the local referential te_it->template mul(*R_it, *Te_it); } } computeForcesByLocalTractionArray(traction_local, type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::afterSolveStep(bool converged) { if (converged) { assembleInternalForce(); } } } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index 4b5b6fd0e..3c80f2bd2 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,404 +1,398 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Philip Mueller * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Apr 01 2021 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_named_argument.hh" #include "boundary_condition.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ #define AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeStructural; } // namespace akantu namespace akantu { struct StructuralMaterial { Real E{0}; Real A{1}; Real I{0}; Real Iz{0}; Real Iy{0}; Real GJ{0}; Real rho{0}; Real t{0}; Real nu{0}; }; class StructuralMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineType = FEEngineTemplate; StructuralMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "structural_mechanics_model"); ~StructuralMechanicsModel() override; /// Init full model void initFullImpl(const ModelOptions & options) override; /// Init boundary FEEngine void initFEEngineBoundary() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from SolverCallback */ /* ------------------------------------------------------------------------ */ /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) const override; /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback to assemble the residual (rhs) void assembleResidual() override; void assembleResidual(const ID & residual_part) override; bool canSplitResidual() const override { return true; } void afterSolveStep(bool converged) override; /// compute kinetic energy Real getKineticEnergy(); /// compute potential energy Real getPotentialEnergy(); /// compute the specified energy Real getEnergy(const ID & energy); - /** * \brief This function computes the an approximation of the lumped mass. * * The mass is computed by looping over all beams and computing their mass. * The mass of a single beam is computed by the (initial) length of the beam, * its cross sectional area and its density. * The beam mass is then equaly distributed among the two nodes. * * For computing the rotational inertia, the function assumes that the mass of * a node is uniformaly distributed inside a disc (2D) or a sphere (3D). The * size of that disc, depends on the volume of the beam. * * Note that the computation of the mass is not unambigius. * The reason for this is, that the units of `StructralMaterial::rho` are not * clear. By default the function assumes that its unit are 'Mass per Volume'. * However, this makes the computed mass different than the consistent mass, * which seams to assume that its units are 'mass per unit length'. * The main difference between thge two are not the values, but that the * first version depends on `StructuralMaterial::A` while the later does not. * By defining the macro `AKANTU_STRUCTURAL_MECHANICS_CONSISTENT_LUMPED_MASS` * the function will compute the mass in a way that is consistent with the * consistent mass matrix. * * \note The lumped mass is not stored inside the DOFManager. * * \param ghost_type Should ghost types be computed. */ void assembleLumpedMassMatrix(); /* ------------------------------------------------------------------------ */ /* Virtual methods from Model */ /* ------------------------------------------------------------------------ */ protected: /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; static UInt getNbDegreeOfFreedom(ElementType type); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /// initialize the model void initModel() override; /// compute the stresses per elements void computeStresses(); /// compute the nodal forces void assembleInternalForce(); /// compute the nodal forces for an element type void assembleInternalForce(ElementType type, GhostType gt); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMassMatrix(); protected: /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMassMatrix(GhostType ghost_type); /// computes rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template void assembleStiffnessMatrix(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the StructuralMechanicsModel::external_force vector AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /** - * Returns a non const reference to the array that stores the lumped - * mass. + * Returns a const reference to the array that stores the lumped mass. * * The returned array has dimension `N x d` where `N` is the number of nodes * and `d`, is the number of degrees of freedom per node. */ - inline Array & getLumpedMass() { - if (this->mass == nullptr) { - AKANTU_EXCEPTION("The pointer to the mass was not allocated."); - }; - return *(this->mass); - }; - - /** - * Returns a constant reference to the array that stores the lumped - * mass. - */ inline const Array & getLumpedMass() const { if (this->mass == nullptr) { AKANTU_EXCEPTION("The pointer to the mass was not allocated."); }; return *(this->mass); }; + // These function is an alias, for compability with the solid mechanics + inline const Array & getMass() const { return this->getLumpedMass(); } + + // Creates the array for storing the mass + bool allocateLumpedMassArray(); + /** * Tests if *this has a lumped mass pointer. */ inline bool hasLumpedMass() const { return (this->mass != nullptr); }; AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt); /** * \brief This function adds the `StructuralMaterial` material to the list of * materials managed by *this. * * It is important that this function might invalidate all references to * structural materials, that were previously optained by `getMaterial()`. * * \param material The new material. * * \return The ID of the material that was added. * * \note The return type is is new. */ UInt addMaterial(StructuralMaterial & material, const ID & name = ""); const StructuralMaterial & getMaterialByElement(const Element & element) const; /** * \brief Returns the ith material of *this. * \param i The ith material */ const StructuralMaterial & getMaterial(UInt material_index) const; const StructuralMaterial & getMaterial(const ID & name) const; /** * \brief Returns the number of the different materials inside *this. */ UInt getNbMaterials() const { return materials.size(); } /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis void computeForcesByGlobalTractionArray(const Array & traction_global, ElementType type); /// Compute Linear load function set in local axis void computeForcesByLocalTractionArray(const Array & tractions, ElementType type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array std::unique_ptr> displacement_rotation; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// forces array std::unique_ptr> internal_force; /// forces array std::unique_ptr> external_force; /** * \brief This is the "lumped" mass array. * * It is a bit special, since it is not a one dimensional array, bit it is * actually a matrix. The number of rows equals the number of nodes. The * number of colums equals the number of degrees of freedoms per node. This * layout makes the thing a bit more simple. * * Note that it is only allocated in case, the "Lumped" mode is enabled. */ std::unique_ptr> mass; /// boundaries array std::unique_ptr> blocked_dofs; /// stress array ElementTypeMapArray stress; ElementTypeMapArray element_material; // Define sets of beams ElementTypeMapArray set_ID; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray rotation_matrix; // /// analysis method check the list in akantu::AnalysisMethod // AnalysisMethod method; /// flag defining if the increment must be computed or not bool increment_flag; bool need_to_reassemble_mass{true}; bool need_to_reassemble_stiffness{true}; + bool need_to_reassemble_lumpedMass{true}; /* ------------------------------------------------------------------------ */ std::vector materials; std::map materials_names_to_id; }; } // namespace akantu #include "structural_mechanics_model_inline_impl.hh" #endif /* AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ */ diff --git a/src/model/structural_mechanics/structural_mechanics_model_mass.cc b/src/model/structural_mechanics/structural_mechanics_model_mass.cc index 7c55c5606..9aad97666 100644 --- a/src/model/structural_mechanics/structural_mechanics_model_mass.cc +++ b/src/model/structural_mechanics/structural_mechanics_model_mass.cc @@ -1,213 +1,219 @@ /** * @file structural_mechanics_model_mass.cc * * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * * @date creation: Mon Jul 07 2014 * @date last modification: Thu Mar 04 2021 * * @brief function handling mass computation * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "integrator_gauss.hh" #include "material.hh" #include "shape_structural.hh" #include "structural_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class ComputeRhoFunctorStruct { public: explicit ComputeRhoFunctorStruct(const StructuralMechanicsModel & model) : model(model){}; void operator()(Matrix & rho, const Element & element) const { const auto & mat = model.getMaterialByElement(element); const Real mat_rho = mat.rho; // This is the density const Real mat_A = mat.A; // This is the cross section const Real mat_mass_per_length = mat_rho * mat_A; // Mass of the beam per unit length rho.set( mat_mass_per_length); // The integrator _recquiers_ mass per unit length } private: const StructuralMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleMassMatrix() { AKANTU_DEBUG_IN(); if (not need_to_reassemble_mass) { return; } if (not getDOFManager().hasMatrix("M")) { getDOFManager().getNewMatrix("M", getMatrixType("M")); } this->getDOFManager().zeroMatrix("M"); assembleMassMatrix(_not_ghost); need_to_reassemble_mass = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleMassMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); ComputeRhoFunctorStruct compute_rho(*this); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_structural)) { fem.assembleFieldMatrix(compute_rho, "M", "displacement", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleLumpedMassMatrix() { + + if (not this->need_to_reassemble_lumpedMass) + return; + allocNodalField(this->mass, nb_degree_of_freedom, "lumped_mass"); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().zeroLumpedMatrix("M"); // Preparing const UInt nb_nodes = mesh.getNbNodes(); const auto & nodes = mesh.getNodes(); const auto & node_it = make_view(nodes, spatial_dimension).begin(); auto & lumped_mass = *(this->mass); lumped_mass.zero(); // Set the lumped mass to zero Array volumes(nb_nodes, 1, 0., "volumes"); /// We now compute the mass and the volume, but not the inertia for (auto type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { const auto & element_material_id = this->element_material(type); const auto & connectivity = mesh.getConnectivity(type); const UInt nb_element = connectivity.size(); if (type != _bernoulli_beam_3 or type != _bernoulli_beam_2) { AKANTU_EXCEPTION( "The lumped mass was not implemented for non Bernoulli Beams"); } // Now iterate through all the connectivity for (auto && data : zip(make_view(connectivity, 2), element_material_id)) { const auto & conn = std::get<0>(data); auto node1 = conn(0); auto node2 = conn(0); const auto & material = this->materials.at(std::get<1>(data)); const Real rho = material.rho; const Real cross_section = material.A; const Vector coord_node_1 = node_it[node1]; const Vector coord_node_2 = node_it[node2]; const Real length = coord_node_1.distance(coord_node_2); const Real volume = length * cross_section; const Real mass = volume * rho; // Now distribute the mass at the right places for (UInt d = 0; d != spatial_dimension; ++d) { lumped_mass(node1, d) += mass / 2.; lumped_mass(node2, d) += mass / 2.; }; // end for(d): volumes(node1) += volume / 2.; // this entry never point to mass volumes(node2) += volume / 2.; } } const Real pi = std::atan(1.) * 4; // We now compute the inertia. if (spatial_dimension == 2) { /* This is the 2D case, so we are assuming that the mass is inside a disc. * Which is given as * \begin{align} * I := m \cdot r^2 * \end{align} * * The radius is obtained by assuming that the volume, that is associated to * a beam, forms a uniformly disc. From this volume we can compute the * radius. */ for (auto && data : zip(make_view(lumped_mass, nb_degree_of_freedom), volumes)) { const Real volume = std::get<1>(data); auto & masses = std::get<0>(data); const Real r2 = volume / pi; // The square of the volume const Real inertia = masses(0) * r2; masses(spatial_dimension) = inertia; } } else if (spatial_dimension == 3) { /* This is essentially the same as for 2D only that we assume here, * that the mass is uniformly distributed inside a sphere. * And thus we have * \begin{align} * I := \frac{2}{5} m \cdot r^2 * \end{align} * * We also have to set three values, for the three axis. * But since we assume a sphere, they are all the same. */ for (auto && data : zip(make_view(lumped_mass, nb_degree_of_freedom), volumes)) { const Real volume = std::get<1>(data); auto & masses = std::get<0>(data); const Real r2 = std::pow((volume * 3.) / (4 * pi), 2. / 3.); const Real inertia = (2. / 5.) * masses(0) * r2; for (UInt d = spatial_dimension; d < masses.size(); ++d) { masses(d) = inertia; } } } else { AKANTU_EXCEPTION("The dimension " << spatial_dimension << " is not known."); } this->getDOFManager().assembleToLumpedMatrix("displacement", lumped_mass, "M"); + this->need_to_reassemble_lumpedMass = false; + return; } } // namespace akantu diff --git a/src/solver/solver_vector.hh b/src/solver/solver_vector.hh index 8592a6320..d05727fa3 100644 --- a/src/solver/solver_vector.hh +++ b/src/solver/solver_vector.hh @@ -1,98 +1,96 @@ /** * @file solver_vector.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue May 26 2020 * * @brief Solver vector interface base class * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_HH_ #define AKANTU_SOLVER_VECTOR_HH_ namespace akantu { class DOFManager; } namespace akantu { class SolverVector { public: SolverVector(DOFManager & dof_manager, const ID & id = "solver_vector") : id(id), _dof_manager(dof_manager) {} SolverVector(const SolverVector & vector, const ID & id = "solver_vector") : id(id), _dof_manager(vector._dof_manager) {} virtual ~SolverVector() = default; // resize the vector to the size of the problem virtual void resize() = 0; // clear the vector virtual void set(Real val) = 0; void zero() { this->set({}); } virtual operator const Array &() const = 0; virtual Int size() const = 0; virtual Int localSize() const = 0; virtual SolverVector & operator+(const SolverVector & y) = 0; virtual SolverVector & operator=(const SolverVector & y) = 0; UInt & release() { return release_; } UInt release() const { return release_; } virtual void printself(std::ostream & stream, int indent = 0) const = 0; - /** - * The default implementation of this function, which should be suitable for - * everyone, converts `*this` into an array and then performs the checks on - * the result. - */ virtual bool isFinite() const = 0; + /// Returns `true` if `*this` is distributed or not. + virtual bool isDistributed() const { return false; } + protected: ID id; /// Underlying dof manager DOFManager & _dof_manager; UInt release_{0}; }; inline std::ostream & operator<<(std::ostream & stream, SolverVector & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_HH_ */ diff --git a/src/solver/solver_vector_default.hh b/src/solver/solver_vector_default.hh index 1a56fdf32..49e11dd24 100644 --- a/src/solver/solver_vector_default.hh +++ b/src/solver/solver_vector_default.hh @@ -1,146 +1,156 @@ /** * @file solver_vector_default.hh * * @author Nicolas Richart * * @date creation: Tue Jan 01 2019 * @date last modification: Sat May 23 2020 * * @brief Solver vector interface to Array * * * @section LICENSE * * Copyright (©) 2018-2021 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 "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_DEFAULT_HH_ #define AKANTU_SOLVER_VECTOR_DEFAULT_HH_ namespace akantu { class DOFManagerDefault; } // namespace akantu namespace akantu { class SolverVectorArray : public SolverVector { public: SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id); SolverVectorArray(const SolverVectorArray & vector, const ID & id); ~SolverVectorArray() override = default; virtual Array & getVector() = 0; virtual const Array & getVector() const = 0; void printself(std::ostream & stream, int indent = 0) const override { std::string space(indent, AKANTU_INDENT); stream << space << "SolverVectorArray [" << std::endl; stream << space << " + id: " << id << std::endl; this->getVector().printself(stream, indent + 1); stream << space << "]" << std::endl; } + + using SolverVector::isDistributed; }; /* -------------------------------------------------------------------------- */ template class SolverVectorArrayTmpl : public SolverVectorArray { public: SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, Array_ & vector, const ID & id = "solver_vector_default") : SolverVectorArray(dof_manager, id), dof_manager(dof_manager), vector(vector) {} template ::value> * = nullptr> SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, const ID & id = "solver_vector_default") : SolverVectorArray(dof_manager, id), dof_manager(dof_manager), vector(0, 1, id + ":vector") {} SolverVectorArrayTmpl(const SolverVectorArrayTmpl & vector, const ID & id = "solver_vector_default") : SolverVectorArray(vector, id), dof_manager(vector.dof_manager), vector(vector.vector) {} operator const Array &() const override { return getVector(); }; virtual operator Array &() { return getVector(); }; SolverVector & operator+(const SolverVector & y) override; SolverVector & operator=(const SolverVector & y) override; void resize() override { static_assert(not std::is_const>::value, "Cannot resize a const Array"); this->vector.resize(this->localSize(), 0.); ++this->release_; } void set(Real val) override { static_assert(not std::is_const>::value, "Cannot clear a const Array"); this->vector.set(val); ++this->release_; } + virtual + bool + isDistributed() + const + override + { return false; } + + public: Array & getVector() override { return vector; } const Array & getVector() const override { return vector; } Int size() const override; Int localSize() const override; virtual Array & getGlobalVector() { return this->vector; } virtual void setGlobalVector(const Array & solution) { this->vector.copy(solution); } bool isFinite() const override { return vector.isFinite(); } protected: DOFManagerDefault & dof_manager; Array_ vector; template friend class SolverVectorArrayTmpl; }; /* -------------------------------------------------------------------------- */ using SolverVectorDefault = SolverVectorArrayTmpl>; /* -------------------------------------------------------------------------- */ template using SolverVectorDefaultWrap = SolverVectorArrayTmpl; template decltype(auto) make_solver_vector_default_wrap(DOFManagerDefault & dof_manager, Array & vector) { return SolverVectorDefaultWrap(dof_manager, vector); } } // namespace akantu /* -------------------------------------------------------------------------- */ #include "solver_vector_default_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLVER_VECTOR_DEFAULT_HH_ */ diff --git a/src/solver/solver_vector_distributed.hh b/src/solver/solver_vector_distributed.hh index 951454c93..a928e0ac8 100644 --- a/src/solver/solver_vector_distributed.hh +++ b/src/solver/solver_vector_distributed.hh @@ -1,61 +1,63 @@ /** * @file solver_vector_distributed.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue Jan 01 2019 * * @brief Solver vector interface for distributed arrays * * * @section LICENSE * * Copyright (©) 2018-2021 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 "solver_vector_default.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ #define AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ namespace akantu { class SolverVectorDistributed : public SolverVectorDefault { public: SolverVectorDistributed(DOFManagerDefault & dof_manager, const ID & id = "solver_vector_mumps"); SolverVectorDistributed(const SolverVectorDefault & vector, const ID & id = "solver_vector_mumps"); Array & getGlobalVector() override; void setGlobalVector(const Array & solution) override; bool isFinite() const override; + virtual bool isDistributed() const override { return true; } + protected: // full vector in case it needs to be centralized on master std::unique_ptr> global_vector; }; } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ */ diff --git a/src/solver/solver_vector_petsc.hh b/src/solver/solver_vector_petsc.hh index 223d59a98..570220445 100644 --- a/src/solver/solver_vector_petsc.hh +++ b/src/solver/solver_vector_petsc.hh @@ -1,208 +1,215 @@ /** * @file solver_vector_petsc.hh * * @author Nicolas Richart * * @date creation: Tue Jan 01 2019 * @date last modification: Fri Jul 24 2020 * * @brief Solver vector interface for PETSc * * * @section LICENSE * * Copyright (©) 2018-2021 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 "dof_manager_petsc.hh" #include "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_PETSC_HH_ #define AKANTU_SOLVER_VECTOR_PETSC_HH_ namespace akantu { class DOFManagerPETSc; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ namespace internal { /* ------------------------------------------------------------------------ */ class PETScVector { public: virtual ~PETScVector() = default; operator Vec &() { return x; } operator const Vec &() const { return x; } Int size() const { PetscInt n; PETSc_call(VecGetSize, x, &n); return n; } Int local_size() const { PetscInt n; PETSc_call(VecGetLocalSize, x, &n); return n; } AKANTU_GET_MACRO_NOT_CONST(Vec, x, auto &); AKANTU_GET_MACRO(Vec, x, const auto &); protected: Vec x{nullptr}; }; } // namespace internal /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class SolverVectorPETSc : public SolverVector, public internal::PETScVector { public: SolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); SolverVectorPETSc(const SolverVectorPETSc & vector, const ID & id = "solver_vector_petsc"); SolverVectorPETSc(Vec x, DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); ~SolverVectorPETSc() override; // resize the vector to the size of the problem void resize() override; void set(Real val) override; operator const Array &() const override; SolverVector & operator+(const SolverVector & y) override; SolverVector & operator=(const SolverVector & y) override; SolverVectorPETSc & operator=(const SolverVectorPETSc & y); /// get values using processors global indexes void getValues(const Array & idx, Array & values) const; /// get values using processors local indexes void getValuesLocal(const Array & idx, Array & values) const; /// adding values to the vector using the global indices void addValues(const Array & gidx, const Array & values, Real scale_factor = 1.); /// adding values to the vector using the local indices void addValuesLocal(const Array & lidx, const Array & values, Real scale_factor = 1.); Int size() const override { return internal::PETScVector::size(); } Int localSize() const override { return internal::PETScVector::local_size(); } void printself(std::ostream & stream, int indent = 0) const override; + virtual + bool + isDistributed() + const + override + { return true; } + protected: void applyModifications(); void updateGhost(); protected: DOFManagerPETSc & dof_manager; // used for the conversion operator Array cache; }; /* -------------------------------------------------------------------------- */ namespace internal { /* ------------------------------------------------------------------------ */ template class PETScWrapedVector : public PETScVector { public: PETScWrapedVector(Array && array) : array(array) { PETSc_call(VecCreateSeqWithArray, PETSC_COMM_SELF, 1, array.size(), array.storage(), &x); } ~PETScWrapedVector() override { PETSc_call(VecDestroy, &x); } private: Array array; }; /* ------------------------------------------------------------------------ */ template class PETScLocalVector : public PETScVector { public: PETScLocalVector(const Vec & g) : g(g) { PETSc_call(VecGetLocalVectorRead, g, x); } PETScLocalVector(const SolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {} ~PETScLocalVector() override { PETSc_call(VecRestoreLocalVectorRead, g, x); PETSc_call(VecDestroy, &x); } private: const Vec & g; }; template <> class PETScLocalVector : public PETScVector { public: PETScLocalVector(Vec & g) : g(g) { PETSc_call(VecGetLocalVectorRead, g, x); } PETScLocalVector(SolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {} ~PETScLocalVector() override { PETSc_call(VecRestoreLocalVectorRead, g, x); PETSc_call(VecDestroy, &x); } private: Vec & g; }; /* ------------------------------------------------------------------------ */ template decltype(auto) make_petsc_wraped_vector(Array && array) { return PETScWrapedVector(std::forward(array)); } template < typename V, std::enable_if_t>::value> * = nullptr> decltype(auto) make_petsc_local_vector(V && vec) { constexpr auto read_only = std::is_const>::value; return PETScLocalVector(vec); } template >::value> * = nullptr> decltype(auto) make_petsc_local_vector(V && vec) { constexpr auto read_only = std::is_const>::value; return PETScLocalVector( dynamic_cast &>(vec)); } } // namespace internal } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_PETSC_HH_ */ diff --git a/test/test_model/test_common/CMakeLists.txt b/test/test_model/test_common/CMakeLists.txt index 391ed3b73..93e2ffacc 100644 --- a/test/test_model/test_common/CMakeLists.txt +++ b/test/test_model/test_common/CMakeLists.txt @@ -1,40 +1,45 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Wed Jan 30 2019 # @date last modification: Tue Feb 26 2019 # # @brief CMakeLists for common part of Models # # # @section LICENSE # # Copyright (©) 2018-2021 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 . # #=============================================================================== add_akantu_test(test_model_solver "Test for the solvers") add_akantu_test(test_non_local_toolbox "Test of the functionalities in the non-local toolbox") add_mesh(dof_manager_mesh mesh.geo 3 1) register_gtest_sources( SOURCES test_dof_manager.cc PACKAGE core) register_gtest_test(test_dof_manager DEPENDS dof_manager_mesh) + +register_gtest_sources( + SOURCES test_dof_mesh_distribute.cc + PACKAGE parallel) +register_gtest_test(test_dof_mesh_distribute DEPENDS dof_manager_mesh) diff --git a/test/test_model/test_common/test_dof_mesh_distribute.cc b/test/test_model/test_common/test_dof_mesh_distribute.cc new file mode 100644 index 000000000..32f566e36 --- /dev/null +++ b/test/test_model/test_common/test_dof_mesh_distribute.cc @@ -0,0 +1,77 @@ +/** + * @file test_dof_mesh_distribute.cc + * + * @author Philip Müller + * + * @brief Test if the DOFManager can handle the distribution event of teh mesh. + * + * This test is based on the file `test_dof_synchronizer.cc`. + * + * @section LICENSE + * + * Copyright (©) 2010-2021 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 "communicator.hh" +#include "dof_synchronizer.hh" +#include "element_synchronizer.hh" +#include "mesh_io.hh" +#include "mesh_partition_scotch.hh" + +#include +#include + +/* -------------------------------------------------------------------------- */ +#include "io_helper.hh" + +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +TEST(TestMesh, TestMeshDistribute) { + + const UInt spatial_dimension = 3; + + int argc = 0; + char ** argv = nullptr; + initialize(argc, argv); + + const auto & comm = akantu::Communicator::getStaticCommunicator(); + Int prank = comm.whoAmI(); + + Mesh mesh1(spatial_dimension); + Mesh mesh2(spatial_dimension); + + if (prank == 0) { + mesh1.read("mesh.msh"); + mesh2.read("mesh.msh"); + } + + /* This should fail since the mesh is distributed after the + * DOFManager is created. */ + DOFManagerDefault dof_manager(mesh1, "test_dof_manager_1__failing"); + EXPECT_THROW(mesh1.distribute(), debug::Exception); + + /* This will succeed since the mesh is distributed before the manager is + * created. */ + mesh2.distribute(); + EXPECT_NO_THROW( + DOFManagerDefault dof_manager2(mesh2, "test_dof_manager_2__succeeding")); + + finalize(); +}