diff --git a/packages/contact_mechanics.cmake b/packages/contact_mechanics.cmake index b814a9bb4..19c2fa2ba 100644 --- a/packages/contact_mechanics.cmake +++ b/packages/contact_mechanics.cmake @@ -1,78 +1,78 @@ #=============================================================================== # @file contact_mechanics.cmake # # @author Mohit Pundir # # @date creation: Fri Sep 03 2010 # @date last modification: Wed Jun 23 2021 # # @brief package description for contact mechanics # # # @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 . # #=============================================================================== package_declare(contact_mechanics DEPENDS model_couplers cohesive_element DESCRIPTION "Use Contact Mechanics package of Akantu") package_declare_sources(contact_mechanics model/contact_mechanics/contact_mechanics_model.hh model/contact_mechanics/contact_mechanics_model.cc model/contact_mechanics/contact_detector.hh model/contact_mechanics/contact_detector.cc model/contact_mechanics/contact_detector_inline_impl.cc - model/contact_mechanics/contact_detector_shared.hh - model/contact_mechanics/contact_detector_shared.cc - model/contact_mechanics/contact_detector_shared_inline_impl.cc model/contact_mechanics/contact_element.hh + model/contact_mechanics/contact_mechanics_shared.hh + model/contact_mechanics/contact_mechanics_shared.cc + model/contact_mechanics/contact_mechanics_shared_inline_impl.cc model/contact_mechanics/geometry_utils.hh model/contact_mechanics/geometry_utils.cc model/contact_mechanics/geometry_utils_inline_impl.cc model/contact_mechanics/resolution.hh model/contact_mechanics/resolution.cc model/contact_mechanics/resolution_utils.hh model/contact_mechanics/resolution_utils.cc model/contact_mechanics/resolutions/resolution_penalty.hh model/contact_mechanics/resolutions/resolution_penalty.cc model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc model/contact_mechanics/surface_selector.hh model/contact_mechanics/surface_selector.cc model/model_couplers/coupler_solid_contact.hh model/model_couplers/coupler_solid_contact_tmpl.hh model/model_couplers/coupler_solid_contact.cc model/model_couplers/coupler_solid_cohesive_contact.hh model/model_couplers/coupler_solid_cohesive_contact.cc model/model_couplers/cohesive_contact_solvercallback.hh model/model_couplers/cohesive_contact_solvercallback.cc ) package_declare_documentation_files(contact_mechanics manual-contactmechanicsmodel.tex manual-contact-detector.tex ) package_declare_documentation(contact_mechanics "This package activates the contact mechanics model") diff --git a/packages/contact_mechanics_internodes.cmake b/packages/contact_mechanics_internodes.cmake index 299acca79..18e820a89 100644 --- a/packages/contact_mechanics_internodes.cmake +++ b/packages/contact_mechanics_internodes.cmake @@ -1,42 +1,46 @@ #=============================================================================== # @file contact_mechanics_internodes.cmake # # @author Moritz Waldleben # # @date creation: Thu Jul 09 2022 # @date last modification: Fri Jul 10 2020 # # @brief package description for contact mechanics internodes # # # @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 . # #=============================================================================== package_declare(contact_mechanics_internodes DEPENDS implicit DESCRIPTION "Use Contact Mechanics Internodes package of Akantu") package_declare_sources(contact_mechanics_internodes model/contact_mechanics_internodes/contact_mechanics_internodes_model.hh model/contact_mechanics_internodes/contact_mechanics_internodes_model.cc model/contact_mechanics_internodes/contact_detector_internodes.hh model/contact_mechanics_internodes/contact_detector_internodes.cc + + model/model_couplers/coupler_solid_contact_internodes.hh + model/model_couplers/coupler_solid_contact_internodes.cc + model/model_couplers/coupler_solid_contact_internodes_tmpl.hh ) diff --git a/python/py_contact_mechanics_internodes_model.cc b/python/py_contact_mechanics_internodes_model.cc index af8aa65c2..c99f7d80b 100644 --- a/python/py_contact_mechanics_internodes_model.cc +++ b/python/py_contact_mechanics_internodes_model.cc @@ -1,115 +1,114 @@ /** * @file py_contact_mechanics_internodes_model.cc * * @author Moritz Waldleben * * @date creation: Thu Jul 09 2022 * @date last modification: Thu Jul 09 2022 * * @brief Contact mechanics internodes 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 /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](ContactMechanicsInternodesModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function_nocopy_detector(func_name) \ def( \ #func_name, \ [](ContactDetectorInternodes & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](ContactMechanicsInternodesModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_contact_mechanics_internodes_model(py::module & mod) { py::class_(mod, "ContactDetectorInternodes", py::multiple_inheritance()) .def(py::init(), py::arg("mesh"), py::arg("id") = "contact_detector_internodes") .def_function_nocopy_detector(getInitialMasterNodeGroup) .def_function_nocopy_detector(getInitialSlaveNodeGroup) .def_function_nocopy_detector(getMasterNodeGroup) .def_function_nocopy_detector(getSlaveNodeGroup) .def_function_nocopy_detector(getMasterRadiuses) .def_function_nocopy_detector(getSlaveRadiuses) .def( "constructInterpolationMatrix", [](ContactDetectorInternodes & self, NodeGroup & ref_node_group, NodeGroup & eval_node_group, Array eval_radiuses) -> decltype(auto) { return self.constructInterpolationMatrix(ref_node_group, eval_node_group, eval_radiuses); }, py::arg("ref_node_group"), py::arg("eval_node_group"), py::arg("eval_radiuses"), py::return_value_policy::reference); /* ------------------------------------------------------------------------ */ py::class_(mod, "ContactMechanicsInternodesModel", py::multiple_inheritance()) .def(py::init, const ModelType>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "contact_mechanics_internodes_model", py::arg("dof_manager") = nullptr, py::arg("model_type") = ModelType::_solid_mechanics_model) .def( "initFull", [](ContactMechanicsInternodesModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_function_nocopy(assembleInternodesMatrix) - .def_function_nocopy(getSolidMechanicsModel) - .def_function_nocopy(getContactDetectorInternodes) + .def_function_nocopy(getContactDetector) .def_function_nocopy(getLambdas) .def_function_nocopy(updateAfterStep) .def("assembleInterfaceMass", [](ContactMechanicsInternodesModel & self, NodeGroup & contact_node_group) -> decltype(auto) { return self.assembleInterfaceMass(contact_node_group); }, py::arg("contact_node_group"), py::return_value_policy::copy); } } // namespace akantu diff --git a/python/py_model_couplers.cc b/python/py_model_couplers.cc index 9ab064d05..606adbad0 100644 --- a/python/py_model_couplers.cc +++ b/python/py_model_couplers.cc @@ -1,121 +1,138 @@ /** - * @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 . - * - */ +* @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 /* -------------------------------------------------------------------------- */ #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>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "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); + + register_coupler_solid_contact( + mod, "CouplerSolidContactInternodes") + .def( + "getSolidMechanicsModel", + [](CouplerSolidContactInternodes & self) -> decltype(auto) { + return self.getSolidMechanicsModel(); + }, + py::return_value_policy::reference) + .def( + "initFull", + [](CouplerSolidContactInternodes & self, + const AnalysisMethod & analysis_method) { + self.initFull(_analysis_method = analysis_method); + }, + py::arg("_analysis_method") = _static); } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh index d6195644e..5aff0c42e 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,200 +1,200 @@ /** * @file contact_detector.hh * * @author Mohit Pundir * * @date creation: Mon Dec 13 2010 * @date last modification: Thu Jun 24 2021 * * @brief Mother class for all detection algorithms * * * @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_bbox.hh" #include "aka_common.hh" #include "aka_grid_dynamic.hh" -#include "contact_detector_shared.hh" #include "contact_element.hh" +#include "contact_mechanics_shared.hh" #include "element_class.hh" #include "element_group.hh" #include "fe_engine.hh" #include "geometry_utils.hh" #include "mesh.hh" #include "mesh_io.hh" #include "parsable.hh" #include "surface_selector.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_HH__ #define __AKANTU_CONTACT_DETECTOR_HH__ namespace akantu { enum class Surface { master, slave }; /* -------------------------------------------------------------------------- */ class ContactDetector : public Parsable, public AbstractContactDetector { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactDetector(Mesh & /*mesh*/, const ID & id = "contact_detector"); ContactDetector(Mesh & /*mesh*/, Array positions, const ID & id = "contact_detector"); ~ContactDetector() override = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// performs all search steps void search(Array & elements, Array & gaps, Array & normals, Array & tangents, Array & projections); /// performs global spatial search to construct spatial grids void globalSearch(SpatialGrid & /*slave_grid*/, SpatialGrid & /*master_grid*/); /// performs local search to find closet master node to a slave node void localSearch(SpatialGrid & /*slave_grid*/, SpatialGrid & /*master_grid*/); /// create contact elements void createContactElements(Array & elements, Array & gaps, Array & normals, Array & tangents, Array & projections); private: /// reads the input file to get contact detection options void parseSection(const ParserSection & section) override; /* ------------------------------------------------------------------------ */ /* Inline Methods */ /* ------------------------------------------------------------------------ */ public: /// checks whether the natural projection is valid or not inline bool checkValidityOfProjection(Vector & projection) const; /// computes the optimal cell size for grid inline void computeCellSpacing(Vector & spacing) const; /// constructs a grid containing nodes lying within bounding box inline void constructGrid(SpatialGrid & grid, BBox & bbox, const Array & nodes_list) const; /// constructs the bounding box based on nodes list inline void constructBoundingBox(BBox & bbox, const Array & nodes_list) const; /// computes the maximum in radius for a given mesh inline void computeMaximalDetectionDistance(); /// constructs the connectivity for a contact element inline Vector constructConnectivity(UInt & slave, const Element & master) const; /// computes normal on an element inline void computeNormalOnElement(const Element & element, Vector & normal) const; /// extracts vectors which forms the plane of element inline void vectorsAlongElement(const Element & el, Matrix & vectors) const; /// computes the gap between slave and its projection on master /// surface inline Real computeGap(const Vector & slave, const Vector & master) const; /// filter boundary elements inline void filterBoundaryElements(const Array & elements, Array & boundary_elements) const; /// checks whether self contact condition leads to a master element /// which is closet but not orthogonally opposite to slave surface inline bool isValidSelfContact(const UInt & slave_node, const Real & gap, const Vector & normal) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// returns the maximum detection distance AKANTU_GET_MACRO(MaximumDetectionDistance, max_dd, Real); AKANTU_SET_MACRO(MaximumDetectionDistance, max_dd, Real); /// returns the bounding box extension AKANTU_GET_MACRO(MaximumBoundingBox, max_bb, Real); AKANTU_SET_MACRO(MaximumBoundingBox, max_bb, Real); /// returns the minimum detection distance AKANTU_GET_MACRO(MinimumDetectionDistance, min_dd, Real); AKANTU_SET_MACRO(MinimumDetectionDistance, min_dd, Real); AKANTU_GET_MACRO_NOT_CONST(SurfaceSelector, *surface_selector, SurfaceSelector &); AKANTU_SET_MACRO(SurfaceSelector, surface_selector, std::shared_ptr); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// maximal detection distance for grid spacing Real max_dd; /// minimal detection distance for grid spacing Real min_dd; /// maximal bounding box extension Real max_bb; /// tolerance for finding natural projection Real projection_tolerance; /// iterations for finding natural projection UInt max_iterations; /// tolerance for extending a master elements on all sides Real extension_tolerance; /// node selector for selecting master and slave nodes std::shared_ptr surface_selector; /// contact pair slave node to closet master node std::vector> contact_pairs; /// type of detection explicit/implicit DetectionType detection_type; }; } // namespace akantu #include "contact_detector_inline_impl.cc" #endif /* __AKANTU_CONTACT_DETECTOR_HH__ */ diff --git a/src/model/contact_mechanics/contact_detector_shared.cc b/src/model/contact_mechanics/contact_detector_shared.cc deleted file mode 100644 index ffd7054fb..000000000 --- a/src/model/contact_mechanics/contact_detector_shared.cc +++ /dev/null @@ -1,11 +0,0 @@ -#include "contact_detector_shared.hh" - -namespace akantu { - -AbstractContactDetector::AbstractContactDetector(akantu::Mesh & mesh, Array initial_positions) - : mesh(mesh), - spatial_dimension(mesh.getSpatialDimension()), - positions(initial_positions) { -} - -} // namespace akantu diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 817c09f72..c3b23651f 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,704 +1,704 @@ /** * @file contact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Jul 28 2021 * * @brief Contact mechanics model * * * @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 "contact_mechanics_model.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id) { + : AbstractContactMechanicsModel(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension - 1, _not_ghost, _ek_regular); this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, id + ":contact_detector"); registerFEEngineObject("ContactFacetsFEEngine", mesh, Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (not this->parser.getLastParsedFile().empty()) { this->instantiateResolutions(); this->initResolutions(); } this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::instantiateResolutions() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); for (const auto & section : model_resolutions) { this->registerNewResolution(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); for (const auto & section : sub_sections) { this->registerNewResolution(section); } if (resolutions.empty()) { AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" << getID()); } are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ParserSection & section) { std::string res_name; std::string res_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); res_name = tmp; /** this can seem weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A contact resolution of type \'" << res_type << "\' in the input file has been defined without a name!"); } Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); res.parseSection(section); return res; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution( const ID & res_name, const ID & res_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == resolutions_names_to_id.end(), "A resolution with this name '" << res_name << "' has already been registered. " << "Please use unique names for resolutions"); UInt res_count = resolutions.size(); resolutions_names_to_id[res_name] = res_count; std::stringstream sstr_res; sstr_res << this->id << ":" << res_count << ":" << res_type; ID res_id = sstr_res.str(); std::unique_ptr resolution = ResolutionFactory::getInstance().allocate(res_type, spatial_dimension, opt_param, *this, res_id); resolutions.push_back(std::move(resolution)); return *(resolutions.back()); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initResolutions() { AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); if (!are_resolutions_instantiated) { instantiateResolutions(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { AKANTU_DEBUG_IN(); getFEEngine("ContactMechanicsModel").initShapeFunctions(_not_ghost); getFEEngine("ContactMechanicsModel").initShapeFunctions(_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver( TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType /*unused*/) { // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->normal_force, spatial_dimension, "normal_force"); this->allocNodalField(this->tangential_force, spatial_dimension, "tangential_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->contact_state, 1, "contact_state"); this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements"); this->allocNodalField(this->normals, spatial_dimension, "normals"); auto surface_dimension = spatial_dimension - 1; this->allocNodalField(this->tangents, surface_dimension * spatial_dimension, "tangents"); this->allocNodalField(this->projections, surface_dimension, "projections"); this->allocNodalField(this->previous_projections, surface_dimension, "previous_projections"); this->allocNodalField(this->previous_tangents, surface_dimension * spatial_dimension, "previous_tangents"); this->allocNodalField(this->tangential_tractions, surface_dimension, "tangential_tractions"); this->allocNodalField(this->previous_tangential_tractions, surface_dimension, "previous_tangential_tractions"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); UInt nb_nodes = mesh.getNbNodes(); this->internal_force->clear(); this->normal_force->clear(); this->tangential_force->clear(); internal_force->resize(nb_nodes, 0.); normal_force->resize(nb_nodes, 0.); tangential_force->resize(nb_nodes, 0.); // assemble the forces due to contact auto assemble = [&](auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements // AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); // assemble(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->resize(nb_nodes); internal_array->zero(); }; resize_arrays(gaps); resize_arrays(normals); resize_arrays(tangents); resize_arrays(projections); resize_arrays(tangential_tractions); resize_arrays(contact_state); resize_arrays(nodal_area); resize_arrays(external_force); this->detector->search(contact_elements, *gaps, *normals, *tangents, *projections); // interpenetration value must be positive for contact mechanics // model to work by default the gap value from detector is negative std::for_each((*gaps).begin(), (*gaps).end(), [](Real & gap) { gap *= -1.; }); if (!contact_elements.empty()) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections (*previous_projections).copy(*projections); // saving previous tangents (*previous_tangents).copy(*tangents); // saving previous tangential traction (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); previous_master_elements->set(ElementNull); for (auto & element : contact_elements) { (*previous_master_elements)[element.slave] = element.master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas(GhostType ghost_type) { UInt nb_nodes = mesh.getNbNodes(); nodal_area->resize(nb_nodes); nodal_area->zero(); external_force->resize(nb_nodes); external_force->zero(); auto & fem_boundary = getFEEngineClassBoundary("ContactMechanicsModel"); fem_boundary.initShapeFunctions(getPositions(), _not_ghost); fem_boundary.initShapeFunctions(getPositions(), _ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); IntegrationPoint quad_point; quad_point.ghost_type = ghost_type; auto & group = mesh.getElementGroup("contact_surface"); UInt nb_degree_of_freedom = external_force->getNbComponent(); for (auto && type : group.elementTypes(spatial_dimension - 1, ghost_type)) { const auto & element_ids = group.getElements(type, ghost_type); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type, ghost_type); UInt nb_elements = element_ids.size(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array dual_before_integ(nb_elements * nb_quad_points, nb_degree_of_freedom, 0.); Array quad_coords(nb_elements * nb_quad_points, spatial_dimension); const auto & normals_on_quad = fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type); auto normals_begin = normals_on_quad.begin(spatial_dimension); decltype(normals_begin) normals_iter; auto quad_coords_iter = quad_coords.begin(spatial_dimension); auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom); quad_point.type = type; Element subelement; subelement.type = type; subelement.ghost_type = ghost_type; for (auto el : element_ids) { subelement.element = el; const auto & element_to_subelement = mesh.getElementToSubelement(type)(el); Vector outside(spatial_dimension); mesh.getBarycenter(subelement, outside); Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector inside_to_outside(spatial_dimension); inside_to_outside = outside - inside; normals_iter = normals_begin + el * nb_quad_points; quad_point.element = el; for (auto q : arange(nb_quad_points)) { quad_point.num_point = q; auto ddot = inside_to_outside.dot(*normals_iter); Vector normal(*normals_iter); if (ddot < 0) { normal *= -1.0; } (*dual_iter) .mul(Matrix::eye(spatial_dimension, 1), normal); ++dual_iter; ++quad_coords_iter; ++normals_iter; } } Array dual_by_shapes(nb_elements * nb_quad_points, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, ghost_type, element_ids); Array dual_by_shapes_integ(nb_elements, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ, nb_degree_of_freedom * nb_nodes_per_element, type, ghost_type, element_ids); this->getDOFManager().assembleElementalArrayLocalArray( dual_by_shapes_integ, *external_force, type, ghost_type, 1., element_ids); } for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); Vector force(std::get<1>(tuple)); area = force.norm(); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Contact Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + resolutions [" << std::endl; for (const auto & resolution : resolutions) { resolution->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { for (auto & resolution : resolutions) { resolution->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep(bool converged) { for (auto & resolution : resolutions) { resolution->afterSolveStep(converged); } } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string & /*unused*/, const std::string & /*unused*/, bool /*unused*/) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["contact_force"] = this->internal_force.get(); real_nodal_fields["normal_force"] = this->normal_force.get(); real_nodal_fields["tangential_force"] = this->tangential_force.get(); real_nodal_fields["blocked_dofs"] = this->blocked_dofs.get(); real_nodal_fields["normals"] = this->normals.get(); real_nodal_fields["tangents"] = this->tangents.get(); real_nodal_fields["gaps"] = this->gaps.get(); real_nodal_fields["areas"] = this->nodal_area.get(); real_nodal_fields["tangential_traction"] = this->tangential_tractions.get(); std::shared_ptr field; if (padding_flag) { field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); } else { field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); } return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldUInt(const std::string & field_name, const std::string & group_name, bool /*padding_flag*/) { std::shared_ptr field; if (field_name == "contact_state") { auto && func = std::make_unique>(); field = mesh.createNodalField(this->contact_state.get(), group_name); field = dumpers::FieldComputeProxy::createFieldCompute(field, std::move(func)); } return field; } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & elements, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & dofs, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index cddce2919..8b1ac3f15 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,372 +1,370 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jun 24 2021 * * @brief Model of Contact Mechanics * * * @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 "boundary_condition.hh" #include "contact_detector.hh" #include "data_accessor.hh" -#include "fe_engine.hh" -#include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ namespace akantu { class Resolution; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -class ContactMechanicsModel : public Model, +class ContactMechanicsModel : public AbstractContactMechanicsModel, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "contact_mechanics_model", std::shared_ptr dof_manager = nullptr, ModelType model_type = ModelType::_contact_mechanics_model); ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// allocate all vectors void initSolver(TimeStepSolverType /*unused*/, NonLinearSolverType /*unused*/) override; /// initialize all internal arrays for resolutions void initResolutions(); /// initialize the modelType void initModel() override; /// call back for the solver, computes the force residual void assembleResidual() override; /// 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; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const 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; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: - void search(); + void search() override; void computeNodalAreas(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Contact Resolution */ /* ------------------------------------------------------------------------ */ public: /// register an empty contact resolution of a given type Resolution & registerNewResolution(const ID & res_name, const ID & res_type, const ID & opt_param); protected: /// register a resolution in the dynamic database Resolution & registerNewResolution(const ParserSection & res_section); /// read the resolution files to instantiate all the resolutions void instantiateResolutions(); /// save the parameters from previous state void savePreviousState(); /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* 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; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: /// contact detection class friend class ContactDetector; /// contact resolution class friend class Resolution; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the ContactMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the ContactMechanicsModel::increment vector \warn only consistent /// if ContactMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the ContactMechanics::normal_force vector (normal forces) AKANTU_GET_MACRO(NormalForce, *normal_force, Array &); /// get the ContactMechanics::tangential_force vector (friction forces) AKANTU_GET_MACRO(TangentialForce, *tangential_force, Array &); /// get the ContactMechanics::traction vector (friction traction) AKANTU_GET_MACRO(TangentialTractions, *tangential_tractions, Array &); /// get the ContactMechanics::previous_tangential_tractions vector AKANTU_GET_MACRO(PreviousTangentialTractions, *previous_tangential_tractions, Array &); /// get the ContactMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the ContactMechanics::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the ContactMechanics::gaps (contact gaps) AKANTU_GET_MACRO(Gaps, *gaps, Array &); /// get the ContactMechanics::normals (normals on slave nodes) AKANTU_GET_MACRO(Normals, *normals, Array &); /// get the ContactMechanics::tangents (tangents on slave nodes) AKANTU_GET_MACRO(Tangents, *tangents, Array &); /// get the ContactMechanics::previous_tangents (tangents on slave nodes) AKANTU_GET_MACRO(PreviousTangents, *previous_tangents, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); /// get the ContactMechanics::previous_projections (previous_projections) AKANTU_GET_MACRO(PreviousProjections, *previous_projections, Array &); /// get the ContactMechanics::projections (projections) AKANTU_GET_MACRO(Projections, *projections, Array &); /// get the ContactMechanics::contact_state vector (no_contact/stick/slip /// state) AKANTU_GET_MACRO(ContactState, *contact_state, Array &); /// get the ContactMechanics::previous_master_elements AKANTU_GET_MACRO(PreviousMasterElements, *previous_master_elements, Array &); /// get contact detector AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &); /// get the contact elements inline const Array & getContactElements() const { return contact_elements; } /// get the current positions of the nodes inline Array & getPositions() { return detector->getPositions(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array std::unique_ptr> displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// contact forces array std::unique_ptr> internal_force; /// external forces array std::unique_ptr> external_force; /// normal force array std::unique_ptr> normal_force; /// friction force array std::unique_ptr> tangential_force; /// friction traction array std::unique_ptr> tangential_tractions; /// previous friction traction array std::unique_ptr> previous_tangential_tractions; /// boundary vector std::unique_ptr> blocked_dofs; /// array to store gap between slave and master std::unique_ptr> gaps; /// array to store normals from master to slave std::unique_ptr> normals; /// array to store tangents on the master element std::unique_ptr> tangents; /// array to store previous tangents on the master element std::unique_ptr> previous_tangents; /// array to store nodal areas std::unique_ptr> nodal_area; /// array to store stick/slip state : std::unique_ptr> contact_state; /// array to store previous projections in covariant basis std::unique_ptr> previous_projections; // array to store projections in covariant basis std::unique_ptr> projections; /// contact detection std::unique_ptr detector; /// list of contact resolutions std::vector> resolutions; /// mapping between resolution name and resolution internal id std::map resolutions_names_to_id; /// array to store contact elements Array contact_elements; /// array to store previous master elements std::unique_ptr> previous_master_elements; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/contact_mechanics/contact_mechanics_shared.cc b/src/model/contact_mechanics/contact_mechanics_shared.cc new file mode 100644 index 000000000..55a310845 --- /dev/null +++ b/src/model/contact_mechanics/contact_mechanics_shared.cc @@ -0,0 +1,19 @@ +#include "contact_mechanics_shared.hh" + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +AbstractContactDetector::AbstractContactDetector(akantu::Mesh & mesh, Array initial_positions) + : mesh(mesh), + spatial_dimension(mesh.getSpatialDimension()), + positions(initial_positions) { +} + +/* -------------------------------------------------------------------------- */ +AbstractContactMechanicsModel::AbstractContactMechanicsModel( + akantu::Mesh & mesh, const akantu::ModelType & type, + std::shared_ptr dof_manager, akantu::UInt dim, + const akantu::ID id) : Model(mesh, type, dof_manager, dim, id) { +} + +} // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector_shared.hh b/src/model/contact_mechanics/contact_mechanics_shared.hh similarity index 67% rename from src/model/contact_mechanics/contact_detector_shared.hh rename to src/model/contact_mechanics/contact_mechanics_shared.hh index 91f060f23..82663cdf8 100644 --- a/src/model/contact_mechanics/contact_detector_shared.hh +++ b/src/model/contact_mechanics/contact_mechanics_shared.hh @@ -1,66 +1,84 @@ -#ifndef __AKANTU_CONTACT_DETECTOR_SHARED_HH__ -#define __AKANTU_CONTACT_DETECTOR_SHARED_HH__ +#ifndef __AKANTU_CONTACT_MECHANICS_SHARED_HH__ +#define __AKANTU_CONTACT_MECHANICS_SHARED_HH__ #include "aka_common.hh" #include "aka_grid_dynamic.hh" #include "fe_engine.hh" +#include "model.hh" + +/* -------------------------------------------------------------------------- */ +/* Shared abstractions across different contact mechanics models */ +/* -------------------------------------------------------------------------- */ namespace akantu { /// Base class for contact detectors class AbstractContactDetector { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ protected: AbstractContactDetector(Mesh & mesh, Array initial_positions); public: virtual ~AbstractContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// Return a new vector with the positions of a node. inline Vector getNodePosition(UInt node) const; protected: /// Fill the matrix with the coordinates of an element. inline void coordinatesOfElement(const Element & el, Matrix & coords) const; /// Compute the minimum and maximum element sizes. /// Make sure to call fillNodesToElements on the mesh before. template inline auto computeElementSizes(vector_type && nodes) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the mesh AKANTU_GET_MACRO(Mesh, mesh, Mesh &) AKANTU_GET_MACRO_NOT_CONST(Positions, positions, Array &); AKANTU_SET_MACRO(Positions, positions, Array); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Mesh Mesh & mesh; /// dimension of the model const UInt spatial_dimension; /// contains the updated positions of the nodes Array positions; }; +/// Base class for contact mechanics models +// TODO: evaluate whether this is _really_ useful or not... +class AbstractContactMechanicsModel : public Model { +public: + AbstractContactMechanicsModel(Mesh & mesh, const ModelType & type, + std::shared_ptr dof_manager, UInt dim, const ID id); + + /// Perform the contact search. + /// This is called by the coupler before the relevant terms are assembled, to + /// update the internal state of the contact model. + virtual void search() = 0; +}; + } // namespace akantu -#include "contact_detector_shared_inline_impl.cc" +#include "contact_mechanics_shared_inline_impl.cc" -#endif // __AKANTU_CONTACT_DETECTOR_SHARED_HH__ +#endif // __AKANTU_CONTACT_MECHANICS_SHARED_HH__ diff --git a/src/model/contact_mechanics/contact_detector_shared_inline_impl.cc b/src/model/contact_mechanics/contact_mechanics_shared_inline_impl.cc similarity index 88% rename from src/model/contact_mechanics/contact_detector_shared_inline_impl.cc rename to src/model/contact_mechanics/contact_mechanics_shared_inline_impl.cc index e38a39ccb..8982031cc 100644 --- a/src/model/contact_mechanics/contact_detector_shared_inline_impl.cc +++ b/src/model/contact_mechanics/contact_mechanics_shared_inline_impl.cc @@ -1,64 +1,64 @@ -#ifndef __AKANTU_CONTACT_DETECTOR_SHARED_INLINE_IMPL_CC__ -#define __AKANTU_CONTACT_DETECTOR_SHARED_INLINE_IMPL_CC__ +#ifndef __AKANTU_CONTACT_MECHANICS_SHARED_INLINE_IMPL_CC__ +#define __AKANTU_CONTACT_MECHANICS_SHARED_INLINE_IMPL_CC__ /* -------------------------------------------------------------------------- */ -#include "contact_detector_shared.hh" +#include "contact_mechanics_shared.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline Vector AbstractContactDetector::getNodePosition(UInt node) const { Vector position(spatial_dimension); for (UInt s : arange(spatial_dimension)) { position(s) = this->positions(node, s); } return position; } /* -------------------------------------------------------------------------- */ void AbstractContactDetector::coordinatesOfElement( const Element & el, Matrix & coords) const { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); const Vector connect = mesh.getConnectivity(el.type, _not_ghost) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s : arange(spatial_dimension)) { coords(s, n) = this->positions(node, s); } } } /* -------------------------------------------------------------------------- */ template inline auto AbstractContactDetector::computeElementSizes(vector_type && nodes) const { struct { Real min_size; Real max_size; } out_values; out_values = { std::numeric_limits::max(), std::numeric_limits::min() }; for (auto node : nodes) { Array elements; this->mesh.getAssociatedElements(node, elements); for (auto element : elements) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); Matrix elem_coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(element, elem_coords); Real elem_size = FEEngine::getElementInradius(elem_coords, element.type); out_values.max_size = std::max(out_values.max_size, elem_size); out_values.min_size = std::min(out_values.min_size, elem_size); } } return out_values; } } // namespace akantu -#endif // __AKANTU_CONTACT_DETECTOR_SHARED_INLINE_IMPL_CC__ +#endif // __AKANTU_CONTACT_MECHANICS_SHARED_INLINE_IMPL_CC__ diff --git a/src/model/contact_mechanics_internodes/contact_detector_internodes.cc b/src/model/contact_mechanics_internodes/contact_detector_internodes.cc index de7888672..9eed3d8c6 100644 --- a/src/model/contact_mechanics_internodes/contact_detector_internodes.cc +++ b/src/model/contact_mechanics_internodes/contact_detector_internodes.cc @@ -1,467 +1,467 @@ /** * @file contact_detector_internodes.cc * * @author Moritz Waldleben * * @date creation: Thu Jul 09 2022 * @date last modification: Thu Jul 17 2022 * * @brief Algorithm to detetect contact nodes * * * @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_detector_internodes.hh" -#include "contact_detector_shared.hh" +#include "contact_mechanics_shared.hh" #include "element_group.hh" #include "node_group.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactDetectorInternodes::ContactDetectorInternodes(Mesh & mesh, const ID & id) : Parsable(ParserType::_contact_detector, id), AbstractContactDetector(mesh, mesh.getNodes()) { const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); this->parseSection(section); auto & initial_master_node_group = mesh.createNodeGroup("initial_contact_master_nodes"); auto & initial_slave_node_group = mesh.createNodeGroup("initial_contact_slave_nodes"); auto & master_node_group = mesh.createNodeGroup("contact_master_nodes"); auto & slave_node_group = mesh.createNodeGroup("contact_slave_nodes"); initial_master_node_group.append( mesh.getElementGroup(id_master_nodes).getNodeGroup()); initial_slave_node_group.append(mesh.getElementGroup(id_slave_nodes).getNodeGroup()); master_node_group.append( mesh.getElementGroup(id_master_nodes).getNodeGroup()); slave_node_group.append(mesh.getElementGroup(id_slave_nodes).getNodeGroup()); // Remember to fill nodesToElements as they are needed for computeElementSizes auto surface_dimension = spatial_dimension - 1; this->mesh.fillNodesToElements(surface_dimension); auto element_sizes = computeElementSizes(arange(mesh.getNbNodes())); this->max_element_size = element_sizes.max_size; } /* -------------------------------------------------------------------------- */ void ContactDetectorInternodes::parseSection(const ParserSection & section) { this->id_master_nodes = section.getParameterValue("master"); this->id_slave_nodes = section.getParameterValue("slave"); this->relative_grid_spacing_factor = section.getParameter("grid_spacing", 3.0); this->relative_penetration_tolerance = section.getParameter("penetration_tolerance", 0.1); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getInitialMasterNodeGroup() { return mesh.getNodeGroup("initial_contact_master_nodes"); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getInitialSlaveNodeGroup() { return mesh.getNodeGroup("initial_contact_slave_nodes"); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getMasterNodeGroup() { return mesh.getNodeGroup("contact_master_nodes"); } /* -------------------------------------------------------------------------- */ NodeGroup & ContactDetectorInternodes::getSlaveNodeGroup() { return mesh.getNodeGroup("contact_slave_nodes"); } /* -------------------------------------------------------------------------- */ void ContactDetectorInternodes::findContactNodes(NodeGroup & master_node_group, NodeGroup & slave_node_group) { // We need to find master->slave and slave->master interpolation nodes and // radii that satisfy the conditions from [1], page 51, equations (2) and (3). Real grid_spacing = this->max_element_size * this->relative_grid_spacing_factor; bool still_isolated_nodes = true; while (still_isolated_nodes) { auto master_grid = constructGrid(master_node_group, grid_spacing); auto slave_grid = constructGrid(slave_node_group, grid_spacing); // computeRadiuses will produce a result satisfying equation (2). auto && nb_slave_nodes_inside_radius = computeRadiuses(master_radiuses, master_node_group, master_grid, slave_node_group, slave_grid); auto && nb_master_nodes_inside_radius = computeRadiuses(slave_radiuses, slave_node_group, slave_grid, master_node_group, master_grid); // Check equation (3): if a node is still isolated, remove it and iterate. still_isolated_nodes = false; if (master_node_group.applyNodeFilter([&](auto && node) { return nb_master_nodes_inside_radius[node] > 0; }) > 0) { still_isolated_nodes = true; } if (slave_node_group.applyNodeFilter([&](auto && node) { return nb_slave_nodes_inside_radius[node] > 0; }) > 0) { still_isolated_nodes = true; } // To keep removals in O(1), removing from a NodeGroup might reorder nodes master_node_group.optimize(); slave_node_group.optimize(); } // Check that equation (3) is satisfied. if (DEBUG_VERIFY_INTERPOLATION_CONDITIONS) { for (UInt eval_node : slave_node_group) { bool ok = false; for (auto entry : enumerate(master_node_group)) { auto j = std::get<0>(entry); auto ref_node = std::get<1>(entry); Real distance = computeDistance(ref_node, eval_node); if (distance <= 0.95 * master_radiuses(j) + 1e-9) { ok = true; } } if (!ok) { AKANTU_EXCEPTION("Node " << eval_node << " is not in contact with any master node"); } } for (UInt eval_node : master_node_group) { bool ok = false; for (auto entry : enumerate(slave_node_group)) { auto j = std::get<0>(entry); auto ref_node = std::get<1>(entry); Real distance = computeDistance(ref_node, eval_node); if (distance <= 0.95 * slave_radiuses(j) + 1e-9) { ok = true; } } if (!ok) { AKANTU_EXCEPTION("Node " << eval_node << " is not in contact with any slave node"); } } } } /* -------------------------------------------------------------------------- */ SpatialGrid ContactDetectorInternodes::constructGrid(const akantu::NodeGroup & node_group, akantu::Real spacing) const { SpatialGrid grid(spatial_dimension); auto spacing_vector = Vector(spatial_dimension, spacing); grid.setSpacing(spacing_vector); for (UInt node : node_group) { grid.insert(node, getNodePosition(node)); } return grid; } /* -------------------------------------------------------------------------- */ Matrix ContactDetectorInternodes::constructInterpolationMatrix( const NodeGroup & ref_node_group, const NodeGroup & eval_node_group, Array eval_radiuses) { auto && phi_eval_eval = constructPhiMatrix(eval_node_group, eval_node_group, eval_radiuses); auto && phi_ref_eval = constructPhiMatrix(ref_node_group, eval_node_group, eval_radiuses); Matrix phi_eval_eval_inv(eval_node_group.size(), eval_node_group.size(), 1.); phi_eval_eval_inv.inverse(phi_eval_eval); auto && interpol_ref_eval = phi_ref_eval * phi_eval_eval_inv; for (UInt i : arange(ref_node_group.size())) { Real row_sum = 0.; for (UInt k : arange(eval_node_group.size())) { row_sum = row_sum + interpol_ref_eval(i, k); } for (UInt j : arange(eval_node_group.size())) { interpol_ref_eval(i, j) = interpol_ref_eval(i, j) / row_sum; } } // extended to 2D Matrix interpol_ref_eval_ext(spatial_dimension*interpol_ref_eval.rows(), spatial_dimension*interpol_ref_eval.cols(), 0.); for (UInt i : arange(interpol_ref_eval.rows())) { for (UInt j : arange(interpol_ref_eval.cols())) { for (int dim = 0; dim < spatial_dimension; dim++) { interpol_ref_eval_ext(spatial_dimension*i+dim, spatial_dimension*j+dim) = interpol_ref_eval(i, j); } } } return interpol_ref_eval_ext; } /* -------------------------------------------------------------------------- */ Matrix ContactDetectorInternodes::constructPhiMatrix(const NodeGroup & ref_node_group, const NodeGroup & eval_node_group, Array & eval_radiuses) { auto && positions_view = make_view(positions, spatial_dimension).begin(); Array distances(eval_node_group.size()); Matrix phi(ref_node_group.size(), eval_node_group.size()); phi.set(0.); for (auto && ref_node_data : enumerate(ref_node_group.getNodes())) { auto i = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); computeDistancesToRefNode(ref_node, eval_node_group, distances); for (UInt j : arange(eval_node_group.size())) { if (distances(j) <= eval_radiuses(j)) { phi(i, j) = computeRadialBasisInterpolation(distances(j), eval_radiuses(j)); } } } return phi; } /* -------------------------------------------------------------------------- */ std::set ContactDetectorInternodes::findPenetratingNodes( const NodeGroup & ref_group, const NodeGroup & eval_group, const Array & eval_radiuses) { const Real penetration_tolerance = -this->max_element_size * this->relative_penetration_tolerance; const auto R_ref_eval = constructInterpolationMatrix(ref_group, eval_group, eval_radiuses); std::set ref_penetration_nodes; for (auto && entry : enumerate(ref_group)) { auto i = std::get<0>(entry); auto ref_node = std::get<1>(entry); // 1) Compute gap: ref gaps = R_ref_eval * eval_positions - ref_positions // (pointing towards the inside of the ref body if penetrating) Vector gap(spatial_dimension); for (auto && inner_entry : enumerate(eval_group)) { auto j = std::get<0>(inner_entry); auto eval_node = std::get<1>(inner_entry); for (UInt s : arange(spatial_dimension)) { gap(s) += R_ref_eval(i * spatial_dimension + s, j * spatial_dimension + s) * positions(eval_node, s); } } gap -= positions(ref_node); // 2) Compute normal at ref_node by averaging element normals auto normal = getInterfaceNormalAtNode(ref_group, ref_node); // 3) Penetration if dot(gap, normal) is negative (with some tolerance) if (normal.dot(gap) < penetration_tolerance) { ref_penetration_nodes.insert(ref_node); } } return ref_penetration_nodes; } /* -------------------------------------------------------------------------- */ Vector ContactDetectorInternodes::getInterfaceNormalAtNode(const NodeGroup & interface_group, UInt node) const { if (spatial_dimension != 2) { AKANTU_EXCEPTION("Only spatial dimensions of 2 are supported at the moment"); } Vector tangent2 {{0, 0, 1}}; Vector normal(3); Array node_elements; mesh.getAssociatedElements(node, node_elements); for (auto && element : node_elements) { // We only consider boundary elements, i.e. elements that only contain nodes in the interface bool element_ok = true; int element_size = 0; auto && connectivity = mesh.getConnectivity(element); for (auto other_node : connectivity) { element_size++; if (interface_group.find(other_node) == -1) { element_ok = false; } } if (!element_ok || element_size != 2) { continue; } auto tangent1 = getNodePosition(connectivity(1)) - getNodePosition(connectivity(0)); // We need it in 3D for the cross product Vector tangent1_3d(3); for (UInt i = 0; i < 2; ++i) { tangent1_3d(i) = tangent1(i); } normal += tangent1_3d.crossProduct(tangent2); } // Normalize the normal normal.normalize(); return normal; } /* -------------------------------------------------------------------------- */ std::map ContactDetectorInternodes::computeRadiuses( Array & attack_radiuses, const NodeGroup & ref_node_group, const SpatialGrid & ref_grid, const NodeGroup & eval_node_group, const SpatialGrid & eval_grid) { Real c = 0.5; Real C = 0.95; std::vector temp_nodes; attack_radiuses.resize(ref_node_group.size()); std::map nb_neighboor_nodes_inside_radiuses; std::map nb_opposite_nodes_inside_radiuses; UInt nb_iter = 0; while (nb_iter < MAX_RADIUS_ITERATIONS) { // maximum number of support nodes UInt f = std::floor(1 / (std::pow(1 - c, 4) * (1 + 4 * c))); for (auto && ref_node_data : enumerate(ref_node_group.getNodes())) { auto j = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); // Compute radius of attack, i.e. distance to closest neighbor node with // a scaling factor of c. // This guarantees that [1], page 51, equation (2) is satisfied. Real attack_radius = std::numeric_limits::max(); for (auto neighboor_node : ref_grid.setToNeighboring(getNodePosition(ref_node), temp_nodes)) { if (neighboor_node != ref_node) { Real distance = computeDistance(ref_node, neighboor_node); attack_radius = std::min(distance / c, attack_radius); } } attack_radiuses(j) = attack_radius; // compute number of neighboor nodes inside attack radius for (auto neighbor_node : ref_grid.setToNeighboring(getNodePosition(ref_node), temp_nodes)) { Real distance = computeDistance(ref_node, neighbor_node); if (distance <= attack_radius) { nb_neighboor_nodes_inside_radiuses[neighbor_node]++; } } // compute number of opposite nodes inside C * attack_radius for (auto opposite_node : eval_grid.setToNeighboring(getNodePosition(ref_node), temp_nodes)) { Real distance = computeDistance(ref_node, opposite_node); if (distance < C * attack_radius) { nb_opposite_nodes_inside_radiuses[opposite_node]++; } } } // maximum number of neighboors inside radius // aka maximum number of supports UInt max_nb_supports = 0; for (const auto & entry : nb_neighboor_nodes_inside_radiuses) { max_nb_supports = std::max(max_nb_supports, entry.second); } // Enforce strict diagonal dominance by rows. // See [1], page 51, equation (4). if (max_nb_supports <= f) { // ok! break; } // correct maximum number of support nodes and then iterate again c = 0.5 * (1 + c); nb_neighboor_nodes_inside_radiuses.clear(); nb_opposite_nodes_inside_radiuses.clear(); nb_iter++; } if (nb_iter == MAX_RADIUS_ITERATIONS) { AKANTU_EXCEPTION("Could not find suitable radii, maximum number of iterations (" << nb_iter << ") was exceeded"); } if (DEBUG_VERIFY_INTERPOLATION_CONDITIONS) { for (UInt ref_node : ref_node_group) { for (auto entry : enumerate(ref_node_group)) { auto j = std::get<0>(entry); auto other_ref_node = std::get<1>(entry); if (ref_node == other_ref_node) { continue; } Real distance = computeDistance(ref_node, other_ref_node); if (distance < c * attack_radiuses(j) - 1e-9) { AKANTU_EXCEPTION("Radius of attack is too small, distance between nodes " << ref_node << " and " << other_ref_node << " is " << distance << " but radius is " << attack_radiuses(j)); } } } } return nb_opposite_nodes_inside_radiuses; } /* -------------------------------------------------------------------------- */ void ContactDetectorInternodes::computeDistancesToRefNode( UInt & ref_node, const NodeGroup & eval_node_group, Array & out_array) { for (auto && eval_node_data : enumerate(eval_node_group)) { auto i = std::get<0>(eval_node_data); auto eval_node = std::get<1>(eval_node_data); out_array(i) = computeDistance(ref_node, eval_node); } } /* -------------------------------------------------------------------------- */ Real ContactDetectorInternodes::computeRadialBasisInterpolation( const Real distance, const Real radius) { /// rescaled radial basis function: Wendland Real ratio = distance / radius; Real phi_of_x = std::pow(1 - ratio, 4) * (1 + 4 * ratio); return phi_of_x; } } // namespace akantu diff --git a/src/model/contact_mechanics_internodes/contact_detector_internodes.hh b/src/model/contact_mechanics_internodes/contact_detector_internodes.hh index 73feff6ae..06581aecc 100644 --- a/src/model/contact_mechanics_internodes/contact_detector_internodes.hh +++ b/src/model/contact_mechanics_internodes/contact_detector_internodes.hh @@ -1,165 +1,165 @@ /** * @file contact_detector_internodes.hh * * @author Moritz Waldleben * * @date creation: Thu Jul 09 2022 * @date last modification: Thu Jul 17 2022 * * @brief Algorithm to detetect contact nodes * * * @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_detector_shared.hh" +#include "contact_mechanics_shared.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_INTERNODES_HH__ #define __AKANTU_CONTACT_DETECTOR_INTERNODES_HH__ namespace akantu { class Mesh; class NodeGroup; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ContactDetectorInternodes : public Parsable, public AbstractContactDetector { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactDetectorInternodes(Mesh & mesh, const ID & id = "contact_detector"); ~ContactDetectorInternodes() override = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// find contact nodes for iteration void findContactNodes(NodeGroup & master_group, NodeGroup & slave_group); /// construct interpolation matrices Matrix constructInterpolationMatrix(const NodeGroup & ref_node_group, const NodeGroup & eval_node_group, Array eval_radiuses); /// construct Phi matrix used for interpolation Matrix constructPhiMatrix(const NodeGroup & ref_node_group, const NodeGroup & eval_node_group, Array & eval_radiuses); /// find nodes in the ref_group that are penetrating the eval_group std::set findPenetratingNodes(const NodeGroup & ref_group, const NodeGroup & eval_group, const Array & eval_radiuses); /// compute interface normal at a node Vector getInterfaceNormalAtNode(const NodeGroup & interface_group, UInt node) const; private: /// reads the input file to get contact detection options void parseSection(const ParserSection & section) override; /// construct a spatial grid with the given nodes and spacing SpatialGrid constructGrid(const NodeGroup & node_group, Real spacing) const; /// compute radius to detect contact nodes std::map computeRadiuses(Array & attack_radiuses, const NodeGroup & ref_node_group, const SpatialGrid & ref_grid, const NodeGroup & eval_node_group, const SpatialGrid & eval_grid); /// radial basis function Real computeRadialBasisInterpolation(Real distance, Real radius); /// distances between a reference node and other nodes /// the out_array must already be allocated with sufficient size void computeDistancesToRefNode(UInt & ref_node, const NodeGroup & eval_node_group, Array & out_array); inline Real computeDistance(UInt node1, UInt node2) const { return getNodePosition(node1).distance(getNodePosition(node2)); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get radiuses of attack of master nodes AKANTU_GET_MACRO_NOT_CONST(MasterRadiuses, master_radiuses, Array &) /// get radiuses of attack of master nodes AKANTU_GET_MACRO_NOT_CONST(SlaveRadiuses, slave_radiuses, Array &) /// get initial master node group NodeGroup & getInitialMasterNodeGroup(); /// get initial slave node group NodeGroup & getInitialSlaveNodeGroup(); /// get master node group NodeGroup & getMasterNodeGroup(); /// get slave node group NodeGroup & getSlaveNodeGroup(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of master node group, i.e nodes on interface ID id_master_nodes{}; /// id of slave nodes group, i.e nodes on interface ID id_slave_nodes{}; /// relative security factory for detection of close nodes Real relative_grid_spacing_factor; /// relative tolerance for detection of penetration Real relative_penetration_tolerance; /// attack radiuses for master nodes Array master_radiuses{0}; /// attack radiuses for master nodes Array slave_radiuses{0}; /// maximum element inradius Real max_element_size; /// maximum number of iterations for radius computation: /// after each iteration, c is replaced by (1+c)/2, so if this value is ever /// reached it means that c is close to 1 and we won't find suitable radii. const UInt MAX_RADIUS_ITERATIONS = 10; /// true to validate the interpolation conditions /// the method should always produce node groups satisfying the conditions /// should only be used to debug the method const bool DEBUG_VERIFY_INTERPOLATION_CONDITIONS = false; }; } // namespace akantu #endif /* __AKANTU_CONTACT_DETECTOR_INTERNODES_HH__ */ diff --git a/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.cc b/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.cc index 5e2fd8ed2..868839c29 100644 --- a/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.cc +++ b/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.cc @@ -1,548 +1,488 @@ /** * @file contact_mechanics_internodes_model.cc * * @author Moritz Waldleben * * @date creation: Thu Jul 09 2022 * @date last modification: Thu Jul 17 2022 * * @brief Model for Contact Mechanics Internodes * * * @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_internodes_model.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsInternodesModel::ContactMechanicsInternodesModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id) { + : AbstractContactMechanicsModel(mesh, model_type, dof_manager, dim, id) { this->detector = std::make_unique( this->mesh, id + ":contact_detector"); - - this->solid = std::make_unique( - this->mesh, spatial_dimension, id + ":solid_mechanics_model", - this->dof_manager); - - setYoungsModulus(this->solid->getMaterial(0).get("E")); } /* -------------------------------------------------------------------------- */ ContactMechanicsInternodesModel::~ContactMechanicsInternodesModel() = default; /* -------------------------------------------------------------------------- */ void ContactMechanicsInternodesModel::initFullImpl( const ModelOptions & options) { // check if static solver if (options.analysis_method != _static) { AKANTU_EXCEPTION(options.analysis_method << " is not a valid analysis method for " "ContactMechanicsInternodesModel"); } - // run contact detection - detector->findContactNodes(detector->getMasterNodeGroup(), detector->getSlaveNodeGroup()); - Model::initFullImpl(options); - - // init solid mechanics model - solid->initFull(options); } /* -------------------------------------------------------------------------- */ void ContactMechanicsInternodesModel::assembleResidual() { - auto & solid_model_solver = aka::as_type(*solid); - solid_model_solver.assembleResidual(); - auto & master_node_group = detector->getMasterNodeGroup(); auto & initial_master_node_group = detector->getInitialMasterNodeGroup(); auto & slave_node_group = detector->getSlaveNodeGroup(); Vector positions_master(master_node_group.size()*spatial_dimension); Vector positions_slave(slave_node_group.size()*spatial_dimension); for (auto && node_data : enumerate(master_node_group.getNodes())) { auto i = std::get<0>(node_data); auto node = std::get<1>(node_data); auto pos = detector->getNodePosition(node); for (int dim = 0; dim < spatial_dimension; dim++) { positions_master(spatial_dimension*i + dim) = pos(dim); } } for (auto && node_data : enumerate(slave_node_group.getNodes())) { auto i = std::get<0>(node_data); auto node = std::get<1>(node_data); auto pos = detector->getNodePosition(node); for (int dim = 0; dim < spatial_dimension; dim++) { positions_slave(spatial_dimension*i + dim) = pos(dim); } } // interpolation matrix auto && R_master_slave_ext = detector->constructInterpolationMatrix( master_node_group, slave_node_group, detector->getSlaveRadiuses()); // calculate differences Array differences(initial_master_node_group.size(), spatial_dimension, 0.); // R_master_slave_ext * positions_slave - positions_master Vector tmp = (R_master_slave_ext * positions_slave - positions_master) * this->E; UInt next_active_node = 0; for (auto && node_data : enumerate(initial_master_node_group.getNodes())) { auto i = std::get<0>(node_data); auto node = std::get<1>(node_data); if (master_node_group.find(node) != -1) { UInt active_node_index = next_active_node++; for (int dim = 0; dim < spatial_dimension; dim++) { differences(i, dim) = tmp(spatial_dimension*active_node_index + dim); } } } this->dof_manager->assembleToResidual("lambdas", differences, 1); } /* -------------------------------------------------------------------------- */ MatrixType ContactMechanicsInternodesModel::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K") { return _unsymmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsInternodesModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { assembleInternodesMatrix(); } - - if (matrix_id == "M") { - solid->assembleMass(); - } } /* -------------------------------------------------------------------------- */ void ContactMechanicsInternodesModel::assembleLumpedMatrix( const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsInternodesModel::assembleInternodesMatrix() { - solid->assembleStiffnessMatrix(); - auto & initial_master_node_group = detector->getInitialMasterNodeGroup(); - auto & initial_slave_node_group = detector->getInitialSlaveNodeGroup(); auto & master_node_group = detector->getMasterNodeGroup(); auto & slave_node_group = detector->getSlaveNodeGroup(); // R matrices, need for calculations of blocks auto && R_master_slave = detector->constructInterpolationMatrix( master_node_group, slave_node_group, detector->getSlaveRadiuses()); auto && R_slave_master = detector->constructInterpolationMatrix( slave_node_group, master_node_group, detector->getMasterRadiuses()); // get interface masses auto && M_master = assembleInterfaceMass(master_node_group); auto && M_slave = assembleInterfaceMass(slave_node_group); // assemble B and B_tilde // B_master = - M_master; auto B_slave = M_slave * R_slave_master; //B_tilde_master = I; //B_tilde_slave = - R_master_slave; TermsToAssemble termsB("displacement", "lambdas"); TermsToAssemble termsB_tilde("lambdas", "displacement"); // Only used with a "diagonal 1" for "inactive" nodes, to make the matrix non-singular TermsToAssemble termsLambdaDiag("lambdas", "lambdas"); // fill in active nodes for (auto && ref_node_data : enumerate(initial_master_node_group.getNodes())) { auto i = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); auto local_i = master_node_group.find(ref_node); if (local_i != -1) { for (auto && master_node_data : enumerate(master_node_group.getNodes())) { auto j = std::get<0>(master_node_data); auto master_node = std::get<1>(master_node_data); for (int dim = 0; dim < spatial_dimension; dim++) { UInt idx_i = local_i*spatial_dimension + dim; UInt idx_j = j*spatial_dimension + dim; UInt global_idx_i = i*spatial_dimension + dim; UInt global_idx_j = master_node*spatial_dimension + dim; termsB(global_idx_j, global_idx_i) = - this->E*M_master(idx_j, idx_i); if (idx_i == idx_j) { termsB_tilde(global_idx_i, global_idx_j) = this->E; } } } } else { for (int s = 0; s < spatial_dimension; s++) { UInt global_idx = i * spatial_dimension + s; termsLambdaDiag(global_idx, global_idx) = this->E * this->E; } } } for (auto && ref_node_data : enumerate(initial_master_node_group.getNodes())) { auto i = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); auto local_i = master_node_group.find(ref_node); if (local_i != UInt(-1)) { for (auto && slave_node_data : enumerate(slave_node_group.getNodes())) { auto j = std::get<0>(slave_node_data); auto slave_node = std::get<1>(slave_node_data); for (int dim = 0; dim < spatial_dimension; dim++) { UInt idx_i = local_i*spatial_dimension + dim; UInt idx_j = j*spatial_dimension + dim; UInt global_idx_i = i*spatial_dimension + dim; UInt global_idx_j = slave_node*spatial_dimension + dim; termsB(global_idx_j, global_idx_i) = this->E * B_slave(idx_j, idx_i); termsB_tilde(global_idx_i, global_idx_j) = - this->E * R_master_slave(idx_i, idx_j); } } } else { for (int s = 0; s < spatial_dimension; s++) { UInt global_idx = i * spatial_dimension + s; termsLambdaDiag(global_idx, global_idx) = this->E * this->E; } } } this->dof_manager->assemblePreassembledMatrix("K", termsB); this->dof_manager->assemblePreassembledMatrix("K", termsB_tilde); this->dof_manager->assemblePreassembledMatrix("K", termsLambdaDiag); } /* -------------------------------------------------------------------------- */ Matrix ContactMechanicsInternodesModel::assembleInterfaceMass( const NodeGroup & contact_node_group) { // get interface mass for a node group this->assembleMatrix("M"); const auto & M = this->dof_manager->getMatrix("M"); auto nb_contact_nodes = contact_node_group.size(); Matrix M_contact(nb_contact_nodes*spatial_dimension, nb_contact_nodes*spatial_dimension); for (auto && ref_node_data : enumerate(contact_node_group.getNodes())) { auto i = std::get<0>(ref_node_data); auto ref_node = std::get<1>(ref_node_data); for (auto && node_data : enumerate(contact_node_group.getNodes())) { auto j = std::get<0>(node_data); auto node = std::get<1>(node_data); for (int dim = 0; dim < spatial_dimension; dim++) { UInt global_idx_ref = ref_node*spatial_dimension + dim; UInt global_idx = node*spatial_dimension + dim; UInt idx_ref = i*spatial_dimension + dim; UInt idx = j*spatial_dimension + dim; // M must be constant otherwise the modifying overload of operator() is // used, which throws if the entry is absent M_contact(idx_ref, idx) = M(global_idx_ref, global_idx); } } } return M_contact; } /* -------------------------------------------------------------------------- */ bool ContactMechanicsInternodesModel::updateAfterStep() { - // Update positions - detector->getPositions().copy(solid->getCurrentPosition()); - auto & initial_master_node_group = detector->getInitialMasterNodeGroup(); auto & master_node_group = detector->getMasterNodeGroup(); auto & slave_node_group = detector->getSlaveNodeGroup(); // Find nodes to remove (if their projected Lagrange multiplier is positive) std::set to_remove_master, to_remove_slave; { // Master // Note: the lambdas are actually lambdas/E, but that's fine for our sign check const auto & lambda_solution = dof_manager->getSolution("lambdas"); Vector lambdas_master(master_node_group.size() * spatial_dimension); // We need to filter out unused lambdas UInt next_used_lambda = 0; for (auto && entry : enumerate(initial_master_node_group)) { UInt i = std::get<0>(entry); UInt maybe_active_node = std::get<1>(entry); if (master_node_group.find(maybe_active_node) != -1) { for (UInt s = 0; s < spatial_dimension; ++s) { lambdas_master(next_used_lambda++) = lambda_solution(i * spatial_dimension + s); } } } for (auto && entry : enumerate(master_node_group)) { auto i = std::get<0>(entry); auto node = std::get<1>(entry); auto normal = detector->getInterfaceNormalAtNode(master_node_group, node); Real dot_product = 0; for (UInt s : arange(spatial_dimension)) { dot_product += lambdas_master(i * spatial_dimension + s) * normal(s); } if (dot_product > 1e-9) { to_remove_master.insert(node); } } // Slave // Interpolate the Lagrange multipliers of the slave auto old_R21 = detector->constructInterpolationMatrix( slave_node_group, master_node_group, detector->getMasterRadiuses()); auto lambdas_slave = old_R21 * lambdas_master; lambdas_slave *= -1; for (auto && entry : enumerate(slave_node_group)) { auto i = std::get<0>(entry); auto node = std::get<1>(entry); auto normal = detector->getInterfaceNormalAtNode(slave_node_group, node); Real dot_product = 0; for (UInt s : arange(spatial_dimension)) { dot_product += lambdas_slave(i * spatial_dimension + s) * normal(s); } if (dot_product > 1e-9) { to_remove_slave.insert(node); } } } // Find nodes to add (if they are penetrating) std::set to_add_master, to_add_slave; { // COMPUTE PENETRATION // Start with all the initial nodes NodeGroup penetration_master_group("penetration_master", mesh); NodeGroup penetration_slave_group("penetration_slave", mesh); penetration_master_group.append(detector->getInitialMasterNodeGroup()); penetration_slave_group.append(detector->getInitialSlaveNodeGroup()); // Find contact nodes and radii with the contact detection algorithm. // This "pollutes" the state of the detector, but it doesn't matter because we won't be using it again in this iteration. detector->findContactNodes(penetration_master_group, penetration_slave_group); // Compute penetrating nodes (i.e. nodes to add to the interface) to_add_master = detector->findPenetratingNodes(penetration_master_group, penetration_slave_group, detector->getSlaveRadiuses()); to_add_slave = detector->findPenetratingNodes(penetration_slave_group, penetration_master_group, detector->getMasterRadiuses()); } // Remove nodes with positive projected Lagrange multipliers bool interface_nodes_changed = false; if (master_node_group.applyNodeFilter([&] (UInt master_node) { bool keep = to_remove_master.count(master_node) == 0; return keep; }) > 0) { interface_nodes_changed = true; } if (slave_node_group.applyNodeFilter([&] (UInt slave_node) { bool keep = to_remove_slave.count(slave_node) == 0; return keep; }) > 0) { interface_nodes_changed = true; } // Add nodes that are penetrating for (auto node : master_node_group) { to_add_master.erase(node); } for (auto node : slave_node_group) { to_add_slave.erase(node); } for (auto node : to_add_master) { master_node_group.add(node); interface_nodes_changed = true; } for (auto node : to_add_slave) { slave_node_group.add(node); interface_nodes_changed = true; } if (interface_nodes_changed) { // Make sure to re-optimize the node groups after a modification! master_node_group.optimize(); slave_node_group.optimize(); } return interface_nodes_changed; } /* -------------------------------------------------------------------------- */ -void ContactMechanicsInternodesModel::solveStep(SolverCallback & callback, - const ID & solver_id) { - ModelSolver::solveStep(callback, solver_id); -} - -/* -------------------------------------------------------------------------- */ -void ContactMechanicsInternodesModel::solveStep(const ID & solver_id) { - - ModelSolver::solveStep(solver_id); -} - -/* -------------------------------------------------------------------------- */ -void ContactMechanicsInternodesModel::predictor() { - auto & solid_callback = aka::as_type(*solid); - solid_callback.predictor(); -} - -/* -------------------------------------------------------------------------- */ -void ContactMechanicsInternodesModel::corrector() { - auto & solid_callback = aka::as_type(*solid); - solid_callback.corrector(); -} - -/* -------------------------------------------------------------------------- */ -void ContactMechanicsInternodesModel::beforeSolveStep() { - auto & solid_callback = aka::as_type(*solid); - solid_callback.beforeSolveStep(); -} - -/* -------------------------------------------------------------------------- */ -void ContactMechanicsInternodesModel::afterSolveStep(bool converged) { - auto & solid_callback = aka::as_type(*solid); - solid_callback.afterSolveStep(converged); +void ContactMechanicsInternodesModel::search() { + detector->findContactNodes(detector->getMasterNodeGroup(), detector->getSlaveNodeGroup()); } void ContactMechanicsInternodesModel::printself(std::ostream & stream, int indent) const { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ TimeStepSolverType ContactMechanicsInternodesModel::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsInternodesModel::getDefaultSolverID( const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsInternodesModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type for " "ContactMechanicsInternodesModel"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsInternodesModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { - - auto & solid_model_solver = aka::as_type(*solid); - solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); - auto nb_initial_master_nodes = detector->getInitialMasterNodeGroup().size(); // allocate lambdas this->lambdas = std::make_unique>(nb_initial_master_nodes, spatial_dimension, id + ":lambdas"); // Note that we always allocate the lambdas for all the initial master nodes. // If we decide to remove some initial master nodes from the interface, we // will fill the system with an "identity matrix" for the removed nodes. // allocate blocked dofs for lambdas this->blocked_dofs = std::make_unique>( nb_initial_master_nodes, spatial_dimension, false, id + ":blocked_dofs"); this->dof_manager->registerDOFs("lambdas", *lambdas, _dst_generic); this->dof_manager->registerBlockedDOFs("lambdas", *blocked_dofs); } } // namespace akantu diff --git a/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.hh b/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.hh index 9307e75fa..f34120bd6 100644 --- a/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.hh +++ b/src/model/contact_mechanics_internodes/contact_mechanics_internodes_model.hh @@ -1,187 +1,174 @@ /** * @file contact_mechanics_internodes_model.cc * * @author Moritz Waldleben * * @date creation: Thu Jul 09 2022 * @date last modification: Thu Jul 17 2022 * * @brief Model for Contact Mechanics Internodes * * * @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_detector_internodes.hh" #include "solid_mechanics_model.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_INTERNODES_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_INTERNODES_MODEL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { /// Implementation of a contact detector using the internodes method. ///

/// The pseudo-code for the internodes method is as follows: ///

    ///
  1. Find contact nodes.
    2. ///
  2. Assemble and solve internodes system.
    3. ///
  3. If all the Lagrange multipliers are positive and all the penetrating /// nodes are part of the interface, we are finished.
    4. ///
  4. Otherwise, remove nodes with positive Lagrange multipliers, add /// missing penetrating nodes, and go back to 1.
    5. ///
///

///

References

///
    ///
  1. Y. Voet et. al.: The INTERNODES method for applications in contact /// mechanics and dedicated preconditioning techniques. ///

    Computers & /// Mathematics with Applications, vol. 127, 2022, pp. 48-64.

    2. /// ///
  2. Y. Voet: On the preconditioning of the INTERNODES matrix for /// applications in contact mechanics. ///

    Master's thesis EPFL, 2021.

    3. ///
-class ContactMechanicsInternodesModel : public Model { +class ContactMechanicsInternodesModel : public AbstractContactMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactMechanicsInternodesModel( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "contact_mechanics_internodes_model", std::shared_ptr dof_manager = nullptr, + // TODO: fix model type ModelType model_type = ModelType::_solid_mechanics_model); ~ContactMechanicsInternodesModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: - /// costum solve step for Internodes - void solveStep(SolverCallback & callback, const ID & solver_id = "") override; - void solveStep(const ID & solver_id = "") override; + /// Perform contact detector search + void search() override; // assemble extended matrix K void assembleInternodesMatrix(); /// assemble an interface matrix, eather master or slave Matrix assembleInterfaceMass(const NodeGroup & contact_node_group); /// update nodes in the interface after a solve step, return true if the /// contact nodes have changed (i.e. we need more iteration) bool updateAfterStep(); protected: // call back for the solver, computes the force residual void assembleResidual() override; // 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; // get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) 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; - /// allocate all vectors void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: TimeStepSolverType getDefaultSolverType() const override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get contact detector - AKANTU_GET_MACRO(ContactDetectorInternodes, *detector, - ContactDetectorInternodes &); - - /// get solid mechanics model - AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &); + AKANTU_GET_MACRO(ContactDetector, *detector, ContactDetectorInternodes &); /// get lambdas from internodes formulation AKANTU_GET_MACRO(Lambdas, *lambdas, Array); /// set representative Young's modulus AKANTU_SET_MACRO(YoungsModulus, E, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// contact detector std::unique_ptr detector; - /// solid mechanics model - std::unique_ptr solid; - /// lambdas array std::unique_ptr> lambdas; /// blocked dofs for lambdas std::unique_ptr> blocked_dofs; /// Young's modulus of one of the materials /// used as a scaling factor for lambdas, to lower condition number Real E; + + template + friend class AbstractCouplerSolidContactTemplate; + template + friend class CouplerSolidContactInternodesTemplate; }; } // namespace akantu #endif /* __AKANTU_CONTACT_MECHANICS_INTERNODES_MODEL_HH__ */ diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc index 3a2099ad0..ef5bd8d42 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,76 +1,79 @@ /** * @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:: +CouplerSolidCohesiveContact:: CouplerSolidContactTemplate(Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) - : Model(mesh, ModelType::_coupler_solid_cohesive_contact, dof_manager, dim, - id) { + : AbstractCouplerSolidContactTemplate( + mesh, ModelType::_coupler_solid_cohesive_contact, dim, id, + dof_manager) { 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( +void AbstractCouplerSolidContactTemplate::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 5f111e9e2..8efe4863d 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,69 +1,71 @@ /** * @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) - : Model(mesh, ModelType::_coupler_solid_contact, dof_manager, dim, id) { + : AbstractCouplerSolidContactTemplate( + mesh, ModelType::_coupler_solid_contact, dim, id, dof_manager) { 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) { +void AbstractCouplerSolidContactTemplate + ::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 6c04b4770..feb6d8fd2 100644 --- a/src/model/model_couplers/coupler_solid_contact.hh +++ b/src/model/model_couplers/coupler_solid_contact.hh @@ -1,276 +1,287 @@ /** * @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, +template +class AbstractCouplerSolidContactTemplate : public Model, public DataAccessor, public DataAccessor { static_assert( std::is_base_of::value, "SolidMechanicsModelType should be derived from SolidMechanicsModel"); + static_assert( + std::is_base_of::value, + "ContactMechanicsModelType should be derived from AbstractContactMechanicsModel"); /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ -public: - CouplerSolidContactTemplate( - Mesh & mesh, UInt dim = _all_dimensions, - const ID & id = "coupler_solid_contact", +protected: + AbstractCouplerSolidContactTemplate( + Mesh & mesh, const ModelType & type, UInt dim, const ID & id, std::shared_ptr dof_manager = nullptr); - ~CouplerSolidContactTemplate() override; + ~AbstractCouplerSolidContactTemplate() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ + /// solid mechanics model + std::unique_ptr solid; + + /// contact mechanics model + std::unique_ptr contact; + 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; } + void assembleResidual() override = 0; /// get the type of matrix needed - MatrixType getMatrixType(const ID & matrix_id) const override; + MatrixType getMatrixType(const ID & matrix_id) const override = 0; /// callback for the solver, this assembles different matrices - void assembleMatrix(const ID & matrix_id) override; + void assembleMatrix(const ID & matrix_id) override = 0; /// 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 &) + AKANTU_GET_MACRO(ContactMechanicsModel, *contact, ContactMechanicsModelType &) /* ------------------------------------------------------------------------ */ /* 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; + UInt step; +}; - /// contact mechanics model - std::unique_ptr contact; +template +class CouplerSolidContactTemplate : public AbstractCouplerSolidContactTemplate< + SolidMechanicsModelType, ContactMechanicsModel> { +public: + CouplerSolidContactTemplate( + Mesh & mesh, UInt dim = _all_dimensions, + const ID & id = "coupler_solid_contact", + std::shared_ptr dof_manager = nullptr); - UInt step; + ~CouplerSolidContactTemplate() override = default; + +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 is called at end of solve + void corrector() override; + +private: + /// assembles the contant internal forces + void assembleInternalForces(); + + /// assembles the contact stiffness matrix + void assembleStiffnessMatrix(); }; using CouplerSolidContact = CouplerSolidContactTemplate; } // namespace akantu #include "coupler_solid_contact_tmpl.hh" #endif /* __COUPLER_SOLID_CONTACT_HH__ */ diff --git a/src/model/model_couplers/coupler_solid_contact_internodes.cc b/src/model/model_couplers/coupler_solid_contact_internodes.cc new file mode 100644 index 000000000..bd7cbe2bc --- /dev/null +++ b/src/model/model_couplers/coupler_solid_contact_internodes.cc @@ -0,0 +1,38 @@ +#include "coupler_solid_contact_internodes.hh" + +namespace akantu { + +template <> +CouplerSolidContactInternodesTemplate::CouplerSolidContactInternodesTemplate( + Mesh & mesh, UInt dim, const ID & id, + std::shared_ptr dof_manager) + : AbstractCouplerSolidContactTemplate( + mesh, ModelType::_coupler_solid_contact, dim, id, dof_manager) { + 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_internodes_model", + this->dof_manager); +} + +/* -------------------------------------------------------------------------- */ +template <> +void AbstractCouplerSolidContactTemplate + ::initFullImpl(const ModelOptions & options) { + Model::initFullImpl(options); + + solid->initFull(_analysis_method = this->method); + contact->initFull(_analysis_method = this->method); + + contact->setYoungsModulus(solid->getMaterial(0).get("E")); +} + +} // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact_internodes.hh b/src/model/model_couplers/coupler_solid_contact_internodes.hh new file mode 100644 index 000000000..8d9c4f97c --- /dev/null +++ b/src/model/model_couplers/coupler_solid_contact_internodes.hh @@ -0,0 +1,50 @@ +/* -------------------------------------------------------------------------- */ +#include "coupler_solid_contact.hh" +#include "contact_mechanics_internodes_model.hh" + +/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_COUPLER_SOLID_CONTACT_INTERNODES_HH__ +#define __AKANTU_COUPLER_SOLID_CONTACT_INTERNODES_HH__ + +/* -------------------------------------------------------------------------- */ +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template +class CouplerSolidContactInternodesTemplate : + public AbstractCouplerSolidContactTemplate { +public: + CouplerSolidContactInternodesTemplate( + Mesh & mesh, UInt dim = _all_dimensions, + const ID & id = "coupler_solid_contact_internodes", + std::shared_ptr dof_manager = nullptr); + + ~CouplerSolidContactInternodesTemplate() override = default; + + /// custom solve step to handle the INTERNODES iterations + void solveStep(const ID & solver_id = "") override; + +protected: + /// callback for the solver to assemble the rhs + void assembleResidual() override; + + /// 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; + + ModelSolverOptions + getDefaultSolverOptions(const TimeStepSolverType & type) const override; +}; + +/* -------------------------------------------------------------------------- */ +using CouplerSolidContactInternodes = + CouplerSolidContactInternodesTemplate; + +} // namespace akantu + +/* -------------------------------------------------------------------------- */ +#include "coupler_solid_contact_internodes_tmpl.hh" + +#endif // __AKANTU_COUPLER_SOLID_CONTACT_INTERNODES_HH__ diff --git a/src/model/model_couplers/coupler_solid_contact_internodes_tmpl.hh b/src/model/model_couplers/coupler_solid_contact_internodes_tmpl.hh new file mode 100644 index 000000000..1144e5a72 --- /dev/null +++ b/src/model/model_couplers/coupler_solid_contact_internodes_tmpl.hh @@ -0,0 +1,90 @@ +/* -------------------------------------------------------------------------- */ +#include "coupler_solid_contact_internodes.hh" + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactInternodesTemplate + ::solveStep(const ID & solver_id) { + const UInt MAX_STEPS = 100; // security value for now, increase later if needed + + UInt step = 0; + while (true) { + // Always search before every step + this->contact->search(); + + // Perform regular solve + Model::solveStep(solver_id); + + if (++step >= MAX_STEPS) { + AKANTU_EXCEPTION("Reached maximum number of internodes steps: " << step); + } + + // Perform any required update, and solve again if needed + auto & current_positions = this->contact->getContactDetector().getPositions(); + current_positions.copy(this->solid->getCurrentPosition()); + if (!this->contact->updateAfterStep()) { + break; + } + } +} + +/* -------------------------------------------------------------------------- */ +template +MatrixType CouplerSolidContactInternodesTemplate + ::getMatrixType(const ID & matrix_id) const { + if (matrix_id == "K") { + return _unsymmetric; + } + // TODO: this wasn't in the original internodes model, but it should be symm? +// if (matrix_id == "M") { +// return _symmetric; +// } + return _mt_not_defined; +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactInternodesTemplate + ::assembleMatrix(const ID & matrix_id) { + if (matrix_id == "K") { + // We need the solid's mass matrix to compute the stiffness matrix. + this->solid->assembleMass(); + this->contact->assembleInternodesMatrix(); + } else if (matrix_id == "M") { + this->solid->assembleMass(); + } +} + +/* -------------------------------------------------------------------------- */ +template +ModelSolverOptions +CouplerSolidContactInternodesTemplate + ::getDefaultSolverOptions(const TimeStepSolverType & type) const { + return this->contact->getDefaultSolverOptions(type); +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactInternodesTemplate + ::assembleResidual() { +// switch (this->method) { + // TODO: this is needed because contact search hasn't been performed before, + // but why does the ContactMechanicsModel coupler not do it then? (_static) +// case _static: +// case _explicit_lumped_mass: { +// auto & current_positions = this->contact->getContactDetector().getPositions(); +// current_positions.copy(this->solid->getCurrentPosition()); +// this->contact->search(); +// break; +// } +// default: +// break; +// } + + this->solid->assembleResidual(); + this->contact->assembleResidual(); +} + +} // namespace akantu \ No newline at end of file diff --git a/src/model/model_couplers/coupler_solid_contact_tmpl.hh b/src/model/model_couplers/coupler_solid_contact_tmpl.hh index 5bfe1a59f..8f772bb39 100644 --- a/src/model/model_couplers/coupler_solid_contact_tmpl.hh +++ b/src/model/model_couplers/coupler_solid_contact_tmpl.hh @@ -1,396 +1,397 @@ /** * @file coupler_solid_contact_tmpl.hh * * @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" #include "dumpable_inline_impl.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template -CouplerSolidContactTemplate< - SolidMechanicsModelType>::~CouplerSolidContactTemplate() = default; +template +AbstractCouplerSolidContactTemplate + ::AbstractCouplerSolidContactTemplate(Mesh & mesh, const ModelType & type, + UInt dim, const ID & id, + std::shared_ptr dof_manager) + : Model(mesh, type, dof_manager, dim, id) { +} /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::initSolver( +template +void AbstractCouplerSolidContactTemplate + ::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); auto & contact_model_solver = aka::as_type(*contact); contact_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); } /* -------------------------------------------------------------------------- */ -template +template std::tuple -CouplerSolidContactTemplate::getDefaultSolverID( +AbstractCouplerSolidContactTemplate + ::getDefaultSolverID( const AnalysisMethod & method) { return solid->getDefaultSolverID(method); } /* -------------------------------------------------------------------------- */ -template +template TimeStepSolverType -CouplerSolidContactTemplate::getDefaultSolverType() +AbstractCouplerSolidContactTemplate + ::getDefaultSolverType() const { return solid->getDefaultSolverType(); } /* -------------------------------------------------------------------------- */ -template +template ModelSolverOptions -CouplerSolidContactTemplate::getDefaultSolverOptions( +AbstractCouplerSolidContactTemplate + ::getDefaultSolverOptions( const TimeStepSolverType & type) const { return solid->getDefaultSolverOptions(type); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleResidual() { - // computes the internal forces - switch (method) { - case _explicit_lumped_mass: { - auto & current_positions = contact->getContactDetector().getPositions(); - current_positions.copy(solid->getCurrentPosition()); - contact->search(); - break; - } - default: - break; - } - - this->assembleInternalForces(); - - auto & internal_force = solid->getInternalForce(); - auto & external_force = solid->getExternalForce(); - - auto & contact_force = contact->getInternalForce(); - - /* ------------------------------------------------------------------------ */ - this->getDOFManager().assembleToResidual("displacement", external_force, 1); - this->getDOFManager().assembleToResidual("displacement", internal_force, 1); - this->getDOFManager().assembleToResidual("displacement", contact_force, 1); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleResidual( - const ID & residual_part) { - AKANTU_DEBUG_IN(); - - // contact->assembleInternalForces(); - - auto & internal_force = solid->getInternalForce(); - auto & external_force = solid->getExternalForce(); - auto & contact_force = contact->getInternalForce(); - - if ("external" == residual_part) { - this->getDOFManager().assembleToResidual("displacement", external_force, 1); - this->getDOFManager().assembleToResidual("displacement", contact_force, 1); - AKANTU_DEBUG_OUT(); - return; - } - - if ("internal" == residual_part) { - this->getDOFManager().assembleToResidual("displacement", internal_force, 1); - AKANTU_DEBUG_OUT(); - return; - } - - AKANTU_CUSTOM_EXCEPTION( - debug::SolverCallbackResidualPartUnknown(residual_part)); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::predictor() { +template +void AbstractCouplerSolidContactTemplate + ::predictor() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::corrector() { +template +void AbstractCouplerSolidContactTemplate + ::corrector() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.corrector(); - - switch (method) { - case _static: - case _implicit_dynamic: { - auto & current_positions = contact->getContactDetector().getPositions(); - current_positions.copy(solid->getCurrentPosition()); - contact->search(); - break; - } - default: - break; - } } /* -------------------------------------------------------------------------- */ -template -MatrixType CouplerSolidContactTemplate::getMatrixType( - const ID & matrix_id) const { - if (matrix_id == "K") { - return _symmetric; - } - if (matrix_id == "M") { - return _symmetric; - } - return _mt_not_defined; -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleMatrix( - const ID & matrix_id) { - if (matrix_id == "K") { - this->assembleStiffnessMatrix(); - } else if (matrix_id == "M") { - solid->assembleMass(); - } -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleLumpedMatrix( +template +void AbstractCouplerSolidContactTemplate + ::assembleLumpedMatrix( const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::beforeSolveStep() { +template +void AbstractCouplerSolidContactTemplate + ::beforeSolveStep() { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.beforeSolveStep(); auto & contact_solver_callback = aka::as_type(*contact); contact_solver_callback.beforeSolveStep(); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::afterSolveStep( - bool converged) { +template +void AbstractCouplerSolidContactTemplate + ::afterSolveStep(bool converged) { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.afterSolveStep(converged); auto & contact_solver_callback = aka::as_type(*contact); contact_solver_callback.afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate< - SolidMechanicsModelType>::assembleInternalForces() { - AKANTU_DEBUG_IN(); - - AKANTU_DEBUG_INFO("Assemble the internal forces"); - - solid->assembleInternalForces(); - contact->assembleInternalForces(); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate< - SolidMechanicsModelType>::assembleStiffnessMatrix() { - AKANTU_DEBUG_IN(); - - AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); - - solid->assembleStiffnessMatrix(true); - - switch (method) { - case _static: - case _implicit_dynamic: { - contact->assembleStiffnessMatrix(); - break; - } - default: - break; - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate< - SolidMechanicsModelType>::assembleMassLumped() { - solid->assembleMassLumped(); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleMass() { - solid->assembleMass(); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleMassLumped( - GhostType ghost_type) { - solid->assembleMassLumped(ghost_type); -} - -/* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::assembleMass( - GhostType ghost_type) { - solid->assembleMass(ghost_type); -} - -/* -------------------------------------------------------------------------- */ -template +template std::shared_ptr -CouplerSolidContactTemplate::createElementalField( +AbstractCouplerSolidContactTemplate + ::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; field = contact->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); if (not field) { field = solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); } return field; } /* -------------------------------------------------------------------------- */ -template +template std::shared_ptr -CouplerSolidContactTemplate::createNodalFieldReal( +AbstractCouplerSolidContactTemplate + ::createNodalFieldReal( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = contact->createNodalFieldReal(field_name, group_name, padding_flag); if (not field) { field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ -template +template std::shared_ptr -CouplerSolidContactTemplate::createNodalFieldUInt( +AbstractCouplerSolidContactTemplate:: + createNodalFieldUInt( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = contact->createNodalFieldUInt(field_name, group_name, padding_flag); if (not field) { field = solid->createNodalFieldUInt(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ -template +template std::shared_ptr -CouplerSolidContactTemplate::createNodalFieldBool( +AbstractCouplerSolidContactTemplate + ::createNodalFieldBool( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = contact->createNodalFieldBool(field_name, group_name, padding_flag); if (not field) { field = solid->createNodalFieldBool(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::dump( - const std::string & dumper_name) { +template +void AbstractCouplerSolidContactTemplate + ::dump(const std::string & dumper_name) { solid->onDump(); Model::dump(dumper_name); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::dump( - const std::string & dumper_name, UInt step) { +template +void AbstractCouplerSolidContactTemplate + ::dump(const std::string & dumper_name, UInt step) { solid->onDump(); Model::dump(dumper_name, step); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::dump( - const std::string & dumper_name, Real time, UInt step) { +template +void AbstractCouplerSolidContactTemplate + ::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); Model::dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::dump() { +template +void AbstractCouplerSolidContactTemplate + ::dump() { solid->onDump(); Model::dump(); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::dump(UInt step) { +template +void AbstractCouplerSolidContactTemplate + ::dump(UInt step) { solid->onDump(); Model::dump(step); } /* -------------------------------------------------------------------------- */ -template -void CouplerSolidContactTemplate::dump(Real time, - UInt step) { +template +void AbstractCouplerSolidContactTemplate + ::dump(Real time, UInt step) { solid->onDump(); Model::dump(time, step); } +/* -------------------------------------------------------------------------- */ +/* Implementations specific to ContactMechanicsModel */ +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +template +MatrixType CouplerSolidContactTemplate::getMatrixType( + const ID & matrix_id) const { + if (matrix_id == "K") { + return _symmetric; + } + if (matrix_id == "M") { + return _symmetric; + } + return _mt_not_defined; +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactTemplate::assembleMatrix( + const ID & matrix_id) { + if (matrix_id == "K") { + this->assembleStiffnessMatrix(); + } else if (matrix_id == "M") { + this->solid->assembleMass(); + } +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactTemplate::corrector() { + AbstractCouplerSolidContactTemplate::corrector(); + + switch (this->method) { + case _static: + case _implicit_dynamic: { + auto & current_positions = this->contact->getContactDetector().getPositions(); + current_positions.copy(this->solid->getCurrentPosition()); + this->contact->search(); + break; + } + default: + break; + } +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactTemplate::assembleResidual() { + // computes the internal forces + switch (this->method) { + case _explicit_lumped_mass: { + auto & current_positions = this->contact->getContactDetector().getPositions(); + current_positions.copy(this->solid->getCurrentPosition()); + this->contact->search(); + break; + } + default: + break; + } + + this->assembleInternalForces(); + + auto & internal_force = this->solid->getInternalForce(); + auto & external_force = this->solid->getExternalForce(); + + auto & contact_force = this->contact->getInternalForce(); + + /* ------------------------------------------------------------------------ */ + this->getDOFManager().assembleToResidual("displacement", external_force, 1); + this->getDOFManager().assembleToResidual("displacement", internal_force, 1); + this->getDOFManager().assembleToResidual("displacement", contact_force, 1); +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactTemplate::assembleResidual( + const ID & residual_part) { + AKANTU_DEBUG_IN(); + + // contact->assembleInternalForces(); + + auto & internal_force = this->solid->getInternalForce(); + auto & external_force = this->solid->getExternalForce(); + auto & contact_force = this->contact->getInternalForce(); + + if ("external" == residual_part) { + this->getDOFManager().assembleToResidual("displacement", external_force, 1); + this->getDOFManager().assembleToResidual("displacement", contact_force, 1); + AKANTU_DEBUG_OUT(); + return; + } + + if ("internal" == residual_part) { + this->getDOFManager().assembleToResidual("displacement", internal_force, 1); + AKANTU_DEBUG_OUT(); + return; + } + + AKANTU_CUSTOM_EXCEPTION( + debug::SolverCallbackResidualPartUnknown(residual_part)); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactTemplate::assembleInternalForces() { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_INFO("Assemble the internal forces"); + + this->solid->assembleInternalForces(); + this->contact->assembleInternalForces(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void CouplerSolidContactTemplate:: + assembleStiffnessMatrix() { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); + + this->solid->assembleStiffnessMatrix(true); + + switch (this->method) { + case _static: + case _implicit_dynamic: { + this->contact->assembleStiffnessMatrix(); + break; + } + default: + break; + } + + AKANTU_DEBUG_OUT(); +} + } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index c59136979..63575195d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,601 +1,606 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Model of Solid Mechanics * * * @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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", std::shared_ptr dof_manager = nullptr, ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(bool need_to_reassemble = false); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); 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 matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); public: /// 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: /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(GhostType ghost_type) override; void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, ElementKind element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, ElementKind kind, GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(ElementKind kind); //! inverse operation of the flatten void inflateInternal(const std::string & field_name, const ElementTypeMapArray & field, ElementKind kind, GhostType ghost_type = _not_ghost); std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; 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; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material (by material name) inline const Material & getMaterial(const Element & element) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /** * @brief Returns the total energy for a given energy type * * Energy types of SolidMechanicsModel expected as argument are: * - `kinetic` * - `external work` * * Other energy types are passed on to the materials. All materials should * define a `potential` energy type. For additional energy types, see material * documentation. */ Real getEnergy(const std::string & energy_id); /// Compute energy for an element type and material index Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// Compute energy for an individual element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// Compute energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, // material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, material_selector, std::shared_ptr); void setMaterialSelector(std::shared_ptr material_selector) { this->material_selector = std::move(material_selector); } /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// mapping between material name and material internal id std::map materials_names_to_id; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// class defining of to choose a material std::shared_ptr material_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// tells if the material are instantiated flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// tells if the material are instantiated bool are_materials_instantiated{false}; friend class Material; - template friend class CouplerSolidContactTemplate; + template + friend class AbstractCouplerSolidContactTemplate; + template + friend class CouplerSolidContactTemplate; + template + friend class CouplerSolidContactInternodesTemplate; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.hh" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */