diff --git a/python/py_contact_mechanics_model.cc b/python/py_contact_mechanics_model.cc index 0325cc68c..931fe5f21 100644 --- a/python/py_contact_mechanics_model.cc +++ b/python/py_contact_mechanics_model.cc @@ -1,110 +1,109 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](ContactMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](ContactMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_contact_mechanics_model(py::module & mod) { py::class_(mod, "ContactMechanicsModelOptions") .def(py::init(), py::arg("analysis_method") = _explicit_contact); py::class_(mod, "ContactMechanicsModel", py::multiple_inheritance()) - .def(py::init, const ModelType>(), + .def(py::init, + const ModelType>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, - py::arg("id") = "contact_mechanics_model", py::arg("memory_id") = 0, - py::arg("dof_manager") = nullptr, + py::arg("id") = "contact_mechanics_model", py::arg("dof_manager") = nullptr, py::arg("model_type") = ModelType::_contact_mechanics_model) .def( "initFull", [](ContactMechanicsModel & self, const ContactMechanicsModelOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = ContactMechanicsModelOptions()) .def( "initFull", [](ContactMechanicsModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_function(search) .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function_nocopy(getExternalForce) .def_function_nocopy(getNormalForce) .def_function_nocopy(getTangentialForce) .def_function_nocopy(getInternalForce) .def_function_nocopy(getGaps) .def_function_nocopy(getNormals) .def_function_nocopy(getNodalArea) .def("dump", py::overload_cast<>(&ContactMechanicsModel::dump)) .def("dump", py::overload_cast(&ContactMechanicsModel::dump)) .def("dump", py::overload_cast( &ContactMechanicsModel::dump)) .def("dump", py::overload_cast( &ContactMechanicsModel::dump)); py::class_(mod, "ContactElement").def(py::init<>()); py::class_(mod, "GeometryUtils") .def_static( "normal", py::overload_cast &, const Element &, Vector &, bool>(&GeometryUtils::normal), py::arg("mesh"), py::arg("positions"), py::arg("element"), py::arg("normal"), py::arg("outward") = true) .def_static( "covariantBasis", py::overload_cast &, const Element &, const Vector &, Vector &, Matrix &>(&GeometryUtils::covariantBasis), py::arg("mesh"), py::arg("positions"), py::arg("element"), py::arg("normal"), py::arg("natural_projection"), py::arg("basis")) .def_static("curvature", &GeometryUtils::curvature) .def_static("contravariantBasis", &GeometryUtils::contravariantBasis, py::arg("covariant_basis"), py::arg("basis")) .def_static("realProjection", py::overload_cast &, const Vector &, const Element &, const Vector &, Vector &>( &GeometryUtils::realProjection), py::arg("mesh"), py::arg("positions"), py::arg("slave"), py::arg("element"), py::arg("normal"), py::arg("projection")) // .def_static("naturalProjection", &GeometryUtils::naturalProjection, // py::arg("mesh"), py::arg("positions"), py::arg("element"), // py::arg("real_projection"), py::arg("projection")) .def_static("isBoundaryElement", &GeometryUtils::isBoundaryElement); } } // namespace akantu diff --git a/python/py_model_couplers.cc b/python/py_model_couplers.cc index c2555eed4..08e814e19 100644 --- a/python/py_model_couplers.cc +++ b/python/py_model_couplers.cc @@ -1,73 +1,71 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def(#func_name, \ [](CouplerSolidContact & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](CouplerSolidContact & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void register_model_couplers(py::module & mod) { py::class_(mod, "CouplerSolidContactOptions") .def(py::init(), py::arg("analysis_method") = _explicit_dynamic_contact); py::class_(mod, "CouplerSolidContact") - .def(py::init, const ModelType>(), + .def(py::init, const ModelType>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, - py::arg("id") = "coupler_solid_contact", py::arg("memory_id") = 0, - py::arg("dof_manager") = nullptr, + py::arg("id") = "coupler_solid_contact", py::arg("dof_manager") = nullptr, py::arg("model_type") = ModelType::_coupler_solid_contact) .def("initFull", [](CouplerSolidContact & self, const CouplerSolidContactOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = CouplerSolidContactOptions()) .def("initFull", [](CouplerSolidContact & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def("setTimeStep", &CouplerSolidContact::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getSolidMechanicsModel", &CouplerSolidContact::getSolidMechanicsModel, py::return_value_policy::reference) .def("getContactMechanicsModel", &CouplerSolidContact::getContactMechanicsModel, py::return_value_policy::reference) .def("dump", py::overload_cast<>(&CouplerSolidContact::dump)) .def("dump", py::overload_cast(&CouplerSolidContact::dump)) .def("dump", py::overload_cast( &CouplerSolidContact::dump)) .def("dump", py::overload_cast( &CouplerSolidContact::dump)); } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 20b8eb9cb..c33124caf 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,353 +1,352 @@ /** * @file contact_detector.cc * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Fri Sep 21 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "contact_detector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -ContactDetector::ContactDetector(Mesh & mesh, const ID & id, - const MemoryID & memory_id) - : ContactDetector(mesh, mesh.getNodes(), id, memory_id) {} +ContactDetector::ContactDetector(Mesh & mesh, const ID & id) + : ContactDetector(mesh, mesh.getNodes(), id) {} /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, Array positions, - const ID & id, const MemoryID & memory_id) - : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), + const ID & id) + : Parsable(ParserType::_contact_detector, id), mesh(mesh), positions(0, mesh.getSpatialDimension()) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->positions.copy(positions); this->parseSection(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::parseSection() { const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); auto type = section.getParameterValue("type"); if (type == "implicit") { this->detection_type = _implicit; } else if (type == "explicit") { this->detection_type = _explicit; } else { AKANTU_ERROR("Unknown detection type : " << type); } this->projection_tolerance = section.getParameterValue("projection_tolerance"); this->max_iterations = section.getParameterValue("max_iterations"); this->extension_tolerance = section.getParameterValue("extension_tolerance"); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(Array & elements, Array & gaps, Array & normals, Array & tangents, Array & projections) { UInt surface_dimension = spatial_dimension - 1; this->mesh.fillNodesToElements(surface_dimension); this->computeMaximalDetectionDistance(); contact_pairs.clear(); SpatialGrid master_grid(spatial_dimension); SpatialGrid slave_grid(spatial_dimension); this->globalSearch(slave_grid, master_grid); this->localSearch(slave_grid, master_grid); this->createContactElements(elements, gaps, normals, tangents, projections); } /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { auto & master_list = surface_selector->getMasterList(); auto & slave_list = surface_selector->getSlaveList(); BBox bbox_master(spatial_dimension); this->constructBoundingBox(bbox_master, master_list); BBox bbox_slave(spatial_dimension); this->constructBoundingBox(bbox_slave, slave_list); auto && bbox_intersection = bbox_master.intersection(bbox_slave); AKANTU_DEBUG_INFO("Intersection BBox " << bbox_intersection); Vector center(spatial_dimension); bbox_intersection.getCenter(center); Vector spacing(spatial_dimension); this->computeCellSpacing(spacing); master_grid.setCenter(center); master_grid.setSpacing(spacing); this->constructGrid(master_grid, bbox_intersection, master_list); slave_grid.setCenter(center); slave_grid.setSpacing(spacing); this->constructGrid(slave_grid, bbox_intersection, slave_list); // search slave grid nodes in contactelement array and if they exits // and still have orthogonal projection on its associated master // facetremove it from the spatial grid or do not consider it for // local search, maybe better option will be to have spatial grid of // type node info and one of the variable of node info should be // facet already exits // so contact elements will be updated based on the above // consideration , this means only those contact elements will be // keep whose slave node is still in intersection bbox and still has // projection in its master facet // also if slave node is already exists in contact element and // orthogonal projection does not exits then search the associated // master facets with the current master facets within a given // radius , this is subjected to computational cost as searching // neighbbor cells can be more effective. } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { // local search // out of these array check each cell for closet node in that cell // and neighbouring cells find the actual orthogonally closet // check the projection of slave node on master facets connected to // the closet master node, if yes update the contact element with // slave node and master node and master surfaces connected to the // master node // these master surfaces will be needed later to update contact // elements /// find the closet master node for each slave node for (auto && cell_id : slave_grid) { /// loop over all the slave nodes of the current cell for (auto && slave_node : slave_grid.getCell(cell_id)) { bool pair_exists = false; Vector pos(spatial_dimension); for (UInt s : arange(spatial_dimension)) pos(s) = this->positions(slave_node, s); Real closet_distance = std::numeric_limits::max(); UInt closet_master_node; /// loop over all the neighboring cells of the current cell for (auto && neighbor_cell : cell_id.neighbors()) { /// loop over the data of neighboring cells from master grid for (auto && master_node : master_grid.getCell(neighbor_cell)) { /// check for self contact if (slave_node == master_node) continue; bool is_valid = true; Array elements; this->mesh.getAssociatedElements(slave_node, elements); for (auto & elem : elements) { if (elem.kind() != _ek_regular) continue; Vector connectivity = const_cast(this->mesh).getConnectivity(elem); auto node_iter = std::find(connectivity.begin(), connectivity.end(), master_node); if (node_iter != connectivity.end()) { is_valid = false; break; } } Vector pos2(spatial_dimension); for (UInt s : arange(spatial_dimension)) pos2(s) = this->positions(master_node, s); Real distance = pos.distance(pos2); if (distance <= closet_distance and is_valid) { closet_master_node = master_node; closet_distance = distance; pair_exists = true; } } } if (pair_exists) contact_pairs.push_back(std::make_pair(slave_node, closet_master_node)); } } } /* -------------------------------------------------------------------------- */ void ContactDetector::createContactElements( Array & contact_elements, Array & gaps, Array & normals, Array & tangents, Array & projections) { auto surface_dimension = spatial_dimension - 1; Real alpha; switch (detection_type) { case _explicit: { alpha = 1.0; break; } case _implicit: { alpha = -1.0; break; } default: AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type"); break; } for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; Vector slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector projection(projections.begin(surface_dimension)[slave_node]); Matrix tangent( tangents.begin(surface_dimension, spatial_dimension)[slave_node]); auto index = GeometryUtils::orthogonalProjection( mesh, positions, slave, elements, gap, projection, normal, tangent, alpha, this->max_iterations, this->projection_tolerance, this->extension_tolerance); // if a valid projection is not found on the patch of elements // index is -1 or if not a valid self contact, the contact element // is not created if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) { gap *= 0; normal *= 0; projection *= 0; tangent *= 0; continue; } // create contact element contact_elements.push_back(ContactElement(slave_node, elements[index])); } contact_pairs.clear(); } /* -------------------------------------------------------------------------- */ /*void ContactDetector::createContactElements(Array & contact_elements, Array & gaps, Array & normals, Array & projections) { auto surface_dimension = spatial_dimension - 1; Real alpha; switch (detection_type) { case _explicit: { alpha = 1.0; break; } case _implicit: { alpha = -1.0; break; } default: AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type"); break; } for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; Vector slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector projection(projections.begin(surface_dimension)[slave_node]); auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements, gap, projection, normal, alpha, this->max_iterations, this->projection_tolerance, this->extension_tolerance); // if a valid projection is not found on the patch of elements // index is -1 or if not a valid self contact, the contact element // is not created if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) { gap *= 0; normal *= 0; continue; } // create contact element contact_elements.push_back(ContactElement(slave_node, elements[index])); } contact_pairs.clear(); }*/ } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh index 8b475a2c5..b68264d24 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,219 +1,215 @@ /** * @file contact_detection.hh * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Tue Oct 23 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_bbox.hh" #include "aka_common.hh" #include "aka_grid_dynamic.hh" -#include "aka_memory.hh" #include "contact_element.hh" #include "element_class.hh" #include "element_group.hh" #include "fe_engine.hh" #include "mesh.hh" #include "mesh_io.hh" #include "parsable.hh" #include "surface_selector.hh" #include "geometry_utils.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_HH__ #define __AKANTU_CONTACT_DETECTOR_HH__ namespace akantu { enum class Surface { master, slave }; /* -------------------------------------------------------------------------- */ -class ContactDetector : private Memory, - public Parsable { +class ContactDetector : public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: - ContactDetector(Mesh &, const ID & id = "contact_detector", - const MemoryID & memory_id = 0); + ContactDetector(Mesh &, const ID & id = "contact_detector"); ContactDetector(Mesh &, Array positions, - const ID & id = "contact_detector", - const MemoryID & memory_id = 0); + const ID & id = "contact_detector"); ~ContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// performs all search steps void search(Array & contact_map, Array & gaps, Array & normals, Array & tangents, Array & projections); /// performs global spatial search to construct spatial grids void globalSearch(SpatialGrid &, SpatialGrid &); /// performs local search to find closet master node to a slave node void localSearch(SpatialGrid &, SpatialGrid &); /// 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(); /* ------------------------------------------------------------------------ */ /* Inline Methods */ /* ------------------------------------------------------------------------ */ public: /// checks whether the natural projection is valid or not inline bool checkValidityOfProjection(Vector &); /// extracts the coordinates of an element inline void coordinatesOfElement(const Element &, Matrix &); /// computes the optimal cell size for grid inline void computeCellSpacing(Vector &); /// constructs a grid containing nodes lying within bounding box inline void constructGrid(SpatialGrid &, BBox &, const Array &); /// constructs the bounding box based on nodes list inline void constructBoundingBox(BBox &, const Array &); /// computes the maximum in radius for a given mesh inline void computeMaximalDetectionDistance(); /// constructs the connectivity for a contact element inline Vector constructConnectivity(UInt &, const Element &); /// computes normal on an element inline void computeNormalOnElement(const Element &, Vector &); /// extracts vectors which forms the plane of element inline void vectorsAlongElement(const Element &, Matrix &); /// computes the gap between slave and its projection on master /// surface inline Real computeGap(Vector &, Vector &); /// filter boundary elements inline void filterBoundaryElements(Array & elements, Array & boundary_elements); /// checks whether self contact condition leads to a master element /// which is closet but not orthogonally opposite to slave surface //inline bool checkValidityOfSelfContact(const UInt &, ContactElement &); /// 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 &, const Real & , const Vector &); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the mesh AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// 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(Positions, positions, Array &); AKANTU_SET_MACRO(Positions, positions, Array); 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; /// Mesh Mesh & mesh; /// dimension of the model UInt spatial_dimension{0}; /// 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; /// contains the updated positions of the nodes Array positions; /// 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_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 2918f5af2..03b661571 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,789 +1,789 @@ /** * @file coontact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Contact mechanics model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "contact_mechanics_model.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "group_manager_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( - Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, + Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id, memory_id) { + : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif 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() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (this->parser.getLastParsedFile() != "") { 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(); // \TODO check if each resolution needs a initResolution() method } /* -------------------------------------------------------------------------- */ 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) { // 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); // TODO : uncomment when developing code for parallelization, // currently it addes the force twice for not ghost elements // hence source of error AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); contact_elements.resize(0); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->clear(); internal_array->resize(nb_nodes, 0.); }; 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); // intepenetration 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.size() != 0) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections previous_projections->clear(); previous_projections->resize(projections->size(), 0.); (*previous_projections).copy(*projections); // saving previous tangents previous_tangents->clear(); previous_tangents->resize(tangents->size(), 0.); (*previous_tangents).copy(*tangents); // saving previous tangential traction previous_tangential_tractions->clear(); previous_tangential_tractions->resize(tangential_tractions->size(), 0.); (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); 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->clear(); external_force->clear(); nodal_area->resize(nb_nodes, 0.); external_force->resize(nb_nodes, 0.); 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); /*switch (spatial_dimension) { case 1: { std::fill((*nodal_area).begin(), (*nodal_area).end(), 1.); break; } case 2: case 3: { this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), mesh.getElementGroup("contact_surface")); 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(); } break; } default: break; } this->external_force->clear();*/ 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; for (Int i = 0; i < indent; i++, space += 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 (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) { 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_IN(); 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); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { 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; - real_nodal_fields["normal_force"] = this->normal_force; - real_nodal_fields["tangential_force"] = this->tangential_force; - real_nodal_fields["blocked_dofs"] = this->blocked_dofs; - real_nodal_fields["normals"] = this->normals; - real_nodal_fields["tangents"] = this->tangents; - real_nodal_fields["gaps"] = this->gaps; - real_nodal_fields["areas"] = this->nodal_area; - real_nodal_fields["contact_state"] = this->contact_state; - real_nodal_fields["tangential_traction"] = this->tangential_tractions; + 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["contact_state"] = this->contact_state.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; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & , const std::string & , bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ 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 30561e567..eac9fab8f 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,383 +1,382 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Contact Mechanics * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, - const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, - std::shared_ptr dof_manager = nullptr, + const ID & id = "contact_mechanics_model", std::shared_ptr dof_manager = nullptr, const 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, NonLinearSolverType) 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) 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; /// 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 converted = 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 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 createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); /* ------------------------------------------------------------------------ */ /* 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: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// 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 Array & getContactElements() { 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 - Array * displacement{nullptr}; + std::unique_ptr> displacement; /// increment of displacement - Array * displacement_increment{nullptr}; + std::unique_ptr> displacement_increment; /// contact forces array - Array * internal_force{nullptr}; + std::unique_ptr> internal_force; /// external forces array - Array * external_force{nullptr}; + std::unique_ptr> external_force; /// normal force array - Array * normal_force{nullptr}; + std::unique_ptr> normal_force; /// friction force array - Array * tangential_force{nullptr}; + std::unique_ptr> tangential_force; /// friction traction array - Array * tangential_tractions{nullptr}; + std::unique_ptr> tangential_tractions; /// previous friction traction array - Array * previous_tangential_tractions{nullptr}; + std::unique_ptr> previous_tangential_tractions; /// boundary vector - Array * blocked_dofs{nullptr}; + std::unique_ptr> blocked_dofs; /// array to store gap between slave and master - Array * gaps{nullptr}; + std::unique_ptr> gaps; /// array to store normals from master to slave - Array * normals{nullptr}; + std::unique_ptr> normals; /// array to store tangents on the master element - Array * tangents{nullptr}; + std::unique_ptr> tangents; /// array to store previous tangents on the master element - Array * previous_tangents{nullptr}; + std::unique_ptr> previous_tangents; /// array to store nodal areas - Array * nodal_area{nullptr}; + std::unique_ptr> nodal_area; /// array to store stick/slip state : - Array * contact_state{nullptr}; + std::unique_ptr> contact_state; /// array to store previous projections in covariant basis - Array * previous_projections{nullptr}; + std::unique_ptr> previous_projections; // array to store projections in covariant basis - Array * projections{nullptr}; + 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 - Array * previous_master_elements{nullptr}; + 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/resolution.cc b/src/model/contact_mechanics/resolution.cc index 08cf26d19..1a422c603 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,225 +1,221 @@ /** * @file resolution.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Implementation of common part of the contact resolution class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution.hh" #include "contact_mechanics_model.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Resolution::Resolution(ContactMechanicsModel & model, const ID & id) - : Memory(id, model.getMemoryID()), - Parsable(ParserType::_contact_resolution, id), fem(model.getFEEngine()), - name(""), model(model) { + : Parsable(ParserType::_contact_resolution, id), id(id), + fem(model.getFEEngine()), model(model) { AKANTU_DEBUG_IN(); spatial_dimension = model.getMesh().getSpatialDimension(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Resolution::~Resolution() = default; /* -------------------------------------------------------------------------- */ void Resolution::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, "Friction Coefficient"); registerParam("is_master_deformable", is_master_deformable, bool(false), _pat_parsable | _pat_readable, "Is master surface deformable"); } /* -------------------------------------------------------------------------- */ void Resolution::printself(std::ostream & stream, int indent) const { - std::string space; - for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) - ; - + std::string space(indent, AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Contact Resolution " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); this->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces() { AKANTU_DEBUG_IN(); for (auto & element : model.getContactElements()) { auto nb_nodes = element.getNbNodes(); Vector local_fn(nb_nodes * spatial_dimension); computeNormalForce(element, local_fn); Vector local_ft(nb_nodes * spatial_dimension); computeTangentialForce(element, local_ft); Vector local_fc(nb_nodes * spatial_dimension); local_fc = local_fn + local_ft; assembleLocalToGlobalArray(element, local_fn, model.getNormalForce()); assembleLocalToGlobalArray(element, local_ft, model.getTangentialForce()); assembleLocalToGlobalArray(element, local_fc, model.getInternalForce()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleLocalToGlobalArray(const ContactElement & element, Vector & local, Array & global) { auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = const_cast(model.getMesh()).getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i - 1]; } return elem_conn; }; auto & surface_selector = model.getContactDetector().getSurfaceSelector(); auto & slave_list = surface_selector.getSlaveList(); auto & master_list = surface_selector.getMasterList(); auto connectivity = get_connectivity(element.slave, element.master); UInt nb_dofs = global.getNbComponent(); UInt nb_nodes = is_master_deformable ? connectivity.size() : 1; Real alpha = is_master_deformable ? 0.5: 1.; for (UInt i : arange(nb_nodes)) { UInt n = connectivity[i]; auto slave_result = std::find(slave_list.begin(), slave_list.end(), n); auto master_result = std::find(master_list.begin(), master_list.end(), n); for (UInt j : arange(nb_dofs)) { UInt offset_node = n * nb_dofs + j; global[offset_node] += alpha*local[i * nb_dofs + j]; } } } /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); auto & global_stiffness = const_cast(model.getDOFManager().getMatrix("K")); for (auto & element : model.getContactElements()) { auto nb_nodes = element.getNbNodes(); Matrix local_kn(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); computeNormalModuli(element, local_kn); assembleLocalToGlobalMatrix(element, local_kn, global_stiffness); Matrix local_kt(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); computeTangentialModuli(element, local_kt); assembleLocalToGlobalMatrix(element, local_kt, global_stiffness); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleLocalToGlobalMatrix(const ContactElement & element, const Matrix & local, SparseMatrix & global) { auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = const_cast(model.getMesh()).getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i - 1]; } return elem_conn; }; auto connectivity = get_connectivity(element.slave, element.master); auto nb_dofs = spatial_dimension; UInt nb_nodes = is_master_deformable ? connectivity.size() : 1; UInt total_nb_dofs = nb_dofs * nb_nodes; std::vector equations; for (UInt i : arange(connectivity.size())) { UInt conn = connectivity[i]; for (UInt j : arange(nb_dofs)) { equations.push_back(conn * nb_dofs + j); } } for (UInt i : arange(total_nb_dofs)) { UInt row = equations[i]; for (UInt j : arange(total_nb_dofs)) { UInt col = equations[j]; global.add(row, col, local(i, j)); } } } /* -------------------------------------------------------------------------- */ void Resolution::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void Resolution::afterSolveStep(__attribute__((unused)) bool converged) {} } // namespace akantu diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh index 0d7f94745..c746ab842 100644 --- a/src/model/contact_mechanics/resolution.hh +++ b/src/model/contact_mechanics/resolution.hh @@ -1,228 +1,234 @@ /** * @file resolution.hh * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Mother class for all contact resolutions * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" -#include "aka_memory.hh" #include "parsable.hh" #include "parser.hh" #include "fe_engine.hh" #include "contact_element.hh" #include "resolution_utils.hh" #include "geometry_utils.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_HH__ #define __AKANTU_RESOLUTION_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class ContactMechanicsModel; } // namespace akantu namespace akantu { /** * Interface of all contact resolutions * Prerequisites for a new resolution * - inherit from this class * - implement the following methods: * \code * * virtual void computeNormalForce(); * virtual void computeTangentialForce(); * virtual void computeNormalModuli(); * virtual void computeTangentialModuli(); * * \endcode * */ -class Resolution : public Memory, - public Parsable { +class Resolution : public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ public: /// instantiate contact resolution with defaults Resolution(ContactMechanicsModel & model, const ID & id = ""); /// Destructor ~Resolution() override; protected: void initialize(); /// computes coordinates of a given element void computeCoordinates(const Element & , Matrix &); /* ------------------------------------------------------------------------ */ /* Functions that resolutions should reimplement for force */ /* ------------------------------------------------------------------------ */ public: /// computes the force vector due to normal traction virtual void computeNormalForce(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } /// computes the tangential force vector due to frictional traction virtual void computeTangentialForce(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Functions that resolutions should reimplement for stiffness */ /* ------------------------------------------------------------------------ */ public: /// compute the normal moduli due to normal traction virtual void computeNormalModuli(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Matrix & ) { AKANTU_TO_IMPLEMENT(); } /// compute the tangent moduli due to tangential traction virtual void computeTangentialModuli(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Matrix & ) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// assemble the residual for this resolution void assembleInternalForces(GhostType ghost_type); /// assemble the stiffness matrix for this resolution void assembleStiffnessMatrix(GhostType ghost_type); private: /// assemble the residual for this resolution void assembleInternalForces(); /// assemble the local array to global array for a contact element void assembleLocalToGlobalArray(const ContactElement & , Vector & , Array & ); /// assemble the local stiffness to global stiffness for a contact element void assembleLocalToGlobalMatrix(const ContactElement &, const Matrix &, SparseMatrix &); public: virtual void beforeSolveStep(); virtual void afterSolveStep(bool converged = true); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + AKANTU_GET_MACRO(ID, id, const ID &); + public: /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - + ID id; + /// friction coefficient : mu Real mu; /// spatial dimension UInt spatial_dimension; /// is master surface deformable bool is_master_deformable; /// Link to the fe engine object in the model FEEngine & fem; /// resolution name std::string name; /// model to which the resolution belong ContactMechanicsModel & model; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Resolution & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { using ResolutionFactory = Factory; /// macaulay bracket to convert positive gap to zero template T macaulay(T var) {return var < 0 ? 0 : var; } template T heaviside(T var) {return var < 0 ? 0 : 1.0; } } // namespace akantu #define INSTANTIATE_RESOLUTION_ONLY(res_name) \ class res_name #define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name) \ [](UInt dim, const ID &, ContactMechanicsModel & model, \ const ID & id) -> std::unique_ptr { \ switch (dim) { \ case 1: \ return std::make_unique(model, id); \ case 2: \ return std::make_unique(model, id); \ case 3: \ return std::make_unique(model, id); \ default: \ AKANTU_EXCEPTION("The dimension " \ << dim << "is not a valid dimension for the contact resolution " \ << #id); \ } \ } #define INSTANTIATE_RESOLUTION(id, res_name) \ INSTANTIATE_RESOLUTION_ONLY(res_name); \ static bool resolution_is_alocated_##id[[gnu::unused]] = \ ResolutionFactory::getInstance().registerAllocator( \ #id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name)) #endif /* __AKANTU_RESOLUTION_HH__ */ diff --git a/src/model/model.hh b/src/model/model.hh index f6468cb27..cd19bb935 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,378 +1,378 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Interface of a model * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_named_argument.hh" #include "fe_engine.hh" #include "mesh.hh" #include "model_options.hh" #include "model_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MODEL_HH_ #define AKANTU_MODEL_HH_ namespace akantu { class SynchronizerRegistry; class Parser; class DumperIOHelper; class DOFManager; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class Model : public ModelSolver, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Normal constructor where the DOFManager is created internally Model(Mesh & mesh, const ModelType & type, UInt dim = _all_dimensions, const ID & id = "model"); /// Model constructor the the dof manager is created externally, for example /// in a ModelCoupler Model(Mesh & mesh, const ModelType & type, std::shared_ptr dof_manager, UInt dim = _all_dimensions, - const ID & id = "model", const MemoryID & memory_id = 0); + const ID & id = "model"); ~Model() override; using FEEngineMap = std::map>; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: virtual void initFullImpl(const ModelOptions & options); public: template std::enable_if_t::value> initFull(pack &&... _pack) { switch (this->model_type) { #ifdef AKANTU_SOLID_MECHANICS case ModelType::_solid_mechanics_model: this->initFullImpl(SolidMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_COHESIVE_ELEMENT case ModelType::_solid_mechanics_model_cohesive: this->initFullImpl(SolidMechanicsModelCohesiveOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_HEAT_TRANSFER case ModelType::_heat_transfer_model: this->initFullImpl(HeatTransferModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_EMBEDDED case ModelType::_embedded_model: this->initFullImpl(EmbeddedInterfaceModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_CONTACT_MECHANICS case ModelType::_contact_mechanics_model: this->initFullImpl(ContactMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_MODEL_COUPLERS case ModelType::_coupler_solid_contact: this->initFullImpl(CouplerSolidContactOptions{ use_named_args, std::forward(_pack)...}); break; case ModelType::_coupler_solid_cohesive_contact: this->initFullImpl(CouplerSolidCohesiveContactOptions{ use_named_args, std::forward(_pack)...}); break; #endif default: this->initFullImpl(ModelOptions{use_named_args, std::forward(_pack)...}); } } template std::enable_if_t::value> initFull(pack &&... _pack) { this->initFullImpl(std::forward(_pack)...); } /// initialize a new solver if needed void initNewSolver(const AnalysisMethod & method); protected: /// get some default values for derived classes virtual std::tuple getDefaultSolverID(const AnalysisMethod & method) = 0; virtual void initModel() = 0; virtual void initFEEngineBoundary(); /// function to print the containt of the class void printself(std::ostream & /*stream*/, int /*indent*/ = 0) const override{}; public: /* ------------------------------------------------------------------------ */ /* Access to the dumpable interface of the boundaries */ /* ------------------------------------------------------------------------ */ /// Dump the data for a given group void dumpGroup(const std::string & group_name); void dumpGroup(const std::string & group_name, const std::string & dumper_name); /// Dump the data for all boundaries void dumpGroup(); /// Set the directory for a given group void setGroupDirectory(const std::string & directory, const std::string & group_name); /// Set the directory for all boundaries void setGroupDirectory(const std::string & directory); /// Set the base name for a given group void setGroupBaseName(const std::string & basename, const std::string & group_name); /// Get the internal dumper of a given group DumperIOHelper & getGroupDumper(const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Function for non local capabilities */ /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(__attribute__((unused)) ElementTypeMapReal & criterion) { AKANTU_TO_IMPLEMENT(); } protected: template void allocNodalField(std::unique_ptr> & array, UInt nb_component, const ID & name) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID &) /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries(){}; /// return the fem object associated with a provided name inline FEEngine & getFEEngine(const ID & name = "") const; /// return the fem boundary object associated with a provided name virtual FEEngine & getFEEngineBoundary(const ID & name = ""); /// register a fem object associated with name template inline void registerFEEngineObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEEngineObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); /// return the fem object associated with a provided name template inline FEEngineClass & getFEEngineClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template inline FEEngineClass & getFEEngineClassBoundary(std::string name = ""); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /* ------------------------------------------------------------------------ */ /* Pack and unpack hexlper functions */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbIntegrationPoints(const Array & elements, const ID & fem_id = ID()) const; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ void setTextModeToDumper(); virtual void addDumpGroupFieldToDumper(const std::string & field_id, std::shared_ptr field, DumperIOHelper & dumper); virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void setBaseName(const std::string & field_id); virtual void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, ElementKind element_kind, bool padding_flag); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual std::shared_ptr createNodalFieldReal(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldUInt(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createElementalField( const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind /*kind*/) { return nullptr; } void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); /* ------------------------------------------------------------------------ */ virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); /* ------------------------------------------------------------------------ */ virtual void dump(); virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend std::ostream & operator<<(std::ostream & /*stream*/, const Model & /*_this*/); ID id; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// Mesh Mesh & mesh; /// Spatial dimension of the problem UInt spatial_dimension; /// the main fem object present in all models FEEngineMap fems; /// the fem object present in all models for boundaries FEEngineMap fems_boundary; /// default fem object std::string default_fem; /// parser to the pointer to use Parser & parser; /// default ElementKind for dumper ElementKind dumper_default_element_kind{_ek_regular}; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "model_inline_impl.hh" #endif /* AKANTU_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 deed306eb..2b16d11a7 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,655 +1,654 @@ /** * @file coupler_solid_cohesive_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics cohesive and conatct mechanics * model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_cohesive_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidCohesiveContact::CouplerSolidCohesiveContact( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject( "CohesiveFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_cohesive_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_cohesive_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModelCohesive(mesh, Model::spatial_dimension, - "solid_mechanics_model_cohesive", 0, - this->dof_manager); + "solid_mechanics_model_cohesive", this->dof_manager); contact = new ContactMechanicsModel( mesh.getMeshFacets(), Model::spatial_dimension, "contact_mechanics_model", - 0, this->dof_manager); + this->dof_manager); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidCohesiveContact::~CouplerSolidCohesiveContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); 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); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initModel() { getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidCohesiveContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::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); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidCohesiveContact::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); } /*switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); }*/ } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidCohesiveContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidCohesiveContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_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::_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 CouplerSolidCohesiveContact::assembleResidual() { // computes the internal forces 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); /*auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn(const_cast(mesh).getConnectivity(master)); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i - 1]; } return elem_conn; }; switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & element : contact->getContactElements()) { for (auto & conn : get_connectivity(element.slave, element.master)) { for (auto dim : arange(spatial_dimension)) { external_force(conn, dim) = contact_force(conn, dim); } } } } default: break; }*/ /* ------------------------------------------------------------------------ */ /*this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; }*/ } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); /*auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn(const_cast(mesh).getConnectivity(master)); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i - 1]; } return elem_conn; }; switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & element : contact->getContactElements()) { for (auto & conn : get_connectivity(element.slave, element.master)) { for (auto dim : arange(spatial_dimension)) { external_force(conn, dim) = contact_force(conn, dim); } } } } default: break; } if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } AKANTU_DEBUG_OUT(); return; }*/ 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(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::predictor() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); switch (method) { case _static: case _explicit_lumped_mass: { auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(solid->getCurrentPosition()); contact->search(); break; } default: break; } /*switch (method) { case _explicit_dynamic_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } contact->search(); break; } default: break; }*/ } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::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; } /*switch (method) { case _implicit_contact: case _explicit_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } contact->search(); break; } default: break; }*/ } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::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(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::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); } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidCohesiveContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); switch (method) { case _static: case _implicit_dynamic: { contact->assembleStiffnessMatrix(); break; } default: break; } /*switch (method) { case _implicit_contact: { contact->assembleStiffnessMatrix(); break; } default: break; }*/ AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; field = solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; if (field_name == "contact_force" or field_name == "normals" or field_name == "normal_force" or field_name == "tangential_force" or field_name == "contact_state" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { field = contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = solid->createNodalFieldBool(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index 2649346c5..09c237ad9 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,468 +1,466 @@ /** * @file coupler_solid_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct mechanics * model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidContact::CouplerSolidContact( - Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, + Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id, memory_id) { + : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidContact", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); - solid = - new SolidMechanicsModel(mesh, Model::spatial_dimension, - "solid_mechanics_model", 0, this->dof_manager); + solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, + "solid_mechanics_model", this->dof_manager); contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, - "contact_mechanics_model", 0, - this->dof_manager); + "contact_mechanics_model", this->dof_manager); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidContact::~CouplerSolidContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); solid->initFull(_analysis_method = this->method); contact->initFull(_analysis_method = this->method); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::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); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidContact::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); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_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::_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 CouplerSolidContact::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); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::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(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::predictor() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::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; } } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::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(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::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); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::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(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; field = solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; if (field_name == "contact_force" or field_name == "normals" or field_name == "normal_force" or field_name == "tangential_force" or field_name == "contact_state" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { field = contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = solid->createNodalFieldBool(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.hh b/src/model/model_couplers/coupler_solid_contact.hh index 3e2c7832b..5d9e01a66 100644 --- a/src/model/model_couplers/coupler_solid_contact.hh +++ b/src/model/model_couplers/coupler_solid_contact.hh @@ -1,281 +1,280 @@ /** * @file coupler_solid_contact_explicit.hh * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu Jan 17 2019 * * @brief class for coupling of solid mechanics and conatct mechanics * model in explicit * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "contact_mechanics_model.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COUPLER_SOLID_CONTACT_HH__ #define __AKANTU_COUPLER_SOLID_CONTACT_HH__ /* ------------------------------------------------------------------------ */ /* Coupling : Solid Mechanics / Contact Mechanics */ /* ------------------------------------------------------------------------ */ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; class DOFManager; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class CouplerSolidContact : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: CouplerSolidContact( Mesh & mesh, UInt spatial_dimension = _all_dimensions, - const ID & id = "coupler_solid_contact", const MemoryID & memory_id = 0, - std::shared_ptr dof_manager = nullptr, + const ID & id = "coupler_solid_contact", std::shared_ptr dof_manager = nullptr, const ModelType model_type = ModelType::_coupler_solid_contact); ~CouplerSolidContact() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// initialize the modelType void initModel() 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(); 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() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// callback for the model to instantiate the matricess when needed void initSolver(TimeStepSolverType, NonLinearSolverType) 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; public: bool isDefaultSolverExplicit() { return method == _explicit_dynamic_contact; } /* ------------------------------------------------------------------------ */ public: // DataAccessor UInt getNbData(const Array &, const SynchronizationTag &) const override { return 0; } void packData(CommunicationBuffer &, const Array &, const SynchronizationTag &) const override {} void unpackData(CommunicationBuffer &, const Array &, const SynchronizationTag &) override {} // DataAccessor nodes UInt getNbData(const Array &, const SynchronizationTag &) const override { return 0; } void packData(CommunicationBuffer &, const Array &, const SynchronizationTag &) const override {} void unpackData(CommunicationBuffer &, const Array &, const SynchronizationTag &) override {} /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// 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 ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, 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 solid mechanics model AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &); /// get the contact mechanics model AKANTU_GET_MACRO(ContactMechanicsModel, *contact, ContactMechanicsModel &); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr 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; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// solid mechanics model SolidMechanicsModel * solid{nullptr}; /// contact mechanics model ContactMechanicsModel * contact{nullptr}; /// Array * displacement{nullptr}; /// Array * displacement_increment{nullptr}; /// external forces array Array * external_force{nullptr}; bool search_for_contact; UInt step; }; } // namespace akantu #endif /* __COUPLER_SOLID_CONTACT_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index be928bd13..11e8e63fc 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,587 +1,582 @@ /** * @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: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) 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(); /// applya 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); 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; - - virtual void dump(const std::string & dumper_name); - - virtual void dump(const std::string & dumper_name, UInt step); - 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: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// 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 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(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for an element Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// compute the energy for an element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// compute the 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, MaterialSelector &); 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>>; /// map a registered internals to be flattened for dump purposes 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; }; /* -------------------------------------------------------------------------- */ 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_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index 656247b01..7b85c7c4c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,726 +1,725 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Fabian Barras * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Solid mechanics model for cohesive elements * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), *model.cohesive_synchronizer) {} /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto *cohesive_nodes_event = dynamic_cast(&nodes_event); if (cohesive_nodes_event == nullptr) { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); auto & nodes_flags = mesh_accessor.getNodesFlags(); auto nb_old_nodes = nodes_flags.size(); nodes_flags.resize(nb_old_nodes + local_nb_new_nodes); for (auto && data : zip(old_nodes, new_nodes)) { UInt old_node; UInt new_node; std::tie(old_node, new_node) = data; nodes_flags(new_node) = nodes_flags(old_node); } model.updateCohesiveSynchronizers(elements_event); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } Vector nb_new_stuff = {nb_new_nodes, elements_event.getList().size()}; const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_stuff, SynchronizerOperation::_sum); if (nb_new_stuff(1) > 0) { mesh.sendEvent(elements_event); } if (nb_new_stuff(0) > 0) { mesh.sendEvent(nodes_event); // mesh.sendEvent(global_ids_updater.getChangedNodeEvent()); } return std::make_tuple(nb_new_stuff(0), nb_new_stuff(1)); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) : SolidMechanicsModel(mesh, dim, id, dof_manager, - ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = aka::as_type(options); this->is_extrinsic = smmc_options.is_extrinsic; inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = aka::as_type(mesh_facets.getElementSynchronizer()); // synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, SynchronizationTag::_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (not are_materials_instantiated) { instantiateMaterials(); } /// find the first cohesive material UInt cohesive_index = UInt(-1); for (auto && material : enumerate(materials)) { if (dynamic_cast(std::get<1>(material).get()) != nullptr) { cohesive_index = std::get<0>(material); break; } } if (cohesive_index == UInt(-1)) { AKANTU_EXCEPTION("No cohesive materials in the material input file"); } material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = aka::as_type(*materials[mat_index]); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initMaterials(); if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.empty()) { continue; } type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (const auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) { continue; } /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) { continue; } Matrix quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; // for (UInt q = 0; q < nb_quad_per_facet; ++q) { // for (UInt s = 0; s < Model::spatial_dimension; ++s) { // el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + // f * nb_quad_per_facet + q, // s) = quad_f(global_facet * nb_quad_per_facet + q, // s); // } // } //} } } } /// loop over non cohesive materials for (auto && material : materials) { if (aka::is_of_type(material)) { continue; } /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = aka::as_type(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { if (not aka::is_of_type(material)) { /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } } this->synchronize(SynchronizationTag::_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); for (auto & mat : materials) { if (aka::is_of_type(mat)) { /// check which not ghost cohesive elements are to be created auto * mat_cohesive = aka::as_type(mat.get()); mat_cohesive->checkInsertion(); } } /// communicate data among processors // this->synchronize(SynchronizationTag::_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onElementsAdded(element_list, event); if (is_extrinsic) { resizeFacetStress(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onNodesAdded(new_nodes, event); const CohesiveNewNodesEvent * cohesive_event; if ((cohesive_event = dynamic_cast(&event)) == nullptr) { return; } const auto & old_nodes = cohesive_event->getOldNodesList(); auto copy = [this, &new_nodes, &old_nodes](auto & arr) { UInt new_node; UInt old_node; auto view = make_view(arr, spatial_dimension); auto begin = view.begin(); for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; auto old_ = begin + old_node; auto new_ = begin + new_node; *new_ = *old_; } }; copy(*displacement); copy(*blocked_dofs); if (velocity) { copy(*velocity); } if (acceleration) { copy(*acceleration); } if (current_position) { copy(*current_position); } if (previous_displacement) { copy(*previous_displacement); } // if (external_force) // copy(*external_force); // if (internal_force) // copy(*internal_force); if (displacement_increment) { copy(*displacement_increment); } copy(getDOFManager().getSolution("displacement")); // this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::afterSolveStep(bool converged) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ if (converged) { for (auto & mat : materials) { if (mat->isFiniteDeformation()) { mat->computeAllCauchyStresses(_not_ghost); } } } SolidMechanicsModel::afterSolveStep(converged); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << "\n"; SolidMechanicsModel::printself(stream, indent + 2); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu