diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc index 73cb90e31..cd1da940a 100644 --- a/src/mesh/element_group.cc +++ b/src/mesh/element_group.cc @@ -1,227 +1,243 @@ /** * @file element_group.cc * * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Dec 09 2020 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "aka_csr.hh" #include "dumpable.hh" #include "dumpable_inline_impl.hh" #include "group_manager.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_utils.hh" #if defined(AKANTU_COHESIVE_ELEMENT) #include "cohesive_element_inserter.hh" #endif #include #include #include /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" namespace akantu { /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const std::string & group_name, const Mesh & mesh, NodeGroup & node_group, UInt dimension, const std::string & id) : mesh(mesh), name(group_name), elements("elements", id), node_group(node_group), dimension(dimension) { AKANTU_DEBUG_IN(); this->registerDumper("paraview_" + group_name, group_name, true); this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), _all_dimensions); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const ElementGroup & /*other*/) = default; /* -------------------------------------------------------------------------- */ void ElementGroup::clear() { elements.free(); } /* -------------------------------------------------------------------------- */ void ElementGroup::clear(ElementType type, GhostType ghost_type) { if (elements.exists(type, ghost_type)) { elements(type, ghost_type).clear(); } } /* -------------------------------------------------------------------------- */ bool ElementGroup::empty() const { return elements.empty(); } /* -------------------------------------------------------------------------- */ void ElementGroup::append(const ElementGroup & other_group) { AKANTU_DEBUG_IN(); node_group.append(other_group.node_group); /// loop on all element types in all dimensions for (auto ghost_type : ghost_types) { for (auto type : other_group.elementTypes(_ghost_type = ghost_type, _element_kind = _ek_not_defined)) { const Array & other_elem_list = other_group.elements(type, ghost_type); UInt nb_other_elem = other_elem_list.size(); Array * elem_list; UInt nb_elem = 0; /// create current type if doesn't exists, otherwise get information if (elements.exists(type, ghost_type)) { elem_list = &elements(type, ghost_type); nb_elem = elem_list->size(); } else { elem_list = &(elements.alloc(0, 1, type, ghost_type)); } /// append new elements to current list elem_list->resize(nb_elem + nb_other_elem); std::copy(other_elem_list.begin(), other_elem_list.end(), elem_list->begin() + nb_elem); } } this->optimize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) { ; } stream << space << "ElementGroup [" << std::endl; stream << space << " + name: " << name << std::endl; stream << space << " + dimension: " << dimension << std::endl; elements.printself(stream, indent + 1); node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void ElementGroup::optimize() { // increasing the locality of data when iterating on the element of a group for (auto ghost_type : ghost_types) { for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { auto & els = elements(type, ghost_type); std::sort(els.begin(), els.end()); auto end = std::unique(els.begin(), els.end()); els.resize(end - els.begin()); } } node_group.optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::fillFromNodeGroup() { CSR node_to_elem; MeshUtils::buildNode2Elements(this->mesh, node_to_elem, this->dimension); std::set seen; for (auto node : node_group) { auto ite = node_to_elem.begin(node); auto ende = node_to_elem.end(node); for (auto && element : node_to_elem.getRow(node)) { if (this->dimension != _all_dimensions && this->dimension != Mesh::getSpatialDimension(element.type)) { continue; } if (seen.find(element) != seen.end()) { continue; } auto nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); auto conn = mesh.getConnectivity(element); UInt count = 0; for (auto && n : conn) { count += (this->node_group.getNodes().find(n) != UInt(-1) ? 1 : 0); } if (count == nb_nodes_per_element) { this->add(element); } seen.insert(element); } } this->optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::addDimension(UInt dimension) { this->dimension = std::max(dimension, this->dimension); } /* -------------------------------------------------------------------------- */ void ElementGroup::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { #if defined(AKANTU_COHESIVE_ELEMENT) if (aka::is_of_type(event)) { // nodes might have changed in the connectivity node_group.clear(); + const auto & mesh_to_mesh_facet = + mesh.getData("mesh_to_mesh_facet"); + for (auto ghost_type : ghost_types) { for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { auto & els = elements(type, ghost_type); + + const auto & mesh_to_mesh_facet_type = + mesh_to_mesh_facet(type, ghost_type); + auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto && conn_it = make_view(mesh.getConnectivity(type, ghost_type), nb_nodes_per_element) .begin(); + + auto && mesh_facet_conn_it = + make_view(mesh.getMeshFacets().getConnectivity(type, ghost_type), + nb_nodes_per_element) + .begin(); + for (auto element : els) { + auto && mesh_facet_conn = + mesh_facet_conn_it[mesh_to_mesh_facet_type(element).element]; auto && conn = conn_it[element]; + conn = mesh_facet_conn; for (auto && n : conn) { node_group.add(n, false); } } } } node_group.optimize(); } #endif } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index 61e4d3df0..1f34f3a85 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,658 +1,665 @@ /** * @file mesh.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Mohit Pundir * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 09 2021 * * @brief class handling meshes * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_global_data_updater.hh" #include "mesh_io.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" #include "periodic_node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "dumper_internal_material_field.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, Communicator & communicator) : GroupManager(*this, id + ":group_manager"), MeshData("mesh_data", id), id(id), connectivities("connectivities", id), ghosts_counters("ghosts_counters", id), normals("normals", id), spatial_dimension(spatial_dimension), size(spatial_dimension, 0.), bbox(spatial_dimension), bbox_local(spatial_dimension), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id) : Mesh(spatial_dimension, id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); this->nodes_flags = std::make_shared>(0, 1, NodeFlag::_normal, id + ":nodes_flags"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id) : Mesh(spatial_dimension, id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, ElementType type, GhostType ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } class FacetGlobalConnectivityAccessor : public DataAccessor { public: FacetGlobalConnectivityAccessor(Mesh & mesh) : global_connectivity("global_connectivity", "facet_connectivity_synchronizer") { global_connectivity.initialize( mesh, _spatial_dimension = _all_dimensions, _with_nb_element = true, _with_nb_nodes_per_element = true, _element_kind = _ek_regular); mesh.getGlobalConnectivity(global_connectivity); } UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override { UInt size = 0; if (tag == SynchronizationTag::_smmc_facets_conn) { UInt nb_nodes = Mesh::getNbNodesPerElementList(elements); size += nb_nodes * sizeof(UInt); } return size; } void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_smmc_facets_conn) { for (const auto & element : elements) { const auto & conns = global_connectivity(element.type, element.ghost_type); for (auto n : arange(conns.getNbComponent())) { buffer << conns(element.element, n); } } } } void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_smmc_facets_conn) { for (const auto & element : elements) { auto & conns = global_connectivity(element.type, element.ghost_type); for (auto n : arange(conns.getNbComponent())) { buffer >> conns(element.element, n); } } } } AKANTU_GET_MACRO(GlobalConnectivity, (global_connectivity), decltype(auto)); protected: ElementTypeMapArray global_connectivity; }; /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (mesh_facets) { AKANTU_DEBUG_OUT(); return *mesh_facets; } mesh_facets = std::make_unique(spatial_dimension, this->nodes, getID() + ":" + id); mesh_facets->mesh_parent = this; mesh_facets->is_mesh_facets = true; mesh_facets->nodes_flags = this->nodes_flags; mesh_facets->nodes_global_ids = this->nodes_global_ids; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); FacetGlobalConnectivityAccessor data_accessor(*mesh_facets); /// communicate mesh_facets->element_synchronizer->synchronizeOnce( data_accessor, SynchronizationTag::_smmc_facets_conn); /// flip facets MeshUtils::flipFacets(*mesh_facets, data_accessor.getGlobalConnectivity(), _ghost); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } auto & mesh_phys_data = this->getData("physical_names"); + auto & mesh_to_mesh_facet = this->getData("mesh_to_mesh_facet"); + mesh_to_mesh_facet.initialize(*this, + _spatial_dimension = spatial_dimension - 1, + _with_nb_element = true); + auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes(spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) { tolerance *= norm_barycenter; } Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view( barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name auto && facet_element = Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}; phys_data(facet_element) = mesh_phys_data(element); + + mesh_to_mesh_facet(element) = facet_element; }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { AKANTU_DEBUG_ASSERT(not is_distributed, "You cannot read a mesh that is already distributed"); MeshIO::read(filename, *this, mesh_io_type); auto types = this->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); auto it = types.begin(); auto last = types.end(); if (it == last) { AKANTU_DEBUG_WARNING( "The mesh contained in the file " << filename << " does not seem to be of the good dimension." << " No element of dimension " << spatial_dimension << " were read."); } this->makeReady(); } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO::write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::makeReady() { this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_flags->resize(nodes->size(), NodeFlag::_normal); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); bbox_local.reset(); for (auto & pos : make_view(*nodes, spatial_dimension)) { // if(!isPureGhostNode(i)) bbox_local += pos; } if (this->is_distributed) { bbox = bbox_local.allSum(*communicator); } else { bbox = bbox_local; } size = bbox.size(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_not_defined, _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) { continue; } auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [&](UInt l) -> UInt { return this->getNodeGlobalId(l); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name, const std::string & group_name) { if (group_name == "all") { return this->getDumper(dumper_name); } return element_groups[group_name]->getDumper(dumper_name); } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & arrays) { ElementTypeMap nb_data_per_elem; for (auto type : arrays.elementTypes(_element_kind = _ek_not_defined)) { UInt nb_elements = this->getNbElement(type); auto & array = arrays(type); nb_data_per_elem(type) = array.getNbComponent() * array.size(); nb_data_per_elem(type) /= nb_elements; } return nb_data_per_elem; } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array); /* -------------------------------------------------------------------------- */ template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind) { std::shared_ptr field; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind); template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind); /* -------------------------------------------------------------------------- */ void Mesh::distributeImpl( Communicator & communicator, const std::function & edge_weight_function [[gnu::unused]], const std::function & vertex_weight_function [[gnu::unused]]) { AKANTU_DEBUG_ASSERT(is_distributed == false, "This mesh is already distribute"); this->communicator = &communicator; this->element_synchronizer = std::make_unique( *this, this->getID() + ":element_synchronizer", true); this->node_synchronizer = std::make_unique( *this, this->getID() + ":node_synchronizer", true); Int psize = this->communicator->getNbProc(); if (psize > 1) { #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize, edge_weight_function, vertex_weight_function); MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition); } else { MeshUtilsDistribution::distributeMeshCentralized(*this, 0); } #else AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool"); #endif } // if (psize > 1) this->is_distributed = true; this->computeBoundingBox(); MeshIsDistributedEvent event(*this, AKANTU_CURRENT_FUNCTION); this->sendEvent(event); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) const { for (const auto & node : node_list) { for (const auto & element : *nodes_to_elements[node]) { elements.push_back(element); } } } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const UInt & node, Array & elements) const { for (const auto & element : *nodes_to_elements[node]) { elements.push_back(element); } } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements(UInt dimension) { Element e; UInt nb_nodes = nodes->size(); this->nodes_to_elements.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) { this->nodes_to_elements[n]->clear(); } else { this->nodes_to_elements[n] = std::make_unique>(); } } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); auto connectivity = connectivities(type, ghost_type); auto conn_it = connectivity.begin(connectivity.getNbComponent()); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (auto node : conn) { nodes_to_elements[node]->insert(e); } } } } } /* -------------------------------------------------------------------------- */ std::tuple Mesh::updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { if (global_data_updater) { return this->global_data_updater->updateData(nodes_event, elements_event); } return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } /* -------------------------------------------------------------------------- */ void Mesh::registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->global_data_updater = std::move(global_data_updater); } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this, "new_numbering", AKANTU_CURRENT_FUNCTION); auto & remove_list = event.getList(); auto & new_numbering = event.getNewNumbering(); for (auto && el : elements) { if (el.ghost_type != _not_ghost) { auto & count = ghosts_counters(el); --count; if (count > 0) { continue; } } remove_list.push_back(el); if (not new_numbering.exists(el.type, el.ghost_type)) { auto nb_element = mesh.getNbElement(el.type, el.ghost_type); auto & numbering = new_numbering.alloc(nb_element, 1, el.type, el.ghost_type); for (auto && pair : enumerate(numbering)) { std::get<1>(pair) = std::get<0>(pair); } } new_numbering(el) = UInt(-1); } auto find_last_not_deleted = [](auto && array, Int start) -> Int { do { --start; } while (start >= 0 and array[start] == UInt(-1)); return start; }; auto find_first_deleted = [](auto && array, Int start) -> Int { auto begin = array.begin(); auto it = std::find_if(begin + start, array.end(), [](auto & el) { return el == UInt(-1); }); return Int(it - begin); }; for (auto ghost_type : ghost_types) { for (auto type : new_numbering.elementTypes(_ghost_type = ghost_type)) { auto & numbering = new_numbering(type, ghost_type); auto last_not_delete = find_last_not_deleted(numbering, numbering.size()); if (last_not_delete < 0) { continue; } auto pos = find_first_deleted(numbering, 0); while (pos < last_not_delete) { std::swap(numbering[pos], numbering[last_not_delete]); last_not_delete = find_last_not_deleted(numbering, last_not_delete); pos = find_first_deleted(numbering, pos + 1); } } } this->sendEvent(event); } } // namespace akantu diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/CMakeLists.txt index 8d03d640b..6669d494b 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/CMakeLists.txt @@ -1,42 +1,50 @@ #=============================================================================== # @file CMakeLists.txt # # @author Fabian Barras # # @date creation: Sun Oct 19 2014 # @date last modification: Fri Jan 15 2016 # # @brief Tests insertion of cohesive elements # # # @section LICENSE # # Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License along # with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_mesh(3d_spherical_inclusion 3d_spherical_inclusion.geo 3 2) register_test(test_cohesive_insertion_along_physical_surfaces SOURCES test_cohesive_insertion_along_physical_surfaces.cc DEPENDS 3d_spherical_inclusion FILES_TO_COPY input_file.dat PACKAGE cohesive_element ) + +add_mesh(group_update_bar bar.geo DIM 2 ORDER 1) +register_test(test_cohesive_group_update + SOURCES test_group_update.cc + DEPENDS group_update_bar + FILES_TO_COPY material.dat + PACKAGE cohesive_element + ) diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/bar.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/bar.geo new file mode 100644 index 000000000..e126ac447 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/bar.geo @@ -0,0 +1,30 @@ +// Lenght of the bar (m) +L = 50e-3; +x0 = -L/2; +xf = L/2; + +// Number of triangular elements (n_el) +n_el = 2; + +// Length of each linear element (h) +h = L/(n_el/2); + + +Point(1) = { x0, 0, 0, h }; +Point(2) = { xf, 0, 0, h }; +Point(3) = { xf, h, 0, h }; +Point(4) = { x0, h, 0, h }; + +Line(1) = {1,2}; +Line(2) = {2,3}; +Line(3) = {3,4}; +Line(4) = {4,1}; + +Line Loop(1) = {1,2,3,4}; +Plane Surface(1) = {1}; +Physical Surface(1) = {1}; +Physical Line("YBlocked") = {1,3}; +Physical Line("right") = {2}; +Physical Line("left") = {4}; + +Transfinite Surface {1} Left; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/material.dat new file mode 100644 index 000000000..840bd7e5c --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/material.dat @@ -0,0 +1,19 @@ +seed = 1.0 +model solid_mechanics_model_cohesive [ + + material elastic [ + name = linear + rho = 2750.0 # Density (kg/m3) + E = 275.0e9 # Young's module (Pa) + nu = 0.3 + finite_deformation = true + ] + + material cohesive_linear [ + name = insertion + sigma_c = 300.0e6 uniform [-1e6, 1e6] # critical stress (Pa) + G_c = 100.0 # Fracture energy (N/m) + beta = 0. + penalty = 1e11 + ] +] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_group_update.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_group_update.cc new file mode 100644 index 000000000..daf419c65 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_group_update.cc @@ -0,0 +1,65 @@ +#include +#include + +using namespace akantu; + +int main(int argc, char * argv[]) { + initialize("material.dat", argc, argv); + + Int spatial_dimension = 2; + Mesh mesh(spatial_dimension); + mesh.read("bar.msh"); + + // Create model + SolidMechanicsModelCohesive model(mesh); + model.initFull(_analysis_method = _static, _is_extrinsic = true); + + // Configure solver + auto & solver = model.getNonLinearSolver("static"); + solver.set("max_iterations", 100); + solver.set("threshold", 1e-10); + solver.set("convergence_type", SolveConvergenceCriteria::_residual); + model.initNewSolver(_explicit_lumped_mass); + + // Dynamic insertion of cohesive elements + model.updateAutomaticInsertion(); + + auto dt_crit = model.getStableTimeStep(); // Critical time step + auto dt = dt_crit * 0.1; // Adopted time step + auto time_simulation = 6.0e-6; // Total time of simulation (s) + auto n_steps = int(time_simulation / dt); // Number of time steps + + // Apply Dirichlet BC to block displacements at y direction on top and + // bottom + model.applyBC(BC::Dirichlet::FixedValue(0., _y), "YBlocked"); + + // Constant velocity boundary condition + // Applied strain rate(s - 1) + auto strain_rate = 1e5; + // Applied velocity at the boundary + auto vel = + strain_rate * (mesh.getUpperBounds()(_x) - mesh.getLowerBounds()(_x)) / 2; + + model.applyBC(BC::Dirichlet::IncrementValue(-vel * dt, _x), "left"); + model.applyBC(BC::Dirichlet::IncrementValue(vel * dt, _x), "right"); + + // VTK plot setup + model.setBaseName("bar"); + model.addDumpFieldVector("displacement"); + // VTK plot setup for Cohesive model + model.setBaseNameToDumper("cohesive elements", "cohesive"); + model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); + + for (auto n : arange(n_steps)) { + // Apply velocity at the extremities + model.applyBC(BC::Dirichlet::IncrementValue(-vel * dt, _x), "left"); + model.applyBC(BC::Dirichlet::IncrementValue(vel * dt, _x), "right"); + + model.dump(); + model.dump("cohesive elements"); + + // Run simulation + model.checkCohesiveStress(); + model.solveStep("explicit_lumped"); + } +}