diff --git a/packages/model_couplers.cmake b/packages/model_couplers.cmake index 05348aeac..cf800957b 100644 --- a/packages/model_couplers.cmake +++ b/packages/model_couplers.cmake @@ -1,48 +1,50 @@ #=============================================================================== # @file model_couplers.cmake # # @author Mohit Pundir # # @date creation: Sun Sep 30 2018 # @date last modification: Sun Sep 28 2018 # # @brief package description for model couplers # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(model_couplers DESCRIPTION "Use Model Couplers package of Akantu" DEPENDS contact_mechanics) package_declare_sources(model_couplers model/model_couplers/coupler_solid_contact.hh model/model_couplers/coupler_solid_contact.cc + model/model_couplers/coupler_solid_cohesive_contact.hh + model/model_couplers/coupler_solid_cohesive_contact.cc ) package_declare_documentation_files(model_couplers # ) package_declare_documentation(model_couplers "This package activates the model couplers within Akantu. " "It has no additional dependencies." ) diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index bbbaeb9b9..16582b562 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,613 +1,614 @@ /** * @file mesh.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief class handling meshes * * @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_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.cc" #include "mesh.hh" #include "mesh_global_data_updater.hh" #include "mesh_io.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" #include "periodic_node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_field.hh" #include "dumper_internal_material_field.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator) : Memory(id, memory_id), GroupManager(*this, id + ":group_manager", memory_id), MeshData("mesh_data", id, memory_id), connectivities("connectivities", id, memory_id), ghosts_counters("ghosts_counters", id, memory_id), normals("normals", id, memory_id), spatial_dimension(spatial_dimension), size(spatial_dimension, 0.), bbox(spatial_dimension), bbox_local(spatial_dimension), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); this->nodes_flags = std::make_shared>(0, 1, NodeFlag::_normal, id + ":nodes_flags"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id, memory_id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (mesh_facets) { AKANTU_DEBUG_OUT(); return *mesh_facets; } mesh_facets = std::make_unique(spatial_dimension, this->nodes, getID() + ":" + id, getMemoryID()); mesh_facets->mesh_parent = this; mesh_facets->is_mesh_facets = true; mesh_facets->nodes_flags = this->nodes_flags; mesh_facets->nodes_global_ids = this->nodes_global_ids; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } auto & mesh_phys_data = this->getData("physical_names"); auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes(spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) tolerance *= norm_barycenter; const auto & element_to_facet = mesh_facets->getElementToSubelement( element.type, element.ghost_type); Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view( barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto facet = std::get<0>(data); if (element_to_facet(facet)[1] == ElementNull) return false; auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name phys_data(Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}) = mesh_phys_data(element); }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { AKANTU_DEBUG_ASSERT(not is_distributed, "You cannot read a mesh that is already distributed"); MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); auto types = this->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); auto it = types.begin(); auto last = types.end(); if (it == last) AKANTU_DEBUG_WARNING( "The mesh contained in the file " << filename << " does not seem to be of the good dimension." << " No element of dimension " << spatial_dimension << " where read."); this->makeReady(); } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::makeReady() { this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_flags->resize(nodes->size(), NodeFlag::_normal); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); bbox_local.reset(); for (auto & pos : make_view(*nodes, spatial_dimension)) { // if(!isPureGhostNode(i)) bbox_local += pos; } if (this->is_distributed) { bbox = bbox_local.allSum(*communicator); } else { bbox = bbox_local; } size = bbox.size(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_not_defined, _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) continue; auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [&](UInt l) -> UInt { return this->getNodeGlobalId(l); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name, const std::string & group_name) { if (group_name == "all") return this->getDumper(dumper_name); else return element_groups[group_name]->getDumper(dumper_name); } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & arrays, const ElementKind & element_kind) { ElementTypeMap nb_data_per_elem; for (auto type : elementTypes(spatial_dimension, _not_ghost, element_kind)) { UInt nb_elements = this->getNbElement(type); auto & array = arrays(type); nb_data_per_elem(type) = array.getNbComponent() * array.size(); nb_data_per_elem(type) /= nb_elements; } return nb_data_per_elem; } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind) { dumper::Field * field = nullptr; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal, element_kind); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); #endif /* -------------------------------------------------------------------------- */ void Mesh::distributeImpl( Communicator & communicator, std::function edge_weight_function, std::function vertex_weight_function) { AKANTU_DEBUG_ASSERT(is_distributed == false, "This mesh is already distribute"); this->communicator = &communicator; this->element_synchronizer = std::make_unique( *this, this->getID() + ":element_synchronizer", this->getMemoryID(), true); this->node_synchronizer = std::make_unique( *this, this->getID() + ":node_synchronizer", this->getMemoryID(), true); Int psize = this->communicator->getNbProc(); if (psize > 1) { #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize, edge_weight_function, vertex_weight_function); MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition); } else { MeshUtilsDistribution::distributeMeshCentralized(*this, 0); } #else if (psize > 1) { AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool"); } #endif } // if (psize > 1) this->is_distributed = true; this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) { for (const auto & node : node_list) for (const auto & element : *nodes_to_elements[node]) elements.push_back(element); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const UInt & node, Array & elements) { for (const auto & element : *nodes_to_elements[node]) elements.push_back(element); } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements() { Element e; UInt nb_nodes = nodes->size(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) this->nodes_to_elements[n]->clear(); else this->nodes_to_elements[n] = std::make_unique>(); } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); Array::const_iterator> conn_it = connectivities(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements(UInt dimension) { Element e; UInt nb_nodes = nodes->size(); + for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) this->nodes_to_elements[n]->clear(); else this->nodes_to_elements[n] = std::make_unique>(); } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); Array::const_iterator> conn_it = connectivities(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ std::tuple Mesh::updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { if (global_data_updater) return this->global_data_updater->updateData(nodes_event, elements_event); else { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } } /* -------------------------------------------------------------------------- */ void Mesh::registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->global_data_updater = std::move(global_data_updater); } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this); auto & remove_list = event.getList(); auto & new_numbering = event.getNewNumbering(); for (auto && el : elements) { if (el.ghost_type != _not_ghost) { auto & count = ghosts_counters(el); --count; if (count > 0) continue; } remove_list.push_back(el); if (not last_element.exists(el.type, el.ghost_type)) { UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); last_element(nb_element - 1, el.type, el.ghost_type); auto & numbering = new_numbering.alloc(nb_element, 1, el.type, el.ghost_type); for (auto && pair : enumerate(numbering)) { std::get<1>(pair) = std::get<0>(pair); } } UInt & pos = last_element(el.type, el.ghost_type); auto & numbering = new_numbering(el.type, el.ghost_type); numbering(el.element) = UInt(-1); numbering(pos) = el.element; --pos; } this->sendEvent(event); } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 3c4f04e3c..3ce57bd55 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,476 +1,480 @@ /** * @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, UInt memory_id) : ContactDetector(mesh, mesh.getNodes(), id, memory_id) { } /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, Array positions, const ID & id, UInt memory_id) : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), mesh(mesh) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->positions = positions; - this->mesh.fillNodesToElements(this->spatial_dimension - 1); + //this->mesh.fillNodesToElements(this->spatial_dimension - 1); + this->parseSection(); this->computeMaximalDetectionDistance(); 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); } surfaces[Surface::master] = section.getParameterValue("master"); surfaces[Surface::slave ] = section.getParameterValue("slave"); two_pass_algorithm = section.getParameterValue("two_pass_algorithm"); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(std::map & contact_map) { 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); //if (two_pass_algorithm) { // this->localSearch(master_grid, slave_grid); //} this->constructContactMap(contact_map); } /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { auto & master_list = mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes(); auto & slave_list = mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes(); 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)) { + + if (slave_node == master_node) + continue; 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) { 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::constructContactMap(std::map & contact_map) { auto surface_dimension = spatial_dimension - 1; std::map previous_contact_map = contact_map; contact_map.clear(); auto get_index = [&](auto & gaps, auto & projections) { UInt index; Real gap_min = std::numeric_limits::max(); UInt counter = 0; for (auto && values : zip(gaps, make_view(projections, surface_dimension))) { auto & gap = std::get<0>(values); auto & projection = std::get<1>(values); /// todo adhoc fix to assign a master element in case the /// projection does not lie in the extended element. As it is /// tolerance based bool is_valid = this->checkValidityOfProjection(projection); if (is_valid and gap <= gap_min) { gap_min = gap; index = counter; } counter++; } if (index >= gaps.size()) { auto gap_min_it = std::min_element(gaps.begin(), gaps.end()); auto index_it = std::find(gaps.begin(), gaps.end(), *gap_min_it); index = *index_it; } return index; }; auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = this->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; }; for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); Array gaps(elements.size(), 1, "gaps"); Array normals(elements.size(), spatial_dimension, "normals"); Array projections(elements.size(), surface_dimension, "projections"); this->computeOrthogonalProjection(slave_node, elements, normals, gaps, projections); auto index = get_index(gaps, projections); auto connectivity = get_connectivity(slave_node, elements[index]); contact_map[slave_node].setMaster(elements[index]); contact_map[slave_node].setGap(gaps[index]); contact_map[slave_node].setNormal(Vector(normals.begin(spatial_dimension)[index], true)); contact_map[slave_node].setProjection(Vector(projections.begin(surface_dimension)[index], true)); contact_map[slave_node].setConnectivity(connectivity); Matrix tangents(surface_dimension, spatial_dimension); this->computeTangentsOnElement(contact_map[slave_node].master, contact_map[slave_node].projection, tangents); switch (spatial_dimension) { case 2: { Vector e_z(3); e_z[0] = 0.; e_z[1] = 0.; e_z[2] = 1.; Vector tangent(3); tangent[0] = tangents(0, 0); tangent[1] = tangents(0, 1); tangent[2] = 0.; auto exp_normal = e_z.crossProduct(tangent); auto & cal_normal = contact_map[slave_node].normal; auto ddot = cal_normal.dot(exp_normal); if (ddot < 0) { tangents *= -1.0; } break; } case 3: { auto tang_trans = tangents.transpose(); auto tang1 = Vector(tang_trans(0)); auto tang2 = Vector(tang_trans(1)); auto tang1_cross_tang2 = tang1.crossProduct(tang2); auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); auto & cal_normal = contact_map[slave_node].normal; auto ddot = cal_normal.dot(exp_normal); if (ddot < 0) { tang_trans(1) *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } contact_map[slave_node].setTangent(tangents); auto search = previous_contact_map.find(slave_node); if (search != previous_contact_map.end()) { auto previous_projection = previous_contact_map[slave_node].getPreviousProjection(); contact_map[slave_node].setPreviousProjection(previous_projection); auto previous_traction = previous_contact_map[slave_node].getTraction(); contact_map[slave_node].setTraction(previous_traction); } else { Vector previous_projection(surface_dimension, 0.); contact_map[slave_node].setPreviousProjection(previous_projection); Vector previous_traction(surface_dimension, 0.); contact_map[slave_node].setTraction(previous_traction); } } contact_pairs.clear(); } /* -------------------------------------------------------------------------- */ void ContactDetector::computeOrthogonalProjection(const UInt & node, const Array & elements, Array & normals, Array & gaps, Array & projections) { Vector query(spatial_dimension); for (UInt s: arange(spatial_dimension)) query(s) = this->positions(node, s); for (auto && values : zip( elements, gaps, make_view(normals , spatial_dimension), make_view(projections, spatial_dimension - 1))) { const auto & element = std::get<0>(values); auto & gap = std::get<1>(values); auto & normal = std::get<2>(values); auto & projection = std::get<3>(values); this->computeNormalOnElement(element, normal); Vector real_projection(spatial_dimension); this->computeProjectionOnElement(element, normal, query, projection, real_projection); gap = this->computeGap(query, real_projection, normal); } } /* -------------------------------------------------------------------------- */ void ContactDetector::computeProjectionOnElement(const Element & element, const Vector & normal, const Vector & query, Vector & natural_projection, Vector & real_projection) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); Matrix coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(element, coords); Vector point(coords(0)); Real alpha = (query - point).dot(normal); real_projection = query - alpha * normal; this->computeNaturalProjection(element, real_projection, natural_projection); } /* -------------------------------------------------------------------------- */ void ContactDetector::computeNaturalProjection(const Element & element, Vector & real_projection, Vector & natural_projection) { const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(this->positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); #define GET_NATURAL_COORDINATE(type) \ ElementClass::inverseMap(real_projection, nodes_coord, natural_projection) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE); #undef GET_NATURAL_COORDINATE } /* -------------------------------------------------------------------------- */ void ContactDetector::computeTangentsOnElement(const Element & el, Vector & projection, Matrix & tangents) { const ElementType & type = el.type; UInt nb_nodes_master = Mesh::getNbNodesPerElement(type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shapes_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix coords(spatial_dimension, nb_nodes_master); coordinatesOfElement(el, coords); tangents.mul(shapes_derivatives, coords); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension -1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); } } // akantu diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index bb333fd51..ef681da0f 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,641 +1,657 @@ /** * @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" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, Array & positions, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id, memory_id), current_positions(positions) { 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, current_positions, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : ContactMechanicsModel(mesh, mesh.getNodes(), dim, id, memory_id, dof_manager, model_type) {} /* -------------------------------------------------------------------------- */ 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() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { auto & dof_manager = this->getDOFManager(); /* -------------------------------------------------------------------------- */ // 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->gaps, 1, "gaps"); this->allocNodalField(this->previous_gaps, 1, "previous_gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->normals, spatial_dimension, "normals"); this->allocNodalField(this->tangents, spatial_dimension, "tangents"); /* -------------------------------------------------------------------------- */ // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", _tsst_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", _tsst_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", _tsst_dynamic_lumped); break; } default: return std::make_tuple("unkown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson_contact; options.integration_scheme_type["displacement"] = _ist_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(); + internal_force->resize(nb_nodes, 0.); + // assemble the forces due to contact auto assemble = [&](auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); assemble(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { //this->contact_map.clear(); this->detector->search(this->contact_map); for (auto & entry : contact_map) { auto & element = entry.second; if (element.gap < 0) element.gap = std::abs(element.gap); else element.gap = -element.gap; } this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search(Array & increment) { this->contact_map.clear(); this->detector->search(this->contact_map); for (auto & entry : contact_map) { auto & element = entry.second; const auto & connectivity = element.connectivity; const auto & normal = element.normal; auto slave_node = connectivity[0]; Vector u_slave(spatial_dimension); auto master_node = connectivity[1]; Vector u_master(spatial_dimension); for (UInt s : arange(spatial_dimension)) { u_slave(s) = increment(slave_node, s); u_master(s) = increment(master_node, s); } // todo check this initial guess auto u = (u_master.norm() > u_slave.norm()) ? u_master : u_slave; Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension); if (element.gap - uv <= 0) { element.gap = std::abs(element.gap - uv); } else { element.gap = 0.0; } } this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleFieldsFromContactMap() { + UInt nb_nodes = mesh.getNbNodes(); + this->previous_gaps->copy(*(this->gaps)); this->gaps->clear(); + + gaps->resize(nb_nodes, 0.); + normals->resize(nb_nodes, 0.); + tangents->resize(nb_nodes, 0.); if (this->contact_map.empty()) return; for (auto & entry : contact_map) { const auto & element = entry.second; auto connectivity = element.connectivity; auto node = connectivity(0); (*gaps)[node] = element.gap; for (UInt i = 0; i < spatial_dimension; ++i) { (*normals)(node, i) = element.normal[i]; (*tangents)(node, i) = element.tangents(0, i); } } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas() { + UInt nb_nodes = mesh.getNbNodes(); + this->nodal_area->clear(); this->external_force->clear(); + nodal_area->resize(nb_nodes, 0.); + external_force->resize(nb_nodes, 0.); + + std::cerr << *external_force <applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), this->detector->getSurfaceId()); this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), this->detector->getSurfaceId()); for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); auto & force = std::get<1>(tuple); for (auto & f : force) area += pow(f, 2); area = sqrt(area); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep() {} /* -------------------------------------------------------------------------- */ 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(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * ContactMechanicsModel::createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { dumper::Field * field = nullptr; return field; } /* -------------------------------------------------------------------------- */ dumper::Field * 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["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["previous_gaps"] = this->previous_gaps; real_nodal_fields["areas"] = this->nodal_area; dumper::Field * field = nullptr; 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 /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { 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 451b072fb..2a6d8b9ad 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,333 +1,333 @@ /** * @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 ModelType model_type = ModelType::_contact_mechanics_model); ContactMechanicsModel( Mesh & mesh, Array & positions, UInt spatial_dimension = _all_dimensions, const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, 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() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: void search(); void search(Array &); void computeNodalAreas(); void assembleFieldsFromContactMap(); /* ------------------------------------------------------------------------ */ /* 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(); - + /* ------------------------------------------------------------------------ */ /* 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: dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * 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 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::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(CurrentPositions, current_positions, Array &); /// get the contat map inline std::map & getContactMap() { return contact_map; } /// inline void setPositions(Array positions) { detector->setPositions(positions); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array Array * displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// contact forces array Array * internal_force{nullptr}; /// external forces array Array * external_force{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; /// array to store gap between slave and master Array * gaps{nullptr}; /// array to store gap from previous iteration Array * previous_gaps{nullptr}; /// array to store normals from master to slave Array * normals{nullptr}; /// array to store tangents on the master element Array * tangents{nullptr}; /// array to store nodal areas Array * nodal_area{nullptr}; /// array of current position used during update residual Array & current_positions; /// 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; /// mapping between slave node its respective contact element std::map contact_map; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc new file mode 100644 index 000000000..26bd28e96 --- /dev/null +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -0,0 +1,570 @@ +/** + * @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 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); + contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, + "contact_mechanics_model", 0, 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); +} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidCohesiveContact::initModel() { + + + getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); + getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); + + getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); + getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); + + //getFEEngine().initShapeFunctions(_not_ghost); + //getFEEngine().initShapeFunctions(_ghost); + +} + +/* -------------------------------------------------------------------------- */ +FEEngine & CouplerSolidCohesiveContact::getFEEngineBoundary(const ID & name) { + return dynamic_cast( + getFEEngineClassBoundary(name)); +} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidCohesiveContact::initSolver(TimeStepSolverType, NonLinearSolverType) { + DOFManager & dof_manager = this->getDOFManager(); +} + +/* -------------------------------------------------------------------------- */ +std::tuple +CouplerSolidCohesiveContact::getDefaultSolverID(const AnalysisMethod & method) { + + switch (method) { + case _explicit_contact: { + return std::make_tuple("explicit_contact", _tsst_static); + } + case _implicit_contact: { + return std::make_tuple("implicit_contact", _tsst_static); + } + case _explicit_dynamic_contact: { + return std::make_tuple("explicit_dynamic_contact", _tsst_dynamic_lumped); + break; + } + default: + return std::make_tuple("unkown", _tsst_not_defined); + } +} + +/* -------------------------------------------------------------------------- */ +TimeStepSolverType CouplerSolidCohesiveContact::getDefaultSolverType() const { + return _tsst_dynamic_lumped; +} + + +/* -------------------------------------------------------------------------- */ +ModelSolverOptions CouplerSolidCohesiveContact::getDefaultSolverOptions( + const TimeStepSolverType & type) const { + ModelSolverOptions options; + + switch (type) { + case _tsst_dynamic_lumped: { + options.non_linear_solver_type = _nls_lumped; + options.integration_scheme_type["displacement"] = _ist_central_difference; + options.solution_type["displacement"] = IntegrationScheme::_acceleration; + break; + } + case _tsst_dynamic: { + options.non_linear_solver_type = _nls_lumped; + options.integration_scheme_type["displacement"] = _ist_central_difference; + options.solution_type["displacement"] = IntegrationScheme::_acceleration; + break; + } + case _tsst_static: { + options.non_linear_solver_type = _nls_newton_raphson_contact; + options.integration_scheme_type["displacement"] = _ist_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(); + auto & contact_map = contact->getContactMap(); + + switch (method) { + case _explicit_dynamic_contact: + case _explicit_contact: { + for (auto & pair: contact_map) { + auto & connectivity = pair.second.connectivity; + for (auto node : connectivity) { + for (auto s : arange(spatial_dimension)) + external_force(node, s) = contact_force(node, s); + } + } + break; + } + 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 & contact_map = contact->getContactMap(); + + + switch (method) { + case _explicit_dynamic_contact: + case _explicit_contact: { + for (auto & pair: contact_map) { + auto & connectivity = pair.second.connectivity; + for (auto node : connectivity) { + for (auto s : arange(spatial_dimension)) + external_force(node, s) = contact_force(node, s); + } + } + break; + } + 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; + } + + AKANTU_CUSTOM_EXCEPTION( + debug::SolverCallbackResidualPartUnknown(residual_part)); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidCohesiveContact::beforeSolveStep() {} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidCohesiveContact::afterSolveStep() {} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidCohesiveContact::predictor() { + + switch (method) { + case _explicit_dynamic_contact: { + Array displacement(0, Model::spatial_dimension); + + Array current_positions(0, Model::spatial_dimension); + auto positions = mesh.getNodes(); + current_positions.copy(positions); + + auto us = this->getDOFManager().getDOFs("displacement"); + //const auto deltas = this->getDOFManager().getSolution("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->setPositions(current_positions); + contact->search(); + break; + } + default: + break; + } + + +} + +/* -------------------------------------------------------------------------- */ +void CouplerSolidCohesiveContact::corrector() { + + switch (method) { + case _implicit_contact: + case _explicit_contact: { + Array displacement(0, Model::spatial_dimension); + + Array current_positions(0, Model::spatial_dimension); + auto positions = mesh.getNodes(); + current_positions.copy(positions); + + auto us = this->getDOFManager().getDOFs("displacement"); + const auto deltas = this->getDOFManager().getSolution("displacement"); + const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); + + for (auto && tuple : zip(make_view(us), + deltas, + make_view(blocked_dofs), + make_view(current_positions))) { + auto & u = std::get<0>(tuple); + const auto & delta = std::get<1>(tuple); + const auto & bld = std::get<2>(tuple); + auto & cp = std::get<3>(tuple); + + if (not bld) + cp += u + delta; + } + + contact->setPositions(current_positions); + contact->search(); + + break; + } + default: + break; + } + + /*auto & internal_force = solid->getInternalForce(); + auto & external_force = solid->getExternalForce(); + + + + std::stringstream filename; + filename << "out" << "-00" << step << ".csv"; + + std::ofstream outfile(filename.str()); + + outfile << "x,gap,residual" << std::endl; + + auto & contact_map = contact->getContactMap(); + for (auto & pair: contact_map) { + auto & connectivity = pair.second.connectivity; + auto node = connectivity(0); + if (pair.second.gap > 0) { + outfile << positions(node, 0) << "," << pair.second.gap << "," + << external_force(node, 1) + internal_force(node, 1) << std::endl; + } + } + + outfile.close(); + step++;*/ +} + +/* -------------------------------------------------------------------------- */ +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 _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 + +/* -------------------------------------------------------------------------- */ +dumper::Field * CouplerSolidCohesiveContact::createElementalField( + const std::string & field_name, const std::string & group_name, + bool padding_flag, const UInt & spatial_dimension, + const ElementKind & kind) { + + dumper::Field * field = nullptr; + + field = solid->createElementalField(field_name, group_name, padding_flag, + spatial_dimension, kind); + + return field; +} + +/* -------------------------------------------------------------------------- */ +dumper::Field * +CouplerSolidCohesiveContact::createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) { + + dumper::Field * field = nullptr; + + if (field_name == "contact_force" or field_name == "normals" 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; +} + +/* -------------------------------------------------------------------------- */ +dumper::Field * CouplerSolidCohesiveContact::createNodalFieldBool( + const std::string & field_name, const std::string & group_name, + __attribute__((unused)) bool padding_flag) { + + dumper::Field * field = nullptr; + + field = solid->createNodalFieldBool(field_name, group_name, padding_flag); + return field; +} + +#else + +/* -------------------------------------------------------------------------- */ +dumper::Field * CouplerSolidCohesiveContact::createElementalField(const std::string &, + const std::string &, + bool, const UInt &, + const ElementKind &) { + return nullptr; +} + +/* ----------------------------------------------------------------------- */ +dumper::Field * CouplerSolidCohesiveContact::createNodalFieldReal(const std::string &, + const std::string &, + bool) { + return nullptr; +} + +/*-------------------------------------------------------------------*/ +dumper::Field * 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/contact_mechanics/contact_mechanics_model.hh b/src/model/model_couplers/coupler_solid_cohesive_contact.hh similarity index 54% copy from src/model/contact_mechanics/contact_mechanics_model.hh copy to src/model/model_couplers/coupler_solid_cohesive_contact.hh index 451b072fb..3840763e6 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.hh @@ -1,333 +1,289 @@ /** - * @file contact_mechanics_model.hh + * @file coupler_solid_cohesive_contact.hh * * @author Mohit Pundir * - * @date creation: Tue Jul 27 2010 - * @date last modification: Wed Feb 21 2018 + * @date creation: Thu Jan 17 2019 + * @date last modification: Thu Jan 17 2019 * - * @brief Model of Contact Mechanics + * @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_detector.hh" +#include "contact_mechanics_model.hh" #include "data_accessor.hh" -#include "fe_engine.hh" #include "model.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "sparse_matrix.hh" +#include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ -#ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ -#define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ +#ifndef __AKANTU_COUPLER_SOLID_COHESIVE_CONTACT_HH__ +#define __AKANTU_COUPLER_SOLID_COHESIVE_CONTACT_HH__ +/* ------------------------------------------------------------------------ */ +/* Coupling : Solid Mechanics / Contact Mechanics */ +/* ------------------------------------------------------------------------ */ namespace akantu { -class Resolution; template class IntegratorGauss; template class ShapeLagrange; +class DOFManager; } // namespace akantu /* -------------------------------------------------------------------------- */ + namespace akantu { /* -------------------------------------------------------------------------- */ -class ContactMechanicsModel : public Model, - public DataAccessor, - public DataAccessor, - public BoundaryCondition { +class CouplerSolidCohesiveContact : public Model, + public DataAccessor, + public DataAccessor, + public BoundaryCondition { /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ + /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ - using MyFEEngineType = FEEngineTemplate; + using MyFEEngineCohesiveType = + FEEngineTemplate; + using MyFEEngineFacetType = + 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 ModelType model_type = ModelType::_contact_mechanics_model); - - ContactMechanicsModel( - Mesh & mesh, Array & positions, + CouplerSolidCohesiveContact(Mesh & mesh, UInt spatial_dimension = _all_dimensions, - const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, + const ID & id = "coupler_solid_cohesive_contact", std::shared_ptr dof_manager = nullptr, - const ModelType model_type = ModelType::_contact_mechanics_model); + const ModelType model_type = ModelType::_coupler_solid_contact); - ~ContactMechanicsModel() override; + ~CouplerSolidCohesiveContact() 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; + + /// get some default values for derived classes + std::tuple + getDefaultSolverID(const AnalysisMethod & method) override; - /// call back for the solver, computes the force residual + + /* ------------------------------------------------------------------------ */ + /* 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; - /// get some default values for derived classes - std::tuple - getDefaultSolverID(const AnalysisMethod & method) override; + /// callback for the solver, this is called at beginning of solve + void predictor() override; - ModelSolverOptions - getDefaultSolverOptions(const TimeStepSolverType & type) const; + /// 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() override; - - /// function to print the containt of the class - void printself(std::ostream & stream, int indent = 0) const override; + + /// callback for the model to instantiate the matricess when needed + void initSolver(TimeStepSolverType, NonLinearSolverType) override; /* ------------------------------------------------------------------------ */ - /* Contact Detection */ + /* Mass matrix for solid mechanics model */ /* ------------------------------------------------------------------------ */ public: - void search(); + /// assemble the lumped mass matrix + void assembleMassLumped(); - void search(Array &); + /// assemble the mass matrix for consistent mass resolutions + void assembleMass(); - void computeNodalAreas(); - - void assembleFieldsFromContactMap(); +protected: + /// assemble the lumped mass matrix for local and ghost elements + void assembleMassLumped(GhostType ghost_type); - /* ------------------------------------------------------------------------ */ - /* 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); + /// assemble the mass matrix for either _ghost or _not_ghost elements + void assembleMass(GhostType ghost_type); + protected: - /// register a resolution in the dynamic database - Resolution & registerNewResolution(const ParserSection & res_section); + /* -------------------------------------------------------------------------- */ + TimeStepSolverType getDefaultSolverType() const override; + /* -------------------------------------------------------------------------- */ + ModelSolverOptions + getDefaultSolverOptions(const TimeStepSolverType & type) const; - /// read the resolution files to instantiate all the resolutions - void instantiateResolutions(); - /* ------------------------------------------------------------------------ */ - /* Solver Interface */ +public: + bool isDefaultSolverExplicit() { + return method == _explicit_dynamic_contact; + } + + /* ------------------------------------------------------------------------ */ public: - /// assembles the contact stiffness matrix - virtual void assembleStiffnessMatrix(); - - /// assembles the contant internal forces - virtual void assembleInternalForces(); + // 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; - /* ------------------------------------------------------------------------ */ - /* Dumpable interface */ - /* ------------------------------------------------------------------------ */ -public: - dumper::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * 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 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 solid mechanics model + AKANTU_GET_MACRO(SolidMechanicsModelCohesive, *solid, SolidMechanicsModelCohesive &); - /// get the ContactMechanics::areas (nodal areas) - AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); - - /// get the ContactMechanics::internal_force vector (internal forces) - AKANTU_GET_MACRO(CurrentPositions, current_positions, Array &); - - /// get the contat map - inline std::map & getContactMap() { - return contact_map; - } - - /// - inline void setPositions(Array positions) { - detector->setPositions(positions); - } + /// get the contact mechanics model + AKANTU_GET_MACRO(ContactMechanicsModel, *contact, ContactMechanicsModel &); + /* ------------------------------------------------------------------------ */ - /* Class Members */ + /* Dumpable interface */ /* ------------------------------------------------------------------------ */ -private: - /// tells if the resolutions are instantiated - bool are_resolutions_instantiated; - - /// displacements array - Array * displacement{nullptr}; +public: + + dumper::Field * createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; - /// increment of displacement - Array * displacement_increment{nullptr}; + dumper::Field * createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; - /// contact forces array - Array * internal_force{nullptr}; + dumper::Field * createElementalField(const std::string & field_name, + const std::string & group_name, + bool padding_flag, + const UInt & spatial_dimension, + const ElementKind & kind) override; - /// external forces array - Array * external_force{nullptr}; + virtual void dump(const std::string & dumper_name); - /// boundary vector - Array * blocked_dofs{nullptr}; + virtual void dump(const std::string & dumper_name, UInt step); - /// array to store gap between slave and master - Array * gaps{nullptr}; + virtual void dump(const std::string & dumper_name, Real time, UInt step); - /// array to store gap from previous iteration - Array * previous_gaps{nullptr}; + void dump() override; - /// array to store normals from master to slave - Array * normals{nullptr}; + virtual void dump(UInt step); - /// array to store tangents on the master element - Array * tangents{nullptr}; + virtual void dump(Real time, UInt step); - /// array to store nodal areas - Array * nodal_area{nullptr}; + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ +private: + /// solid mechanics model + SolidMechanicsModelCohesive * solid{nullptr}; - /// array of current position used during update residual - Array & current_positions; + /// contact mechanics model + ContactMechanicsModel * contact{nullptr}; - /// contact detection - std::unique_ptr detector; + /// + Array * displacement{nullptr}; - /// list of contact resolutions - std::vector> resolutions; + /// + Array * displacement_increment{nullptr}; + + /// external forces array + Array * external_force{nullptr}; - /// mapping between resolution name and resolution internal id - std::map resolutions_names_to_id; + bool search_for_contact; - /// mapping between slave node its respective contact element - std::map contact_map; + UInt step; }; } // namespace akantu -/* ------------------------------------------------------------------------ */ -/* inline functions */ -/* ------------------------------------------------------------------------ */ -#include "parser.hh" -#include "resolution.hh" -/* ------------------------------------------------------------------------ */ - -#endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ +#endif /* __COUPLER_SOLID_CONTACT_HH__ */