diff --git a/python/py_aka_error.cc b/python/py_aka_error.cc index 3a45390a9..ab2520407 100644 --- a/python/py_aka_error.cc +++ b/python/py_aka_error.cc @@ -1,70 +1,95 @@ /** * @file py_aka_error.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue May 07 2019 * @date last modification: Tue Sep 29 2020 * * @brief pybind11 interface to aka_error * * * @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_error.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_error(py::module & mod) { - mod.def("setDebugLevel", &debug::setDebugLevel); - mod.def("getDebugLevel", &debug::getDebugLevel); - mod.def("printBacktrace", - [](bool flag) { debug::debugger.printBacktrace(flag); }); + py::module mod_debug = mod.def_submodule("debug"); + + mod.def("setDebugLevel", [](DebugLevel lvl) { + debug::setDebugLevel(lvl); + PyErr_WarnEx(PyExc_DeprecationWarning, + "setDebugLevel() is deprecated, it has moved in the " + "submodule debug", + 1); + }); + mod.def("getDebugLevel", []() { + PyErr_WarnEx(PyExc_DeprecationWarning, + "getDebugLevel() is deprecated, it has moved in the " + "submodule debug", + 1); + return debug::getDebugLevel(); + }); + + mod.def("printBacktrace", [](bool flag) { + debug::debugger.printBacktrace(flag); + PyErr_WarnEx(PyExc_DeprecationWarning, + "printBacktrace() is deprecated, it has moved in the " + "submodule debug", + 1); + }); + + mod_debug.def("setDebugLevel", &debug::setDebugLevel); + mod_debug.def("getDebugLevel", &debug::getDebugLevel); + mod_debug.def("printBacktrace", + [](bool flag) { debug::debugger.printBacktrace(flag); }); py::enum_(mod, "DebugLevel") .value("dblError", dblError) .value("dblException", dblException) .value("dblCritical", dblCritical) .value("dblMajor", dblMajor) .value("dblWarning", dblWarning) .value("dblInfo", dblInfo) .value("dblTrace", dblTrace) .value("dblAccessory", dblAccessory) .value("dblDebug", dblDebug) .value("dblDump", dblDump) .value("dblTest", dblTest) .export_values(); } } // namespace akantu diff --git a/python/py_group_manager.cc b/python/py_group_manager.cc index 02bd999da..349b757ec 100644 --- a/python/py_group_manager.cc +++ b/python/py_group_manager.cc @@ -1,178 +1,181 @@ /** * @file py_group_manager.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Mon Dec 02 2019 * * @brief pybind11 interface to GroupManager, ElementGroup and NodeGroup * * * @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 /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_group_manager(py::module & mod) { /* ------------------------------------------------------------------------ */ py::class_(mod, "NodeGroup") .def( "getNodes", [](NodeGroup & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("__len__", &NodeGroup::size) .def( "__iter__", [](const NodeGroup & self) { return py::make_iterator(self.begin(), self.end()); }, py::keep_alive<0, 1>()) .def("__contains__", [](const NodeGroup & self, UInt node) { return self.find(node) != UInt(-1); }) .def("getName", &NodeGroup::getName) .def("clear", &NodeGroup::clear) .def("empty", &NodeGroup::empty) .def("append", &NodeGroup::append) .def("add", &NodeGroup::add, py::arg("node"), py::arg("check_for_duplicate") = true) .def("remove", &NodeGroup::add); /* ------------------------------------------------------------------------ */ py::class_(mod, "ElementGroup") .def( "getNodeGroup", [](ElementGroup & self) -> decltype(auto) { return self.getNodeGroup(); }, py::return_value_policy::reference) .def("getName", &ElementGroup::getName) .def( "getElements", [](ElementGroup & self) -> decltype(auto) { return self.getElements(); }, py::return_value_policy::reference) .def( "getNodeGroup", [](ElementGroup & self) -> decltype(auto) { return self.getNodeGroup(); }, py::return_value_policy::reference) .def("__len__", [](const ElementGroup & self) { return self.size(); }) .def("clear", [](ElementGroup & self) { self.clear(); }) .def("empty", &ElementGroup::empty) .def("append", &ElementGroup::append) + .def("optimize", &ElementGroup::optimize) .def( "add", [](ElementGroup & self, const Element & element, bool add_nodes, bool check_for_duplicate) { self.add(element, add_nodes, check_for_duplicate); }, py::arg("element"), py::arg("add_nodes") = false, py::arg("check_for_duplicate") = true) .def("fillFromNodeGroup", &ElementGroup::fillFromNodeGroup) .def("addDimension", &ElementGroup::addDimension); /* ------------------------------------------------------------------------ */ py::class_(mod, "GroupManager") .def( "getElementGroup", [](GroupManager & self, const std::string & name) -> decltype(auto) { return self.getElementGroup(name); }, py::return_value_policy::reference) .def("iterateElementGroups", [](GroupManager & self) -> decltype(auto) { std::vector> groups; for (auto & group : self.iterateElementGroups()) { groups.emplace_back(group); } return groups; }) .def("iterateNodeGroups", [](GroupManager & self) -> decltype(auto) { std::vector> groups; for (auto & group : self.iterateNodeGroups()) { groups.emplace_back(group); } return groups; }) .def("createNodeGroup", &GroupManager::createNodeGroup, py::return_value_policy::reference) .def( "createElementGroup", [](GroupManager & self, const std::string & id, - UInt spatial_dimension, bool b) -> decltype(auto) { - return self.createElementGroup(id, spatial_dimension, b); + UInt spatial_dimension, bool replace_group) -> decltype(auto) { + return self.createElementGroup(id, spatial_dimension, + replace_group); }, - py::return_value_policy::reference) + py::arg("id"), py::arg("spatial_dimension"), + py::arg("replace_group") = false, py::return_value_policy::reference) .def("createGroupsFromMeshDataUInt", &GroupManager::createGroupsFromMeshData) .def("createElementGroupFromNodeGroup", &GroupManager::createElementGroupFromNodeGroup, py::arg("name"), py::arg("node_group"), py::arg("dimension") = _all_dimensions) .def( "getNodeGroup", [](GroupManager & self, const std::string & name) -> decltype(auto) { return self.getNodeGroup(name); }, py::return_value_policy::reference) .def( "nodeGroups", [](GroupManager & self) { std::vector groups; for (auto & g : self.iterateNodeGroups()) { groups.push_back(&g); } return groups; }, py::return_value_policy::reference) .def( "elementGroups", [](GroupManager & self) { std::vector groups; for (auto & g : self.iterateElementGroups()) { groups.push_back(&g); } return groups; }, py::return_value_policy::reference) .def("createBoundaryGroupFromGeometry", &GroupManager::createBoundaryGroupFromGeometry); } } // namespace akantu diff --git a/python/py_mesh.cc b/python/py_mesh.cc index a73250d8f..4f2972024 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,238 +1,240 @@ /** * @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); + .def("isPeriodicMaster", &Mesh::isPeriodicMaster) + .def("initMeshFacets", &Mesh::initMeshFacets, + py::arg("id") = "mesh_facets", py::return_value_policy::reference); /* ------------------------------------------------------------------------ */ 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_couplers.cc b/python/py_model_couplers.cc index 7cc6f6424..9ab064d05 100644 --- a/python/py_model_couplers.cc +++ b/python/py_model_couplers.cc @@ -1,122 +1,121 @@ /** * @file py_model_couplers.cc * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Thu Jun 20 2019 * @date last modification: Thu Jun 24 2021 * * @brief Model Coupler python binding * * * @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 /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { template auto register_coupler_solid_contact(py::module & mod, const std::string & name) -> py::class_ { return py::class_(mod, name.c_str(), py::multiple_inheritance()) - .def(py::init, - const ModelType>(), + .def(py::init>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "coupler_solid_contact", - py::arg("dof_manager") = nullptr, - py::arg("model_type") = ModelType::_coupler_solid_contact) + py::arg("dof_manager") = nullptr) .def("applyBC", [](CouplerSolidContact_ & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](CouplerSolidContact_ & self, BC::Neumann::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &CouplerSolidContact_::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getContactMechanicsModel", &CouplerSolidContact_::getContactMechanicsModel, py::return_value_policy::reference); } } // namespace /* -------------------------------------------------------------------------- */ void register_model_couplers(py::module & mod) { register_coupler_solid_contact(mod, "CouplerSolidContact") .def( "getSolidMechanicsModel", [](CouplerSolidContact & self) -> decltype(auto) { return self.getSolidMechanicsModel(); }, py::return_value_policy::reference) .def( "initFull", [](CouplerSolidContact & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method") = _explicit_lumped_mass); register_coupler_solid_contact( mod, "CouplerSolidCohesiveContact") .def( "initFull", [](CouplerSolidCohesiveContact & self, const AnalysisMethod & analysis_method, bool is_extrinsic) { self.initFull(_analysis_method = analysis_method, _is_extrinsic = is_extrinsic); }, py::arg("_analysis_method") = _explicit_lumped_mass, py::arg("_is_extrinsic") = false) .def("checkCohesiveStress", [](CouplerSolidCohesiveContact & self) { return self.checkCohesiveStress(); }) .def( "getSolidMechanicsModelCohesive", [](CouplerSolidCohesiveContact & self) -> decltype(auto) { return self.getSolidMechanicsModelCohesive(); }, py::return_value_policy::reference); } } // namespace akantu diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc index 6d4ce5f3e..ba64b18c5 100644 --- a/python/py_solid_mechanics_model_cohesive.cc +++ b/python/py_solid_mechanics_model_cohesive.cc @@ -1,95 +1,96 @@ /** * @file py_solid_mechanics_model_cohesive.cc * * @author Nicolas Richart * * @date creation: Tue Jul 21 2020 * @date last modification: Tue Sep 29 2020 * * @brief pybind11 interface to SolidMechanicsModelCohesive * * * @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 /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) void register_solid_mechanics_model_cohesive(py::module & mod) { py::class_(mod, "CohesiveElementInserter") .def("setLimit", &CohesiveElementInserter::setLimit); py::class_( mod, "SolidMechanicsModelCohesiveOptions") .def(py::init(), py::arg("analysis_method") = _explicit_lumped_mass, py::arg("is_extrinsic") = false); py::class_( mod, "SolidMechanicsModelCohesive") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model") .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method, bool is_extrinsic) { self.initFull(_analysis_method = analysis_method, _is_extrinsic = is_extrinsic); }, - py::arg("_analysis_method"), py::arg("_is_extrinsic") = false) + py::arg("_analysis_method") = _explicit_lumped_mass, + py::arg("_is_extrinsic") = false) .def("checkCohesiveStress", &SolidMechanicsModelCohesive::checkCohesiveStress) .def("getElementInserter", &SolidMechanicsModelCohesive::getElementInserter, py::return_value_policy::reference) .def("updateAutomaticInsertion", &SolidMechanicsModelCohesive::updateAutomaticInsertion); } } // namespace akantu diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index 176754b69..4ec7b92a6 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,374 +1,368 @@ /** * @file aka_error.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Sep 06 2010 * @date last modification: Wed Feb 24 2021 * * @brief handling of errors * * * @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_error.hh" #include "aka_common.hh" #include "aka_config.hh" #include "aka_iterators.hh" #include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include #include #include #if (defined(READLINK_COMMAND) || defined(ADDR2LINE_COMMAND)) && \ (!defined(_WIN32)) #include #include #endif #include #include #include #include #include #include #include #include #include #ifdef AKANTU_USE_MPI #include #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace debug { - // static void printBacktraceAndExit(int) { std::terminate(); } - - // /* ------------------------------------------------------------------------ - // */ void initSignalHandler() { std::signal(SIGSEGV, &printBacktraceAndExit); - // } - /* ------------------------------------------------------------------------ */ std::string demangle(const char * symbol) { int status; std::string result; char * demangled_name; if ((demangled_name = abi::__cxa_demangle(symbol, nullptr, nullptr, &status)) != nullptr) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; } /* ------------------------------------------------------------------------ */ #if (defined(READLINK_COMMAND) || defined(ADDR2LINK_COMMAND)) && \ (!defined(_WIN32)) std::string exec(const std::string & cmd) { FILE * pipe = popen(cmd.c_str(), "r"); if (pipe == nullptr) { return ""; } char buffer[1024]; std::string result; while (feof(pipe) == 0) { if (fgets(buffer, 128, pipe) != nullptr) { result += buffer; } } result = result.substr(0, result.size() - 1); pclose(pipe); return result; } #endif auto getBacktrace() -> std::vector { std::vector backtrace_lines; #if not defined(_WIN32) #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::string me; char buf[1024]; /* The manpage says it won't null terminate. Let's zero the buffer. */ memset(buf, 0, sizeof(buf)); /* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */ if (readlink("/proc/self/exe", buf, sizeof(buf) - 1) != 0) { me = std::string(buf); } std::ifstream inmaps; inmaps.open("/proc/self/maps"); std::map addr_map; std::string line; while (inmaps.good()) { std::getline(inmaps, line); std::stringstream sstr(line); size_t first = line.find('-'); std::stringstream sstra(line.substr(0, first)); size_t addr; sstra >> std::hex >> addr; std::string lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; if (not lib.empty() and (addr_map.find(lib) == addr_map.end())) { addr_map[lib] = addr; } } if (not me.empty()) { addr_map[me] = 0; } #endif /// \todo for windows this part could be coded using CaptureStackBackTrace /// and SymFromAddr const size_t max_depth = 100; size_t stack_depth; void * stack_addrs[max_depth]; char ** stack_strings; size_t i; stack_depth = backtrace(stack_addrs, max_depth); stack_strings = backtrace_symbols(stack_addrs, stack_depth); /// -1 to remove the call to the printBacktrace function for (i = 1; i < stack_depth; i++) { std::string bt_line(stack_strings[i]); size_t first; size_t second; if ((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) { std::string location = bt_line.substr(0, first); #if defined(READLINK_COMMAND) std::string location_cmd = std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + std::string(" -f ") + location; location = exec(location_cmd); #endif std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str()); size_t f = bt_line.find('['); size_t s = bt_line.find(']'); std::string address = bt_line.substr(f + 1, s - f - 1); std::stringstream sstra(address); size_t addr; sstra >> std::hex >> addr; std::string trace = location + " [" + call + "]"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) auto it = addr_map.find(location); if (it != addr_map.end()) { std::stringstream syscom; syscom << BOOST_PP_STRINGIZE(ADDR2LINE_COMMAND) << " 0x" << std::hex << (addr - it->second) << " -i -e " << location; std::string line = exec(syscom.str()); trace += " (" + line + ")"; } else { #endif std::stringstream sstr_addr; sstr_addr << std::hex << addr; trace += " (0x" + sstr_addr.str() + ")"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) } #endif backtrace_lines.push_back(trace); } else { backtrace_lines.push_back(bt_line); } } free(stack_strings); #endif return backtrace_lines; } /* ------------------------------------------------------------------------ */ void printBacktrace(const std::vector & backtrace) { auto w = size_t(std::floor(std::log10(double(backtrace.size()))) + 1); std::cerr << "BACKTRACE : " << backtrace.size() << " stack frames.\n"; for (auto && data : enumerate(backtrace)) { std::cerr << " [" << std::setw(w) << (std::get<0>(data) + 1) << "] " << std::get<1>(data) << "\n"; } std::cerr << "END BACKTRACE" << std::endl; } /* ------------------------------------------------------------------------ */ namespace { void terminate_handler() { auto eptr = std::current_exception(); auto * t = abi::__cxa_current_exception_type(); auto name = (t != nullptr) ? demangle(t->name()) : std::string("unknown"); try { if (eptr) { std::rethrow_exception(eptr); } else { printBacktrace(); std::cerr << AKANTU_LOCATION << "!! Execution terminated for unknown reasons !!" << std::endl; } } catch (Exception & e) { printBacktrace(e.backtrace()); std::cerr << "!! Uncaught akantu::Exception of type " << name << " !!\nwhat(): \"" << e.what() << "\"" << std::endl; } catch (std::exception & e) { std::cerr << "!! Uncaught exception of type " << name << " !!\nwhat(): \"" << e.what() << "\"" << std::endl; } catch (...) { std::cerr << "!! Something strange of type \"" << name << "\" was thrown.... !!" << std::endl; } if (debugger.printBacktrace()) { std::cerr << "Random generator seed: " << RandomGenerator::seed() << std::endl; printBacktrace(); } } } // namespace /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() noexcept { cout = &std::cerr; level = dblWarning; parallel_context = ""; file_open = false; print_backtrace = false; // initSignalHandler(); std::set_terminate(terminate_handler); } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if (file_open) { dynamic_cast(cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { if (status != 0) { std::terminate(); } std::exit(0); } /*------------------------------------------------------------------------- */ void Debugger::throwException(const std::string & info, const std::string & file, unsigned int line, __attribute__((unused)) bool silent, __attribute__((unused)) const std::string & location, const std::string & module) const noexcept(false) { #if !defined(AKANTU_NDEBUG) if (not silent) { printMessage("###", dblWarning, info + " " + location, module); } #endif debug::Exception ex(info, file, line); ex.setModule(module); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info, const std::string & module) const { if (testLevel(level, module)) { double timestamp = std::chrono::duration_cast>( std::chrono::system_clock::now().time_since_epoch()) .count(); *(cout) << parallel_context << "{" << (size_t)timestamp << "} " << prefix << " " << info << std::endl; } } /* ------------------------------------------------------------------------ */ void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ const DebugLevel & Debugger::getDebugLevel() const { return this->level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if (file_open) { dynamic_cast(cout)->close(); delete cout; } auto * fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } std::ostream & Debugger::getOutputStream() { return *cout; } /* ------------------------------------------------------------------------ */ void Debugger::setParallelContext(int rank, int size) { std::stringstream sstr; UInt pad = std::ceil(std::log10(size)); sstr << "<" << getpid() << ">[R" << std::setfill(' ') << std::right << std::setw(pad) << rank << "|S" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } const DebugLevel & getDebugLevel() { return debugger.getDebugLevel(); } /* ------------------------------------------------------------------------ */ void exit(int status) { Debugger::exit(status); } } // namespace debug } // namespace akantu diff --git a/src/common/aka_error.hh b/src/common/aka_error.hh index fd65dfacd..89ecef3c8 100644 --- a/src/common/aka_error.hh +++ b/src/common/aka_error.hh @@ -1,422 +1,421 @@ /** * @file aka_error.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Tue Feb 09 2021 * * @brief error management and internal exceptions * * * @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 #include #include #include #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ERROR_HH_ #define AKANTU_ERROR_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ enum DebugLevel { dbl0 = 0, dblError = 0, dblAssert = 0, dbl1 = 1, dblException = 1, dblCritical = 1, dbl2 = 2, dblMajor = 2, dbl3 = 3, dblCall = 3, dblSecondary = 3, dblHead = 3, dbl4 = 4, dblWarning = 4, dbl5 = 5, dblInfo = 5, dbl6 = 6, dblIn = 6, dblOut = 6, dbl7 = 7, dbl8 = 8, dblTrace = 8, dbl9 = 9, dblAccessory = 9, dbl10 = 10, dblDebug = 42, dbl100 = 100, dblDump = 100, dblTest = 1337 }; /* -------------------------------------------------------------------------- */ #define AKANTU_LOCATION \ "(" << std::string(__func__) << "(): " << std::string(__FILE__) << ":" \ << std::to_string(__LINE__) \ << ")" // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay) /* -------------------------------------------------------------------------- */ namespace debug { void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel(); void initSignalHandler(); std::string demangle(const char * symbol); template std::string demangle() { return demangle(typeid(T).name()); } template std::string demangle(const T & t) { return demangle(typeid(t).name()); } auto exec(const std::string & cmd) -> std::string; auto getBacktrace() -> std::vector; void printBacktrace(const std::vector & backtrace = getBacktrace()); void exit(int status) __attribute__((noreturn)); /* ------------------------------------------------------------------------ */ /// exception class that can be thrown by akantu class Exception : public std::exception { /* ---------------------------------------------------------------------- */ /* Constructors/Destructors */ /* ---------------------------------------------------------------------- */ protected: explicit Exception(const std::string & info = "") : _info(info) {} public: //! full constructor Exception(const std::string & info, const std::string & file, unsigned int line) : _info(info), _file(file), _line(line) {} /* ---------------------------------------------------------------------- */ /* Methods */ /* ---------------------------------------------------------------------- */ public: const char * what() const noexcept override { return _info.c_str(); } virtual std::string info() const noexcept { std::stringstream stream; stream << debug::demangle(typeid(*this).name()) << " : " << _info << " [" << _file << ":" << _line << "]"; return stream.str(); } public: void setInfo(const std::string & info) { _info = info; } void setFile(const std::string & file) { _file = file; } void setLine(unsigned int line) { _line = line; } void setModule(const std::string & module) { _module = module; } void setBacktrace(const std::vector & backtrace) { backtrace_ = backtrace; } decltype(auto) backtrace() const { return backtrace_; } /* ---------------------------------------------------------------------- */ /* Class Members */ /* ---------------------------------------------------------------------- */ protected: /// exception description and additionals std::string _info; private: /// file it is thrown from std::string _file; /// line it is thrown from unsigned int _line{0}; /// module in which exception was raised std::string _module{"core"}; std::vector backtrace_; }; class CriticalError : public Exception {}; class AssertException : public Exception {}; class NotImplementedException : public Exception {}; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Exception & _this) { stream << _this.what(); return stream; } - /* -------------------------------------------------------------------------- - */ + /* ------------------------------------------------------------------------ */ class Debugger { public: Debugger() noexcept; virtual ~Debugger(); Debugger(const Debugger &) = default; Debugger & operator=(const Debugger &) = default; Debugger(Debugger &&) noexcept = default; Debugger & operator=(Debugger &&) noexcept = default; static void exit(int status) __attribute__((noreturn)); void throwException(const std::string & info, const std::string & file, unsigned int line, bool /*silent*/, const std::string & /*location*/, const std::string & module) const noexcept(false) __attribute__((noreturn)); /*----------------------------------------------------------------------- */ template void throwCustomException(Except ex, const std::string & info, const std::string & file, unsigned int line, const std::string & module) const noexcept(false) __attribute__((noreturn)); /*----------------------------------------------------------------------- */ template void throwCustomException(Except ex, const std::string & file, unsigned int line, const std::string & module_) const noexcept(false) __attribute__((noreturn)); void printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info, const std::string & module_) const; void setOutStream(std::ostream & out) { cout = &out; } std::ostream & getOutStream() { return *cout; } public: void setParallelContext(int rank, int size); void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel() const; void setLogFile(const std::string & filename); std::ostream & getOutputStream(); inline bool testLevel(const DebugLevel & level, const std::string & module = "core") const { auto level_reached = (this->level >= (level)); auto correct_module = (level <= dblCritical) or (modules_to_debug.empty()) or (modules_to_debug.find(module) != modules_to_debug.end()); return level_reached and correct_module; } void printBacktrace(bool on_off) { this->print_backtrace = on_off; } bool printBacktrace() const { return this->print_backtrace; } void addModuleToDebug(const std::string & id) { this->modules_to_debug.insert(id); } void removeModuleToDebug(const std::string & id) { auto it = this->modules_to_debug.find(id); if (it != this->modules_to_debug.end()) { this->modules_to_debug.erase(it); } } void listModules() { for (const auto & module_ : modules_to_debug) { (*cout) << module_ << std::endl; } } private: std::string parallel_context; std::ostream * cout; bool file_open; DebugLevel level; bool print_backtrace; std::set modules_to_debug; }; extern Debugger debugger; // NOLINT } // namespace debug /* -------------------------------------------------------------------------- */ #define AKANTU_STRINGIZE_(str) #str #define AKANTU_STRINGIZE(str) AKANTU_STRINGIZE_(str) /* -------------------------------------------------------------------------- */ #define AKANTU_DEBUG_MODULE AKANTU_STRINGIZE(AKANTU_MODULE) /* -------------------------------------------------------------------------- */ #define AKANTU_STRINGSTREAM_IN(_str, _sstr) \ ; \ do { \ std::stringstream _dbg_s_info; \ _dbg_s_info << _sstr; /* NOLINT */ \ (_str) = _dbg_s_info.str(); \ } while (false) /* -------------------------------------------------------------------------- */ #define AKANTU_EXCEPTION(info) AKANTU_EXCEPTION_(info, false) #define AKANTU_SILENT_EXCEPTION(info) AKANTU_EXCEPTION_(info, true) #define AKANTU_EXCEPTION_(info, silent) \ do { \ std::stringstream _dbg_str; \ _dbg_str << info; /* NOLINT */ \ std::stringstream _dbg_loc; \ _dbg_loc << AKANTU_LOCATION; \ ::akantu::debug::debugger.throwException(_dbg_str.str(), __FILE__, \ __LINE__, silent, _dbg_loc.str(), \ AKANTU_DEBUG_MODULE); \ } while (false) #define AKANTU_CUSTOM_EXCEPTION_INFO(ex, info) \ do { \ std::stringstream _dbg_str; \ _dbg_str << info; /* NOLINT */ \ ::akantu::debug::debugger.throwCustomException( \ ex, _dbg_str.str(), __FILE__, __LINE__, AKANTU_DEBUG_MODULE); \ } while (false) #define AKANTU_CUSTOM_EXCEPTION(ex) \ do { \ ::akantu::debug::debugger.throwCustomException(ex, __FILE__, __LINE__, \ AKANTU_DEBUG_MODULE); \ } while (false) /* -------------------------------------------------------------------------- */ #ifdef AKANTU_NDEBUG #define AKANTU_DEBUG_TEST(level) (false) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest, AKANTU_DEBUG_MODULE)) #define AKANTU_DEBUG(level, info) #define AKANTU_DEBUG_(pref, level, info) #define AKANTU_DEBUG_IN() #define AKANTU_DEBUG_OUT() #define AKANTU_DEBUG_INFO(info) #define AKANTU_DEBUG_WARNING(info) #define AKANTU_DEBUG_TRACE(info) #define AKANTU_DEBUG_ASSERT(test, info) #define AKANTU_ERROR(info) \ AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::CriticalError(), info) /* -------------------------------------------------------------------------- */ #else #define AKANTU_DEBUG(level, info) AKANTU_DEBUG_(" ", level, info) #define AKANTU_DEBUG_(pref, level, info) \ do { \ std::string _dbg_str; \ AKANTU_STRINGSTREAM_IN(_dbg_str, \ info << " " << AKANTU_LOCATION); /* NOLINT */ \ ::akantu::debug::debugger.printMessage(pref, level, _dbg_str, \ AKANTU_DEBUG_MODULE); \ } while (false) #define AKANTU_DEBUG_TEST(level) \ (::akantu::debug::debugger.testLevel(level, AKANTU_DEBUG_MODULE)) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest)) #define AKANTU_DEBUG_IN() \ AKANTU_DEBUG_( \ "==>", ::akantu::dblIn, \ __func__ \ << "()") // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay, // bugprone-lambda-function-name) #define AKANTU_DEBUG_OUT() \ AKANTU_DEBUG_( \ "<==", ::akantu::dblOut, \ __func__ \ << "()") // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay, // bugprone-lambda-function-name) #define AKANTU_DEBUG_INFO(info) AKANTU_DEBUG_("---", ::akantu::dblInfo, info) #define AKANTU_DEBUG_WARNING(info) \ AKANTU_DEBUG_("/!\\", ::akantu::dblWarning, info) #define AKANTU_DEBUG_TRACE(info) AKANTU_DEBUG_(">>>", ::akantu::dblTrace, info) #define AKANTU_DEBUG_ASSERT(test, info) \ do { \ if (not(test)) \ AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::AssertException(), \ "assert [" << #test << "] " \ << info); /* NOLINT */ \ } while (false) #define AKANTU_ERROR(info) \ do { \ AKANTU_DEBUG_("!!! ", ::akantu::dblError, info); \ AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::CriticalError(), \ info); /* NOLINT */ \ } while (false) #endif // AKANTU_NDEBUG #define AKANTU_TO_IMPLEMENT() \ AKANTU_CUSTOM_EXCEPTION_INFO( \ ::akantu::debug::NotImplementedException(), \ __func__ \ << " : not implemented yet !") // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay, // bugprone-lambda-function-name) /* -------------------------------------------------------------------------- */ namespace debug { /* ------------------------------------------------------------------------ */ template void Debugger::throwCustomException(Except ex, const std::string & info, const std::string & file, unsigned int line, const std::string & module_) const noexcept(false) { ex.setInfo(info); ex.setFile(file); ex.setLine(line); ex.setModule(module_); if (::akantu::debug::debugger.printBacktrace()) { ex.setBacktrace(::akantu::debug::getBacktrace()); } throw ex; } /* ------------------------------------------------------------------------ */ template void Debugger::throwCustomException(Except ex, const std::string & file, unsigned int line, const std::string & module_) const noexcept(false) { ex.setFile(file); ex.setLine(line); ex.setModule(module_); if (::akantu::debug::debugger.printBacktrace()) { ex.setBacktrace(::akantu::debug::getBacktrace()); } throw ex; } } // namespace debug } // namespace akantu #endif /* AKANTU_ERROR_HH_ */ diff --git a/src/io/mesh_io/mesh_io_msh.cc b/src/io/mesh_io/mesh_io_msh.cc index 7832bb016..2681903be 100644 --- a/src/io/mesh_io/mesh_io_msh.cc +++ b/src/io/mesh_io/mesh_io_msh.cc @@ -1,1125 +1,1133 @@ /** * @file mesh_io_msh.cc * * @author Dana Christen * @author Mauro Corrado * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Oct 29 2020 * * @brief Read/Write for MSH files generated by gmsh * * * @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 . * */ /* ----------------------------------------------------------------------------- Version (Legacy) 1.0 $NOD number-of-nodes node-number x-coord y-coord z-coord ... $ENDNOD $ELM number-of-elements elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list ... $ENDELM ----------------------------------------------------------------------------- Version 2.1 $MeshFormat version-number file-type data-size $EndMeshFormat $Nodes number-of-nodes node-number x-coord y-coord z-coord ... $EndNodes $Elements number-of-elements elm-number elm-type number-of-tags < tag > ... node-number-list ... $EndElements $PhysicalNames number-of-names physical-dimension physical-number "physical-name" ... $EndPhysicalNames $NodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... node-number value ... ... $EndNodeData $ElementData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number value ... ... $EndElementData $ElementNodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number number-of-nodes-per-element value ... ... $ElementEndNodeData ----------------------------------------------------------------------------- elem-type 1: 2-node line. 2: 3-node triangle. 3: 4-node quadrangle. 4: 4-node tetrahedron. 5: 8-node hexahedron. 6: 6-node prism. 7: 5-node pyramid. 8: 3-node second order line 9: 6-node second order triangle 10: 9-node second order quadrangle 11: 10-node second order tetrahedron 12: 27-node second order hexahedron 13: 18-node second order prism 14: 14-node second order pyramid 15: 1-node point. 16: 8-node second order quadrangle 17: 20-node second order hexahedron 18: 15-node second order prism 19: 13-node second order pyramid 20: 9-node third order incomplete triangle 21: 10-node third order triangle 22: 12-node fourth order incomplete triangle 23: 15-node fourth order triangle 24: 15-node fifth order incomplete triangle 25: 21-node fifth order complete triangle 26: 4-node third order edge 27: 5-node fourth order edge 28: 6-node fifth order edge 29: 20-node third order tetrahedron 30: 35-node fourth order tetrahedron 31: 56-node fifth order tetrahedron -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIOMSH::MeshIOMSH() { canReadSurface = true; canReadExtendedData = true; _msh_nodes_per_elem[_msh_not_defined] = 0; _msh_nodes_per_elem[_msh_segment_2] = 2; _msh_nodes_per_elem[_msh_triangle_3] = 3; _msh_nodes_per_elem[_msh_quadrangle_4] = 4; _msh_nodes_per_elem[_msh_tetrahedron_4] = 4; _msh_nodes_per_elem[_msh_hexahedron_8] = 8; _msh_nodes_per_elem[_msh_prism_1] = 6; _msh_nodes_per_elem[_msh_pyramid_1] = 1; _msh_nodes_per_elem[_msh_segment_3] = 3; _msh_nodes_per_elem[_msh_triangle_6] = 6; _msh_nodes_per_elem[_msh_quadrangle_9] = 9; _msh_nodes_per_elem[_msh_tetrahedron_10] = 10; _msh_nodes_per_elem[_msh_hexahedron_27] = 27; _msh_nodes_per_elem[_msh_hexahedron_20] = 20; _msh_nodes_per_elem[_msh_prism_18] = 18; _msh_nodes_per_elem[_msh_prism_15] = 15; _msh_nodes_per_elem[_msh_pyramid_14] = 14; _msh_nodes_per_elem[_msh_point] = 1; _msh_nodes_per_elem[_msh_quadrangle_8] = 8; _msh_to_akantu_element_types[_msh_not_defined] = _not_defined; _msh_to_akantu_element_types[_msh_segment_2] = _segment_2; _msh_to_akantu_element_types[_msh_triangle_3] = _triangle_3; _msh_to_akantu_element_types[_msh_quadrangle_4] = _quadrangle_4; _msh_to_akantu_element_types[_msh_tetrahedron_4] = _tetrahedron_4; _msh_to_akantu_element_types[_msh_hexahedron_8] = _hexahedron_8; _msh_to_akantu_element_types[_msh_prism_1] = _pentahedron_6; _msh_to_akantu_element_types[_msh_pyramid_1] = _not_defined; _msh_to_akantu_element_types[_msh_segment_3] = _segment_3; _msh_to_akantu_element_types[_msh_triangle_6] = _triangle_6; _msh_to_akantu_element_types[_msh_quadrangle_9] = _not_defined; _msh_to_akantu_element_types[_msh_tetrahedron_10] = _tetrahedron_10; _msh_to_akantu_element_types[_msh_hexahedron_27] = _not_defined; _msh_to_akantu_element_types[_msh_hexahedron_20] = _hexahedron_20; _msh_to_akantu_element_types[_msh_prism_18] = _not_defined; _msh_to_akantu_element_types[_msh_prism_15] = _pentahedron_15; _msh_to_akantu_element_types[_msh_pyramid_14] = _not_defined; _msh_to_akantu_element_types[_msh_point] = _point_1; _msh_to_akantu_element_types[_msh_quadrangle_8] = _quadrangle_8; _akantu_to_msh_element_types[_not_defined] = _msh_not_defined; _akantu_to_msh_element_types[_segment_2] = _msh_segment_2; _akantu_to_msh_element_types[_segment_3] = _msh_segment_3; _akantu_to_msh_element_types[_triangle_3] = _msh_triangle_3; _akantu_to_msh_element_types[_triangle_6] = _msh_triangle_6; _akantu_to_msh_element_types[_tetrahedron_4] = _msh_tetrahedron_4; _akantu_to_msh_element_types[_tetrahedron_10] = _msh_tetrahedron_10; _akantu_to_msh_element_types[_quadrangle_4] = _msh_quadrangle_4; _akantu_to_msh_element_types[_quadrangle_8] = _msh_quadrangle_8; _akantu_to_msh_element_types[_hexahedron_8] = _msh_hexahedron_8; _akantu_to_msh_element_types[_hexahedron_20] = _msh_hexahedron_20; _akantu_to_msh_element_types[_pentahedron_6] = _msh_prism_1; _akantu_to_msh_element_types[_pentahedron_15] = _msh_prism_15; _akantu_to_msh_element_types[_point_1] = _msh_point; #if defined(AKANTU_STRUCTURAL_MECHANICS) _akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2; _akantu_to_msh_element_types[_bernoulli_beam_3] = _msh_segment_2; _akantu_to_msh_element_types[_discrete_kirchhoff_triangle_18] = _msh_triangle_3; #endif std::map::iterator it; for (it = _akantu_to_msh_element_types.begin(); it != _akantu_to_msh_element_types.end(); ++it) { UInt nb_nodes = _msh_nodes_per_elem[it->second]; std::vector tmp(nb_nodes); for (UInt i = 0; i < nb_nodes; ++i) { tmp[i] = i; } switch (it->first) { case _tetrahedron_10: tmp[8] = 9; tmp[9] = 8; break; case _pentahedron_6: tmp[0] = 2; tmp[1] = 0; tmp[2] = 1; tmp[3] = 5; tmp[4] = 3; tmp[5] = 4; break; case _pentahedron_15: tmp[0] = 2; tmp[1] = 0; tmp[2] = 1; tmp[3] = 5; tmp[4] = 3; tmp[5] = 4; tmp[6] = 8; tmp[8] = 11; tmp[9] = 6; tmp[10] = 9; tmp[11] = 10; tmp[12] = 14; tmp[14] = 12; break; case _hexahedron_20: tmp[9] = 11; tmp[10] = 12; tmp[11] = 9; tmp[12] = 13; tmp[13] = 10; tmp[17] = 19; tmp[18] = 17; tmp[19] = 18; break; default: // nothing to change break; } _read_order[it->first] = tmp; } } /* -------------------------------------------------------------------------- */ MeshIOMSH::~MeshIOMSH() = default; /* -------------------------------------------------------------------------- */ namespace { struct File { std::string filename; std::ifstream infile; std::string line; size_t current_line{0}; size_t first_node_number{std::numeric_limits::max()}, last_node_number{0}; bool has_physical_names{false}; std::unordered_map node_tags; std::unordered_map element_tags; - double version{0}; + int version{0}; int size_of_size_t{0}; Mesh & mesh; MeshAccessor mesh_accessor; std::multimap, int> entity_tag_to_physical_tags; File(const std::string & filename, Mesh & mesh) : filename(filename), mesh(mesh), mesh_accessor(mesh) { infile.open(filename.c_str()); if (not infile.good()) { AKANTU_EXCEPTION("Cannot open file " << filename); } } ~File() { infile.close(); } auto good() { return infile.good(); } std::stringstream get_line() { std::string tmp_str; if (infile.eof()) { AKANTU_EXCEPTION("Reached the end of the file " << filename); } std::getline(infile, tmp_str); line = trim(tmp_str); ++current_line; return std::stringstream(line); } template void read_line(Ts &&... ts) { auto && sstr = get_line(); (void)std::initializer_list{ (sstr >> std::forward(ts), 0)...}; } }; } // namespace /* -------------------------------------------------------------------------- */ template void MeshIOMSH::populateReaders2(File & file, Readers & readers) { readers["$NOD"] = readers["$Nodes"] = [&](const std::string & /*unused*/) { UInt nb_nodes; file.read_line(nb_nodes); Array & nodes = file.mesh_accessor.getNodes(); nodes.resize(nb_nodes); file.mesh_accessor.setNbGlobalNodes(nb_nodes); size_t index; Vector coord(3); /// for each node, read the coordinates for (auto && data : enumerate(make_view(nodes, nodes.getNbComponent()))) { file.read_line(index, coord(0), coord(1), coord(2)); if (index > std::numeric_limits::max()) { AKANTU_EXCEPTION( "There are more nodes in this files than the index type of akantu " "can handle, consider recompiling with a bigger index type"); } file.first_node_number = std::min(file.first_node_number, index); file.last_node_number = std::max(file.last_node_number, index); for (auto && coord_data : zip(std::get<1>(data), coord)) { std::get<0>(coord_data) = std::get<1>(coord_data); } file.node_tags[index] = std::get<0>(data); } }; readers["$ELM"] = readers["$Elements"] = [&](const std::string & /*unused*/) { UInt nb_elements; file.read_line(nb_elements); Int index; UInt msh_type; ElementType akantu_type; for (UInt i = 0; i < nb_elements; ++i) { auto && sstr_elem = file.get_line(); sstr_elem >> index; sstr_elem >> msh_type; /// get the connectivity vector depending on the element type akantu_type = this->_msh_to_akantu_element_types[MSHElementType(msh_type)]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << msh_type << " at line " << file.current_line); continue; } Element elem{akantu_type, 0, _not_ghost}; auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type); auto node_per_element = connectivity.getNbComponent(); auto & read_order = this->_read_order[akantu_type]; /// read tags informations - if (file.version < 2) { + if (file.version < 2000) { Int tag0; Int tag1; Int nb_nodes; // reg-phys, reg-elem, number-of-nodes sstr_elem >> tag0 >> tag1 >> nb_nodes; auto & data0 = file.mesh_accessor.template getData("tag_0", akantu_type); data0.push_back(tag0); auto & data1 = file.mesh_accessor.template getData("tag_1", akantu_type); data1.push_back(tag1); - } else if (file.version < 4) { + } else if (file.version < 4000) { UInt nb_tags; sstr_elem >> nb_tags; for (UInt j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; auto & data = file.mesh_accessor.template getData( "tag_" + std::to_string(j), akantu_type); data.push_back(tag); } } Vector local_connect(node_per_element); for (UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= file.last_node_number, "Node number not in range : line " << file.current_line); local_connect(read_order[j]) = file.node_tags[node_index]; } connectivity.push_back(local_connect); elem.element = connectivity.size() - 1; file.element_tags[index] = elem; } }; readers["$Periodic"] = [&](const std::string & /*unused*/) { UInt nb_periodic_entities; file.read_line(nb_periodic_entities); file.mesh_accessor.getNodesFlags().resize(file.mesh.getNbNodes(), NodeFlag::_normal); for (UInt p = 0; p < nb_periodic_entities; ++p) { // dimension slave-tag master-tag UInt dimension; file.read_line(dimension); // transformation file.get_line(); // nb nodes UInt nb_nodes; file.read_line(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { // slave master auto && sstr = file.get_line(); // The info in the mesh seem inconsistent so they are ignored for now. continue; if (dimension == file.mesh.getSpatialDimension() - 1) { UInt slave; UInt master; sstr >> slave; sstr >> master; file.mesh_accessor.addPeriodicSlave(file.node_tags[slave], file.node_tags[master]); } } } // mesh_accessor.markMeshPeriodic(); }; } /* -------------------------------------------------------------------------- */ template void MeshIOMSH::populateReaders4(File & file, Readers & readers) { static std::map entity_type{ {0, "points"}, {1, "curve"}, {2, "surface"}, {3, "volume"}, }; readers["$Entities"] = [&](const std::string & /*unused*/) { size_t num_entity[4]; file.read_line(num_entity[0], num_entity[1], num_entity[2], num_entity[3]); for (auto entity_dim : arange(4)) { for (auto _ [[gnu::unused]] : arange(num_entity[entity_dim])) { auto && sstr = file.get_line(); int tag; double min_x; double min_y; double min_z; double max_x; double max_y; double max_z; size_t num_physical_tags; sstr >> tag >> min_x >> min_y >> min_z; - if (entity_dim > 0 or file.version < 4.1) { + if (entity_dim > 0 or file.version < 4001) { sstr >> max_x >> max_y >> max_z; } sstr >> num_physical_tags; for (auto _ [[gnu::unused]] : arange(num_physical_tags)) { int phys_tag; sstr >> phys_tag; std::string physical_name; if (this->physical_names.find(phys_tag) == this->physical_names.end()) { physical_name = "msh_block_" + std::to_string(phys_tag); } else { physical_name = this->physical_names[phys_tag]; } if (not file.mesh.elementGroupExists(physical_name)) { file.mesh.createElementGroup(physical_name, entity_dim); } else { file.mesh.getElementGroup(physical_name).addDimension(entity_dim); } file.entity_tag_to_physical_tags.insert( std::make_pair(std::make_pair(tag, entity_dim), phys_tag)); } } } }; readers["$Nodes"] = [&](const std::string & /*unused*/) { size_t num_blocks; size_t num_nodes; - if (file.version >= 4.1) { + if (file.version >= 4001) { file.read_line(num_blocks, num_nodes, file.first_node_number, file.last_node_number); } else { file.read_line(num_blocks, num_nodes); } auto & nodes = file.mesh_accessor.getNodes(); nodes.reserve(num_nodes); file.mesh_accessor.setNbGlobalNodes(num_nodes); if (num_nodes > std::numeric_limits::max()) { AKANTU_EXCEPTION( "There are more nodes in this files than the index type of akantu " "can handle, consider recompiling with a bigger index type"); } size_t node_id{0}; for (auto block [[gnu::unused]] : arange(num_blocks)) { int entity_dim; int entity_tag; int parametric; size_t num_nodes_in_block; Vector pos(3); Vector real_pos(nodes.getNbComponent()); - if (file.version >= 4.1) { + if (file.version >= 4001) { file.read_line(entity_dim, entity_tag, parametric, num_nodes_in_block); if (parametric) { AKANTU_EXCEPTION( "Akantu does not support parametric nodes in msh files"); } for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { size_t tag; file.read_line(tag); file.node_tags[tag] = node_id; ++node_id; } for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { file.read_line(pos(_x), pos(_y), pos(_z)); for (auto && data : zip(real_pos, pos)) { std::get<0>(data) = std::get<1>(data); } nodes.push_back(real_pos); } } else { file.read_line(entity_tag, entity_dim, parametric, num_nodes_in_block); for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { size_t tag; file.read_line(tag, pos(_x), pos(_y), pos(_z)); - if (file.version < 4.1) { - file.first_node_number = std::min(file.first_node_number, tag); - file.last_node_number = std::max(file.last_node_number, tag); - } + file.first_node_number = std::min(file.first_node_number, tag); + file.last_node_number = std::max(file.last_node_number, tag); for (auto && data : zip(real_pos, pos)) { std::get<0>(data) = std::get<1>(data); } nodes.push_back(real_pos); file.node_tags[tag] = node_id; ++node_id; } } } }; readers["$Elements"] = [&](const std::string & /*unused*/) { size_t num_blocks; size_t num_elements; file.read_line(num_blocks, num_elements); for (auto block [[gnu::unused]] : arange(num_blocks)) { int entity_dim; int entity_tag; int element_type; size_t num_elements_in_block; - if (file.version >= 4.1) { + if (file.version >= 4001) { file.read_line(entity_dim, entity_tag, element_type, num_elements_in_block); } else { file.read_line(entity_tag, entity_dim, element_type, num_elements_in_block); } /// get the connectivity vector depending on the element type auto && akantu_type = this->_msh_to_akantu_element_types[(MSHElementType)element_type]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << element_type << " at line " << file.current_line); continue; } Element elem{akantu_type, 0, _not_ghost}; auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type); Vector local_connect(connectivity.getNbComponent()); auto && read_order = this->_read_order[akantu_type]; auto & data0 = file.mesh_accessor.template getData("tag_0", akantu_type); data0.resize(data0.size() + num_elements_in_block, 0); + auto range = file.entity_tag_to_physical_tags.equal_range( + std::make_pair(entity_tag, entity_dim)); + auto & physical_data = file.mesh_accessor.template getData( "physical_names", akantu_type); physical_data.resize(physical_data.size() + num_elements_in_block, ""); for (auto _ [[gnu::unused]] : arange(num_elements_in_block)) { auto && sstr_elem = file.get_line(); size_t elem_tag; sstr_elem >> elem_tag; for (auto && c : arange(connectivity.getNbComponent())) { size_t node_tag; sstr_elem >> node_tag; + AKANTU_DEBUG_ASSERT(node_tag >= file.first_node_number, + "Node number not in range : line " + << file.current_line); + AKANTU_DEBUG_ASSERT(node_tag <= file.last_node_number, "Node number not in range : line " << file.current_line); node_tag = file.node_tags[node_tag]; local_connect(read_order[c]) = node_tag; } connectivity.push_back(local_connect); elem.element = connectivity.size() - 1; file.element_tags[elem_tag] = elem; - auto range = file.entity_tag_to_physical_tags.equal_range( - std::make_pair(entity_tag, entity_dim)); bool first = true; for (auto it = range.first; it != range.second; ++it) { auto phys_it = this->physical_names.find(it->second); if (first) { data0(elem.element) = it->second; // for compatibility with version 2 if (phys_it != this->physical_names.end()) { physical_data(elem.element) = phys_it->second; } first = false; } if (phys_it != this->physical_names.end()) { file.mesh.getElementGroup(phys_it->second).add(elem, true, false); } } } } for (auto && element_group : file.mesh.iterateElementGroups()) { element_group.getNodeGroup().optimize(); } }; } /* -------------------------------------------------------------------------- */ void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { File file(filename, mesh); std::map> readers; readers["$MeshFormat"] = [&](const std::string & /*unused*/) { auto && sstr = file.get_line(); + double version; int format; - sstr >> file.version >> format; + sstr >> version >> format; + + int major = std::trunc(version); + int minor = std::round(10 * (version - major)); + file.version = major * 1000 + minor; if (format != 0) { AKANTU_ERROR("This reader can only read ASCII files."); } - if (file.version > 2) { + if (file.version > 2000) { sstr >> file.size_of_size_t; if (file.size_of_size_t > int(sizeof(UInt))) { AKANTU_DEBUG_INFO("The size of the indexes in akantu might be to small " "to read this file (akantu " << sizeof(UInt) << " vs. msh file " << file.size_of_size_t << ")"); } } - if (file.version < 4) { + if (file.version < 4000) { this->populateReaders2(file, readers); } else { this->populateReaders4(file, readers); } }; auto && read_data = [&](auto && entity_tags, auto && get_data, auto && read_data) { auto read_data_tags = [&](auto x) { UInt number_of_tags{0}; file.read_line(number_of_tags); std::vector tags(number_of_tags); for (auto && tag : tags) { file.read_line(tag); } return tags; }; auto && string_tags = read_data_tags(std::string{}); auto && real_tags [[gnu::unused]] = read_data_tags(double{}); auto && int_tags = read_data_tags(int{}); for (auto & s : string_tags) { s = trim(s, '"'); } auto id = string_tags[0]; auto size = int_tags[2]; auto nb_component = int_tags[1]; auto & data = get_data(id, size, nb_component); for (auto n [[gnu::unused]] : arange(size)) { auto && sstr = file.get_line(); size_t tag; sstr >> tag; const auto & entity = entity_tags[tag]; read_data(entity, sstr, data, nb_component); } }; readers["$NodeData"] = [&](const std::string & /*unused*/) { /* $NodeData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... nodeTag(size_t) value(double) ... ... $EndNodeData */ read_data( file.node_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> Array & { auto & data = file.mesh.template getNodalData(id, nb_component); data.resize(size); return data; }, [&](auto && node, auto && sstr, auto && data, auto && nb_component [[gnu::unused]]) { for (auto c : arange(nb_component)) { sstr >> data(node, c); } }); }; readers["$ElementData"] = [&](const std::string & /*unused*/) { /* $ElementData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... elementTag(size_t) value(double) ... ... $EndElementData */ read_data( file.element_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> ElementTypeMapArray & { file.mesh.template getElementalData(id); return file.mesh.template getElementalData(id); }, [&](auto && element, auto && sstr, auto && data, auto && nb_component) { if (not data.exists(element.type)) { data.alloc(mesh.getNbElement(element.type), nb_component, element.type, element.ghost_type); } auto & data_array = data(element.type); for (auto c : arange(nb_component)) { sstr >> data_array(element.element, c); } }); }; readers["$ElementNodeData"] = [&](const std::string & /*unused*/) { /* $ElementNodeData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... elementTag(size_t) value(double) ... ... $EndElementNodeData */ read_data( file.element_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> ElementTypeMapArray & { file.mesh.template getElementalData(id); auto & data = file.mesh.template getElementalData(id); data.isNodal(true); return data; }, [&](auto && element, auto && sstr, auto && data, auto && nb_component) { int nb_nodes_per_element; sstr >> nb_nodes_per_element; if (not data.exists(element.type)) { data.alloc(mesh.getNbElement(element.type), nb_component * nb_nodes_per_element, element.type, element.ghost_type); } auto & data_array = data(element.type); for (auto c : arange(nb_component)) { sstr >> data_array(element.element, c); } }); }; readers["$PhysicalNames"] = [&](const std::string & /*unused*/) { file.has_physical_names = true; int num_of_phys_names; file.read_line(num_of_phys_names); /// the format line for (auto k [[gnu::unused]] : arange(num_of_phys_names)) { int phys_name_id; int phys_dim; std::string phys_name; file.read_line(phys_dim, phys_name_id, std::quoted(phys_name)); this->physical_names[phys_name_id] = phys_name; } }; readers["Unsupported"] = [&](const std::string & _block) { std::string block = _block.substr(1); AKANTU_DEBUG_WARNING("Unsupported block_kind " << block << " at line " << file.current_line); auto && end_block = "$End" + block; while (file.line != end_block) { file.get_line(); } }; while (file.good()) { std::string block; file.read_line(block); auto && it = readers.find(block); if (it != readers.end()) { it->second(block); std::string end_block; file.read_line(end_block); block = block.substr(1); if (end_block != "$End" + block) { AKANTU_EXCEPTION("The reader failed to properly read the block " << block << ". Expected a $End" << block << " at line " << file.current_line); } } else if (not block.empty()) { readers["Unsupported"](block); } } // mesh.updateTypesOffsets(_not_ghost); - if (file.version < 4) { + if (file.version < 4000) { this->constructPhysicalNames("tag_0", mesh); if (file.has_physical_names) { mesh.createGroupsFromMeshData("physical_names"); } } MeshUtils::fillElementToSubElementsData(mesh); } /* -------------------------------------------------------------------------- */ void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) { std::ofstream outfile; const Array & nodes = mesh.getNodes(); outfile.open(filename.c_str()); outfile << "$MeshFormat" << "\n"; outfile << "2.2 0 8" << "\n"; outfile << "$EndMeshFormat" << "\n"; outfile << std::setprecision(std::numeric_limits::digits10); outfile << "$Nodes" << "\n"; outfile << nodes.size() << "\n"; outfile << std::uppercase; for (UInt i = 0; i < nodes.size(); ++i) { Int offset = i * nodes.getNbComponent(); outfile << i + 1; for (UInt j = 0; j < nodes.getNbComponent(); ++j) { outfile << " " << nodes.storage()[offset + j]; } for (UInt p = nodes.getNbComponent(); p < 3; ++p) { outfile << " " << 0.; } outfile << "\n"; ; } outfile << std::nouppercase; outfile << "$EndNodes" << "\n"; outfile << "$Elements" << "\n"; Int nb_elements = 0; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const Array & connectivity = mesh.getConnectivity(type, _not_ghost); nb_elements += connectivity.size(); } outfile << nb_elements << "\n"; std::map element_to_msh_element; UInt element_idx = 1; auto element = ElementNull; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const auto & connectivity = mesh.getConnectivity(type, _not_ghost); element.type = type; UInt * tag[2] = {nullptr, nullptr}; if (mesh.hasData("tag_0", type, _not_ghost)) { const auto & data_tag_0 = mesh.getData("tag_0", type, _not_ghost); tag[0] = data_tag_0.storage(); } if (mesh.hasData("tag_1", type, _not_ghost)) { const auto & data_tag_1 = mesh.getData("tag_1", type, _not_ghost); tag[1] = data_tag_1.storage(); } for (auto && data : enumerate(make_view(connectivity, connectivity.getNbComponent()))) { element.element = std::get<0>(data); const auto & conn = std::get<1>(data); element_to_msh_element.insert(std::make_pair(element, element_idx)); outfile << element_idx << " " << _akantu_to_msh_element_types[type] << " 2"; /// \todo write the real data in the file for (UInt t = 0; t < 2; ++t) { if (tag[t] != nullptr) { outfile << " " << tag[t][element.element]; } else { outfile << " 0"; } } for (auto && c : conn) { outfile << " " << c + 1; } outfile << "\n"; element_idx++; } } outfile << "$EndElements" << "\n"; if (mesh.hasData(MeshDataType::_nodal)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { auto type = mesh.getTypeCode(tag, MeshDataType::_nodal); if (type != MeshDataTypeCode::_real) { AKANTU_DEBUG_WARNING( "The field " << tag << " is ignored by the MSH writer, msh files do not support " << type << " data"); continue; } auto && data = mesh.getNodalData(tag); outfile << "$NodeData" << "\n"; outfile << "1" << "\n"; outfile << "\"" << tag << "\"\n"; outfile << "1\n0.0" << "\n"; outfile << "3\n0" << "\n"; outfile << data.getNbComponent() << "\n"; outfile << data.size() << "\n"; for (auto && d : enumerate(make_view(data, data.getNbComponent()))) { outfile << std::get<0>(d) + 1; for (auto && v : std::get<1>(d)) { outfile << " " << v; } outfile << "\n"; } outfile << "$EndNodeData" << "\n"; } } if (mesh.hasData(MeshDataType::_elemental)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { auto && data = mesh.getElementalData(tag); auto type = mesh.getTypeCode(tag, MeshDataType::_elemental); if (type != MeshDataTypeCode::_real) { AKANTU_DEBUG_WARNING( "The field " << tag << " is ignored by the MSH writer, msh files do not support " << type << " data"); continue; } if (data.isNodal()) { continue; } auto size = data.size(); if (size == 0) { continue; } auto && nb_components = data.getNbComponents(); auto nb_component = nb_components(*(data.elementTypes().begin())); outfile << "$ElementData" << "\n"; outfile << "1" << "\n"; outfile << "\"" << tag << "\"\n"; outfile << "1\n0.0" << "\n"; outfile << "3\n0" << "\n"; outfile << nb_component << "\n"; outfile << size << "\n"; Element element; for (auto type : data.elementTypes()) { element.type = type; for (auto && _ : enumerate(make_view(data(type), nb_components(type)))) { element.element = std::get<0>(_); outfile << element_to_msh_element[element]; for (auto && v : std::get<1>(_)) { outfile << " " << v; } outfile << "\n"; } } outfile << "$EndElementData" << "\n"; } } outfile.close(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc index 380d2a8cd..849ddeb67 100644 --- a/src/mesh/element_group.cc +++ b/src/mesh/element_group.cc @@ -1,209 +1,205 @@ /** * @file element_group.cc * * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Dec 09 2020 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * @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_csr.hh" #include "dumpable.hh" #include "dumpable_inline_impl.hh" #include "group_manager.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_utils.hh" #include #include #include #include "element_group.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" namespace akantu { /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const std::string & group_name, const Mesh & mesh, NodeGroup & node_group, UInt dimension, const std::string & id) : mesh(mesh), name(group_name), elements("elements", id), node_group(node_group), dimension(dimension) { AKANTU_DEBUG_IN(); this->registerDumper("paraview_" + group_name, group_name, true); this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), _all_dimensions); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const ElementGroup & /*other*/) = default; /* -------------------------------------------------------------------------- */ void ElementGroup::clear() { elements.free(); } /* -------------------------------------------------------------------------- */ void ElementGroup::clear(ElementType type, GhostType ghost_type) { if (elements.exists(type, ghost_type)) { elements(type, ghost_type).clear(); } } /* -------------------------------------------------------------------------- */ bool ElementGroup::empty() const { return elements.empty(); } /* -------------------------------------------------------------------------- */ void ElementGroup::append(const ElementGroup & other_group) { AKANTU_DEBUG_IN(); node_group.append(other_group.node_group); /// loop on all element types in all dimensions for (auto ghost_type : ghost_types) { for (auto type : other_group.elementTypes(_ghost_type = ghost_type, _element_kind = _ek_not_defined)) { const Array & other_elem_list = other_group.elements(type, ghost_type); UInt nb_other_elem = other_elem_list.size(); Array * elem_list; UInt nb_elem = 0; /// create current type if doesn't exists, otherwise get information if (elements.exists(type, ghost_type)) { elem_list = &elements(type, ghost_type); nb_elem = elem_list->size(); } else { elem_list = &(elements.alloc(0, 1, type, ghost_type)); } /// append new elements to current list elem_list->resize(nb_elem + nb_other_elem); std::copy(other_elem_list.begin(), other_elem_list.end(), elem_list->begin() + nb_elem); - - /// remove duplicates - std::sort(elem_list->begin(), elem_list->end()); - Array::iterator<> end = - std::unique(elem_list->begin(), elem_list->end()); - elem_list->resize(end - elem_list->begin()); } } + this->optimize(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) { ; } stream << space << "ElementGroup [" << std::endl; stream << space << " + name: " << name << std::endl; stream << space << " + dimension: " << dimension << std::endl; elements.printself(stream, indent + 1); node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void ElementGroup::optimize() { // increasing the locality of data when iterating on the element of a group for (auto ghost_type : ghost_types) { for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { Array & els = elements(type, ghost_type); std::sort(els.begin(), els.end()); Array::iterator<> end = std::unique(els.begin(), els.end()); els.resize(end - els.begin()); } } node_group.optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::fillFromNodeGroup() { CSR node_to_elem; MeshUtils::buildNode2Elements(this->mesh, node_to_elem, this->dimension); std::set seen; Array::const_iterator<> itn = this->node_group.begin(); Array::const_iterator<> endn = this->node_group.end(); for (; itn != endn; ++itn) { CSR::iterator ite = node_to_elem.begin(*itn); CSR::iterator ende = node_to_elem.end(*itn); for (; ite != ende; ++ite) { const Element & elem = *ite; if (this->dimension != _all_dimensions && this->dimension != Mesh::getSpatialDimension(elem.type)) { continue; } if (seen.find(elem) != seen.end()) { continue; } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(elem.type); Array::const_iterator> conn_it = this->mesh.getConnectivity(elem.type, elem.ghost_type) .begin(nb_nodes_per_element); const Vector & conn = conn_it[elem.element]; UInt count = 0; for (UInt n = 0; n < conn.size(); ++n) { count += (this->node_group.getNodes().find(conn(n)) != UInt(-1) ? 1 : 0); } if (count == nb_nodes_per_element) { this->add(elem); } seen.insert(elem); } } this->optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::addDimension(UInt dimension) { this->dimension = std::max(dimension, this->dimension); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc index 68d0f9f81..3a2099ad0 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,76 +1,76 @@ /** * @file coupler_solid_cohesive_contact.cc * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Mon Jan 21 2019 * @date last modification: Wed Jun 23 2021 * * @brief class for coupling of solid mechanics cohesive and conatct mechanics * model * * * @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 "coupler_solid_cohesive_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template <> CouplerSolidContactTemplate:: CouplerSolidContactTemplate(Mesh & mesh, UInt dim, const ID & id, - std::shared_ptr dof_manager, - ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id) { + std::shared_ptr dof_manager) + : Model(mesh, ModelType::_coupler_solid_cohesive_contact, dof_manager, dim, + id) { this->mesh.registerDumper("coupler_solid_cohesive_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_cohesive_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); this->registerDataAccessor(*this); solid = std::make_unique( mesh, Model::spatial_dimension, "solid_mechanics_model_cohesive", this->dof_manager); contact = std::make_unique(mesh.getMeshFacets(), Model::spatial_dimension, "contact_mechanics_model"); } /* -------------------------------------------------------------------------- */ template <> void CouplerSolidContactTemplate::initFullImpl( const ModelOptions & options) { Model::initFullImpl(options); const auto & cscc_options = aka::as_type(options); solid->initFull(_analysis_method = cscc_options.analysis_method, _is_extrinsic = cscc_options.is_extrinsic); contact->initFull(_analysis_method = cscc_options.analysis_method); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index 5dbedb121..5f111e9e2 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,69 +1,69 @@ /** * @file coupler_solid_contact.cc * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Mon Jan 21 2019 * @date last modification: Wed Jun 23 2021 * * @brief class for coupling of solid mechanics and conatct mechanics * model * * * @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 "coupler_solid_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template <> CouplerSolidContactTemplate::CouplerSolidContactTemplate( Mesh & mesh, UInt dim, const ID & id, - std::shared_ptr dof_manager, ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id) { + std::shared_ptr dof_manager) + : Model(mesh, ModelType::_coupler_solid_contact, dof_manager, dim, id) { this->mesh.registerDumper("coupler_solid_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); this->registerDataAccessor(*this); solid = std::make_unique(mesh, Model::spatial_dimension, "solid_mechanics_model", this->dof_manager); contact = std::make_unique( mesh, Model::spatial_dimension, "contact_mechanics_model", this->dof_manager); } /* -------------------------------------------------------------------------- */ template <> void CouplerSolidContactTemplate::initFullImpl( const ModelOptions & options) { Model::initFullImpl(options); solid->initFull(_analysis_method = this->method); contact->initFull(_analysis_method = this->method); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.hh b/src/model/model_couplers/coupler_solid_contact.hh index d82ffbec3..6c04b4770 100644 --- a/src/model/model_couplers/coupler_solid_contact.hh +++ b/src/model/model_couplers/coupler_solid_contact.hh @@ -1,280 +1,276 @@ /** * @file coupler_solid_contact.hh * * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Sat Jun 26 2021 * * @brief class for coupling of solid mechanics and conatct mechanics * model in explicit * * * @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 "contact_mechanics_model.hh" #include "solid_mechanics_model.hh" #if defined(AKANTU_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive.hh" #endif /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COUPLER_SOLID_CONTACT_HH__ #define __AKANTU_COUPLER_SOLID_CONTACT_HH__ /* ------------------------------------------------------------------------ */ /* Coupling : Solid Mechanics / Contact Mechanics */ /* ------------------------------------------------------------------------ */ namespace akantu { /* -------------------------------------------------------------------------- */ template class CouplerSolidContactTemplate : public Model, public DataAccessor, public DataAccessor { static_assert( std::is_base_of::value, "SolidMechanicsModelType should be derived from SolidMechanicsModel"); /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ public: CouplerSolidContactTemplate( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "coupler_solid_contact", - std::shared_ptr dof_manager = nullptr, - ModelType model_type = std::is_same::value - ? ModelType::_coupler_solid_cohesive_contact - : ModelType::_coupler_solid_contact); + std::shared_ptr dof_manager = nullptr); ~CouplerSolidContactTemplate() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); #if defined(AKANTU_COHESIVE_ELEMENT) template ::value> * = nullptr> UInt checkCohesiveStress() { return solid->checkCohesiveStress(); } #endif template inline void applyBC(const FunctorType & func) { solid->applyBC(func); } template inline void applyBC(const FunctorType & func, const std::string & group_name) { solid->applyBC(func, group_name); } template inline void applyBC(const FunctorType & func, const ElementGroup & element_group) { solid->applyBC(func, element_group); } protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() const override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) const override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// callback for the model to instantiate the matricess when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /* ------------------------------------------------------------------------ */ /* Mass matrix for solid mechanics model */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass; } /* ------------------------------------------------------------------------ */ public: // DataAccessor UInt getNbData(const Array & /*elements*/, const SynchronizationTag & /*tag*/) const override { return 0; } void packData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) const override {} void unpackData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) override {} // DataAccessor nodes UInt getNbData(const Array & /*nodes*/, const SynchronizationTag & /*tag*/) const override { return 0; } void packData(CommunicationBuffer & /*buffer*/, const Array & /*nodes*/, const SynchronizationTag & /*tag*/) const override {} void unpackData(CommunicationBuffer & /*buffer*/, const Array & /*nodes*/, const SynchronizationTag & /*tag*/) override {} /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the solid mechanics model #if defined(AKANTU_COHESIVE_ELEMENT) template ::value> * = nullptr> SolidMechanicsModelCohesive & getSolidMechanicsModelCohesive() { return *solid; } #endif template ::value> * = nullptr> SolidMechanicsModelType & getSolidMechanicsModel() { return *solid; } /// get the contact mechanics model AKANTU_GET_MACRO(ContactMechanicsModel, *contact, ContactMechanicsModel &) /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldUInt(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; void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// solid mechanics model std::unique_ptr solid; /// contact mechanics model std::unique_ptr contact; UInt step; }; using CouplerSolidContact = CouplerSolidContactTemplate; } // namespace akantu #include "coupler_solid_contact_tmpl.hh" #endif /* __COUPLER_SOLID_CONTACT_HH__ */ diff --git a/src/model/solid_mechanics/material_selector.hh b/src/model/solid_mechanics/material_selector.hh index 5cdc2f3f4..6550a4057 100644 --- a/src/model/solid_mechanics/material_selector.hh +++ b/src/model/solid_mechanics/material_selector.hh @@ -1,159 +1,159 @@ /** * @file material_selector.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Fri Apr 09 2021 * * @brief class describing how to choose a material for a given element * * * @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 "element.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_SELECTOR_HH_ #define AKANTU_MATERIAL_SELECTOR_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModel; /** * main class to assign same or different materials for different * elements */ class MaterialSelector { public: MaterialSelector() = default; virtual ~MaterialSelector() = default; virtual inline UInt operator()(const Element & element) { if (fallback_selector) { return (*fallback_selector)(element); } return fallback_value; } inline void setFallback(UInt f) { fallback_value = f; } inline void setFallback(const std::shared_ptr & fallback_selector) { this->fallback_selector = fallback_selector; } inline std::shared_ptr & getFallbackSelector() { return this->fallback_selector; } inline UInt getFallbackValue() const { return this->fallback_value; } protected: UInt fallback_value{0}; std::shared_ptr fallback_selector; }; /* -------------------------------------------------------------------------- */ /** * class that assigns the first material to regular elements by default */ class DefaultMaterialSelector : public MaterialSelector { public: explicit DefaultMaterialSelector( const ElementTypeMapArray & material_index) : material_index(material_index) {} UInt operator()(const Element & element) override { if (not material_index.exists(element.type, element.ghost_type)) { return MaterialSelector::operator()(element); } const auto & mat_indexes = material_index(element.type, element.ghost_type); if (element.element < mat_indexes.size()) { auto && tmp_mat = mat_indexes(element.element); if (tmp_mat != UInt(-1)) { return tmp_mat; } } return MaterialSelector::operator()(element); } private: const ElementTypeMapArray & material_index; }; /* -------------------------------------------------------------------------- */ /** * Use elemental data to assign materials */ template class ElementDataMaterialSelector : public MaterialSelector { public: ElementDataMaterialSelector(const ElementTypeMapArray & element_data, const SolidMechanicsModel & model, UInt first_index = 1) : element_data(element_data), model(model), first_index(first_index) {} inline T elementData(const Element & element) { DebugLevel dbl = debug::getDebugLevel(); debug::setDebugLevel(dblError); - T data = element_data(element.type, element.ghost_type)(element.element); + T data = element_data(element); debug::setDebugLevel(dbl); return data; } inline UInt operator()(const Element & element) override; protected: /// list of element with the specified data (i.e. tag value) const ElementTypeMapArray & element_data; /// the model that the materials belong const SolidMechanicsModel & model; /// first material index: equal to 1 if none specified UInt first_index; }; /* -------------------------------------------------------------------------- */ /** * class to use mesh data information to assign different materials * where name is the tag value: tag_0, tag_1 */ template class MeshDataMaterialSelector : public ElementDataMaterialSelector { public: MeshDataMaterialSelector(const std::string & name, const SolidMechanicsModel & model, UInt first_index = 1); }; } // namespace akantu #endif /* AKANTU_MATERIAL_SELECTOR_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc index 74e2a9301..d96583923 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc @@ -1,169 +1,167 @@ /** * @file material_selector_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Dec 11 2015 * @date last modification: Fri Apr 09 2021 * * @brief Material selector for cohesive elements * * * @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 "material_selector_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DefaultMaterialCohesiveSelector::DefaultMaterialCohesiveSelector( const SolidMechanicsModelCohesive & model) : facet_material(model.getFacetMaterial()), mesh(model.getMesh()) { // backward compatibility v3: to get the former behavior back when the user // creates its own selector this->fallback_selector = std::make_shared(model.getMaterialByElement()); } /* -------------------------------------------------------------------------- */ UInt DefaultMaterialCohesiveSelector::operator()(const Element & element) { if (Mesh::getKind(element.type) == _ek_cohesive) { try { const Array & cohesive_el_to_facet = mesh.getMeshFacets().getSubelementToElement(element.type, element.ghost_type); bool third_dimension = (mesh.getSpatialDimension() == 3); const Element & facet = cohesive_el_to_facet(element.element, UInt(third_dimension)); if (facet_material.exists(facet.type, facet.ghost_type)) { return facet_material(facet.type, facet.ghost_type)(facet.element); } return fallback_value; } catch (...) { return fallback_value; } } else if (Mesh::getSpatialDimension(element.type) == mesh.getSpatialDimension() - 1) { return facet_material(element.type, element.ghost_type)(element.element); } else { return MaterialSelector::operator()(element); } } /* -------------------------------------------------------------------------- */ MeshDataMaterialCohesiveSelector::MeshDataMaterialCohesiveSelector( const SolidMechanicsModelCohesive & model) : model(model), mesh_facets(model.getMeshFacets()), material_index(mesh_facets.getData("physical_names")) { third_dimension = (model.getSpatialDimension() == 3); // backward compatibility v3: to get the former behavior back when the user // creates its own selector this->fallback_selector = std::make_shared>("physical_names", model); } /* -------------------------------------------------------------------------- */ UInt MeshDataMaterialCohesiveSelector::operator()(const Element & element) { if (Mesh::getKind(element.type) == _ek_cohesive or Mesh::getSpatialDimension(element.type) == mesh_facets.getSpatialDimension() - 1) { Element facet; if (Mesh::getKind(element.type) == _ek_cohesive) { facet = mesh_facets.getSubelementToElement(element.type, element.ghost_type)( element.element, UInt(third_dimension)); } else { facet = element; } try { std::string material_name = this->material_index(facet); return this->model.getMaterialIndex(material_name); } catch (...) { return fallback_value; } } return MaterialSelector::operator()(element); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ MaterialCohesiveRulesSelector::MaterialCohesiveRulesSelector( const SolidMechanicsModelCohesive & model, const MaterialCohesiveRules & rules, ID mesh_data_id) // what we have here is the name of model and also // the name of different materials : model(model), mesh_data_id(std::move(mesh_data_id)), mesh(model.getMesh()), mesh_facets(model.getMeshFacets()), spatial_dimension(model.getSpatialDimension()), rules(rules) { // cohesive fallback this->fallback_selector = std::make_shared(model); // non cohesive fallback this->fallback_selector->setFallback( - std::make_shared>(mesh_data_id, - model)); + std::make_shared>( + this->mesh_data_id, model)); } /* -------------------------------------------------------------------------- */ UInt MaterialCohesiveRulesSelector::operator()(const Element & element) { if (mesh_facets.getSpatialDimension(element.type) == (spatial_dimension - 1)) { const std::vector & element_to_subelement = mesh_facets.getElementToSubelement(element.type, element.ghost_type)(element.element); - // Array & facets_check = model.getFacetsCheck(); - const Element & el1 = element_to_subelement[0]; const Element & el2 = element_to_subelement[1]; ID id1 = mesh.getData(mesh_data_id, el1.type, el1.ghost_type)(el1.element); ID id2 = id1; if (el2 != ElementNull) { id2 = mesh.getData(mesh_data_id, el2.type, el2.ghost_type)(el2.element); } auto rit = rules.find(std::make_pair(id1, id2)); if (rit == rules.end()) { rit = rules.find(std::make_pair(id2, id1)); } if (rit != rules.end()) { return model.getMaterialIndex(rit->second); } } return MaterialSelector::operator()(element); } } // namespace akantu diff --git a/test/test_model/patch_tests/patch_test_linear_fixture.py b/test/test_model/patch_tests/patch_test_linear_fixture.py index 2cd441a6f..33fa4366e 100644 --- a/test/test_model/patch_tests/patch_test_linear_fixture.py +++ b/test/test_model/patch_tests/patch_test_linear_fixture.py @@ -1,137 +1,137 @@ #!/usr/bin/env python3 """ patch_test_linear_fixture.py: linear patch test in python""" __author__ = "Guillaume Anciaux and Nicolas Richart" __credits__ = [ "Guillaume Anciaux ", "Nicolas Richart ", ] __copyright__ = "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __license__ = "LGPLv3" import unittest import numpy as np import akantu class TestPatchTestLinear(unittest.TestCase): alpha = np.array([[0.01, 0.02, 0.03, 0.04], [0.05, 0.06, 0.07, 0.08], [0.09, 0.10, 0.11, 0.12]]) gradient_tolerance = 1e-13 result_tolerance = 1e-13 dofs_tolerance = 1e-15 def __init__(self, test_name, elem_type_str, functor=None): self.test_name = test_name self.functor = functor self.elem_type = akantu.getElementTypes()[elem_type_str] self.elem_type_str = elem_type_str self.dim = akantu.Mesh.getSpatialDimension(self.elem_type) super().__init__(test_name) def _internal_call(self): self.functor(self) def __getattr__(self, key): if key == self.test_name: return self._internal_call def setUp(self): self.mesh = akantu.Mesh(self.dim, self.elem_type_str) self.mesh.read(str(self.elem_type_str) + ".msh") akantu.MeshUtils.buildFacets(self.mesh) self.mesh.createBoundaryGroupFromGeometry() self.model = self.model_type(self.mesh) def tearDown(self): del self.model del self.mesh def initModel(self, method, material_file): akantu.parseInput(material_file) - akantu.setDebugLevel(akantu.dblError) + akantu.debug.setDebugLevel(akantu.dblError) self.model.initFull(method) self.applyBC() if method != akantu._static: self.model.setTimeStep(0.8 * self.model.getStableTimeStep()) def applyBC(self): boundary = self.model.getBlockedDOFs() for eg in self.mesh.iterateElementGroups(): nodes = eg.getNodeGroup().getNodes() boundary[nodes, :] = True def applyBConDOFs(self, dofs): coordinates = self.mesh.getNodes() for eg in self.mesh.iterateElementGroups(): nodes = eg.getNodeGroup().getNodes().flatten() dofs[nodes] = self.setLinearDOF(dofs[nodes], coordinates[nodes]) def prescribed_gradient(self, dof): gradient = self.alpha[:dof.shape[1], 1:self.dim + 1] return gradient def checkGradient(self, gradient, dofs): pgrad = self.prescribed_gradient(dofs).T gradient = gradient.reshape(gradient.shape[0], *pgrad.shape) diff = gradient[:] - pgrad norm = np.abs(pgrad).max() gradient_error = np.abs(diff).max() / norm self.assertAlmostEqual(0, gradient_error, delta=self.gradient_tolerance) def checkResults(self, presult_func, results, dofs): presult = presult_func(self.prescribed_gradient(dofs)).flatten() remaining_size = np.prod(np.array(results.shape[1:])) results = results.reshape((results.shape[0], remaining_size)) for result in results: diff = result - presult norm = np.abs(result).max() if norm == 0: result_error = np.abs(diff).max() else: result_error = np.abs(diff).max() / norm self.assertAlmostEqual(0., result_error, delta=self.result_tolerance) def setLinearDOF(self, dof, coord): nb_dofs = dof.shape[1] dof[:] = np.einsum('ik,ak->ai', self.alpha[:nb_dofs, 1:self.dim + 1], coord) for i in range(0, nb_dofs): dof[:, i] += self.alpha[i, 0] return dof def checkDOFs(self, dofs): coordinates = self.mesh.getNodes() ref_dofs = np.zeros_like(dofs) self.setLinearDOF(ref_dofs, coordinates) diff = dofs - ref_dofs dofs_error = np.abs(diff).max() self.assertAlmostEqual(0., dofs_error, delta=self.dofs_tolerance) @classmethod def TYPED_TEST(cls, functor, label): for type_name, _type in akantu.getElementTypes().items(): if type_name == "_point_1": continue method_name = type_name + '_' + label test_case = cls(method_name, type_name, functor) test_case.setUp() functor(test_case) test_case.tearDown()