diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index c3bb62df9..bb2c6741f 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,661 +1,658 @@ /** * @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 * * * 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.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 /* -------------------------------------------------------------------------- */ #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, 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 & 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); }, _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); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER 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); #endif /* -------------------------------------------------------------------------- */ 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 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) 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); - Array::const_iterator> conn_it = - connectivities(type, ghost_type) - .begin(Mesh::getNbNodesPerElement(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 (UInt n = 0; n < conn.size(); ++n) { - nodes_to_elements[conn(n)]->insert(e); + 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/src/mesh_utils/cohesive_element_inserter_helper.cc b/src/mesh_utils/cohesive_element_inserter_helper.cc index bf1b923c2..01ad01b92 100644 --- a/src/mesh_utils/cohesive_element_inserter_helper.cc +++ b/src/mesh_utils/cohesive_element_inserter_helper.cc @@ -1,949 +1,949 @@ /** * @file cohesive_element_inserter_helper.cc * * @author Nicolas Richart * * @date creation jeu sep 03 2020 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 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 "cohesive_element_inserter_helper.hh" /* -------------------------------------------------------------------------- */ #include "data_accessor.hh" #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ CohesiveElementInserterHelper::CohesiveElementInserterHelper( Mesh & mesh, const ElementTypeMapArray & facet_insertion) : doubled_nodes(0, 2, "doubled_nodes"), mesh(mesh), mesh_facets(mesh.getMeshFacets()) { auto spatial_dimension = mesh_facets.getSpatialDimension(); for (auto gt : ghost_types) { for (auto type : mesh_facets.elementTypes(_ghost_type = gt)) { nb_new_facets(type, gt) = mesh_facets.getNbElement(type, gt); } } std::array nb_facet_to_insert{0, 0}; // creates a vector of the facets to insert std::vector potential_facets_to_double; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { const auto & f_insertion = facet_insertion(type_facet, gt_facet); auto & counter = nb_facet_to_insert[gt_facet == _not_ghost ? 0 : 1]; for (auto && data : enumerate(f_insertion)) { if (std::get<1>(data)) { UInt el = std::get<0>(data); potential_facets_to_double.push_back({type_facet, el, gt_facet}); ++counter; } } } } // Defines a global order of insertion oof cohesive element to ensure node // doubling appens in the smae order, this is necessary for the global node // numbering if (mesh_facets.isDistributed()) { const auto & comm = mesh_facets.getCommunicator(); ElementTypeMapArray global_orderings; global_orderings.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = -1); auto starting_index = nb_facet_to_insert[0]; comm.exclusiveScan(starting_index); // define the local numbers for facet to insert for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { for (auto data : zip(facet_insertion(type_facet, gt_facet), global_orderings(type_facet, gt_facet))) { if (std::get<0>(data)) { std::get<1>(data) = starting_index; ++starting_index; } } } } // retreives the oorder number from neighoring processors auto && synchronizer = mesh_facets.getElementSynchronizer(); SimpleElementDataAccessor data_accessor( global_orderings, SynchronizationTag::_ce_insertion_order); synchronizer.synchronizeOnce(data_accessor, SynchronizationTag::_ce_insertion_order); // sort the facets to double based on this global ordering std::sort(potential_facets_to_double.begin(), potential_facets_to_double.end(), [&global_orderings](auto && el1, auto && el2) { return global_orderings(el1) < global_orderings(el2); }); } for (auto d : arange(spatial_dimension)) { facets_to_double_by_dim[d] = std::make_unique>( 0, 2, "facets_to_double_" + std::to_string(d)); } auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; MeshAccessor mesh_accessor(mesh_facets); auto & elements_to_subelements = mesh_accessor.getElementToSubelement(); for (auto && facet_to_double : potential_facets_to_double) { auto gt_facet = facet_to_double.ghost_type; auto type_facet = facet_to_double.type; auto & elements_to_facets = elements_to_subelements(type_facet, gt_facet); auto & elements_to_facet = elements_to_facets(facet_to_double.element); if (elements_to_facet[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || elements_to_facet[1].kind() == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } auto new_facet = nb_new_facets(type_facet, gt_facet)++; facets_to_double.push_back(Vector{ facet_to_double, Element{type_facet, new_facet, gt_facet}}); /// update facet_to_element vector auto & element_to_update = elements_to_facet[1]; auto el = element_to_update.element; const auto & facets_to_elements = mesh_facets.getSubelementToElement( element_to_update.type, element_to_update.ghost_type); auto facets_to_element = Vector( make_view(facets_to_elements, facets_to_elements.getNbComponent()) .begin()[el]); auto facet_to_update = std::find(facets_to_element.begin(), facets_to_element.end(), facet_to_double); AKANTU_DEBUG_ASSERT(facet_to_update != facets_to_element.end(), "Facet not found"); facet_to_update->element = new_facet; /// update elements connected to facet const auto & first_facet_list = elements_to_facet; elements_to_facets.push_back(first_facet_list); /// set new and original facets as boundary facets elements_to_facets(new_facet)[0] = elements_to_facets(new_facet)[1]; elements_to_facets(new_facet)[1] = ElementNull; elements_to_facets(facet_to_double.element)[1] = ElementNull; } } /* -------------------------------------------------------------------------- */ inline auto CohesiveElementInserterHelper::hasElement(const Vector & nodes_element, const Vector & nodes) -> bool { // if one of the nodes of "nodes" is not in "nodes_element" it stops auto it = std::mismatch(nodes.begin(), nodes.end(), nodes_element.begin(), [&](auto && node, auto && /*node2*/) -> bool { auto it = std::find(nodes_element.begin(), nodes_element.end(), node); return (it != nodes_element.end()); }); // true if all "nodes" where found in "nodes_element" return (it.first == nodes.end()); } /* -------------------------------------------------------------------------- */ inline auto CohesiveElementInserterHelper::removeElementsInVector( const std::vector & elem_to_remove, std::vector & elem_list) -> bool { if (elem_list.size() <= elem_to_remove.size()) { return true; } auto el_it = elem_to_remove.begin(); auto el_last = elem_to_remove.end(); std::vector::iterator el_del; UInt deletions = 0; for (; el_it != el_last; ++el_it) { el_del = std::find(elem_list.begin(), elem_list.end(), *el_it); if (el_del != elem_list.end()) { elem_list.erase(el_del); ++deletions; } } AKANTU_DEBUG_ASSERT(deletions == 0 || deletions == elem_to_remove.size(), "Not all elements have been erased"); return deletions == 0; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, const std::vector * facet_list) { AKANTU_DEBUG_IN(); for (const auto & element : element_list) { if (element.type == _not_defined) { continue; } Vector connectivity = mesh.getConnectivity(element); if (element.kind() == _ek_cohesive) { AKANTU_DEBUG_ASSERT( facet_list != nullptr, "Provide a facet list in order to update cohesive elements"); const Vector facets = mesh_facets.getSubelementToElement(element); auto facet_nb_nodes = connectivity.size() / 2; /// loop over cohesive element's facets for (const auto & facet : enumerate(facets)) { /// skip facets if not present in the list if (std::find(facet_list->begin(), facet_list->end(), std::get<1>(facet)) == facet_list->end()) { continue; } auto n = std::get<0>(facet); auto begin = connectivity.begin() + static_cast(n * facet_nb_nodes); auto end = begin + facet_nb_nodes; auto it = std::find(begin, end, old_node); AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element"); *it = new_node; } } else { auto it = std::find(connectivity.begin(), connectivity.end(), old_node); AKANTU_DEBUG_ASSERT(it != connectivity.end(), "Node not found in current element"); /// update connectivity *it = new_node; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateSubelementToElement(UInt dim, bool facet_mode) { auto & facets_to_double = *facets_to_double_by_dim[dim]; auto & facets_to_subfacets = elementsOfDimToElementsOfDim( dim + static_cast(facet_mode), dim); for (auto && data : zip(make_view(facets_to_double, 2), facets_to_subfacets)) { const auto & old_subfacet = std::get<0>(data)[0]; const auto & new_subfacet = std::get<0>(data)[1]; auto & facet_to_subfacets = std::get<1>(data); MeshAccessor mesh_accessor(mesh_facets); for (auto & facet : facet_to_subfacets) { Vector subfacets = mesh_accessor.getSubelementToElement(facet); auto && subfacet_to_update_it = std::find(subfacets.begin(), subfacets.end(), old_subfacet); AKANTU_DEBUG_ASSERT(subfacet_to_update_it != subfacets.end(), "Subfacet not found"); *subfacet_to_update_it = new_subfacet; } } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateElementToSubelement(UInt dim, bool facet_mode) { auto & facets_to_double = *facets_to_double_by_dim[dim]; auto & facets_to_subfacets = elementsOfDimToElementsOfDim( dim + static_cast(facet_mode), dim); MeshAccessor mesh_accessor(mesh_facets); // resize the arrays mesh_accessor.getElementToSubelement().initialize( mesh_facets, _spatial_dimension = dim, _with_nb_element = true); for (auto && data : zip(make_view(facets_to_double, 2), facets_to_subfacets)) { const auto & new_facet = std::get<0>(data)[1]; mesh_accessor.getElementToSubelement(new_facet) = std::get<1>(data); } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateQuadraticSegments(UInt dim) { AKANTU_DEBUG_IN(); auto spatial_dimension = mesh.getSpatialDimension(); auto & facets_to_double = *facets_to_double_by_dim[dim]; MeshAccessor mesh_accessor(mesh_facets); auto & connectivities = mesh_accessor.getConnectivities(); /// this ones matter only for segments in 3D Array> * element_to_subfacet_double = nullptr; Array> * facet_to_subfacet_double = nullptr; if (dim == spatial_dimension - 2) { element_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 2, dim); facet_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 1, dim); } const auto & element_to_subelement = mesh_facets.getElementToSubelement(); std::vector middle_nodes; for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double[0]; if (old_facet.type != _segment_3) { continue; } auto node = connectivities(old_facet, 2); if (not mesh.isPureGhostNode(node)) { middle_nodes.push_back(node); } } auto n = doubled_nodes.size(); doubleNodes(middle_nodes); for (auto && data : enumerate(make_view(facets_to_double, 2))) { auto facet = std::get<0>(data); auto & old_facet = std::get<1>(data)[0]; if (old_facet.type != _segment_3) { continue; } auto old_node = connectivities(old_facet, 2); if (mesh.isPureGhostNode(old_node)) { continue; } auto new_node = doubled_nodes(n, 1); auto & new_facet = std::get<1>(data)[1]; connectivities(new_facet, 2) = new_node; if (dim == spatial_dimension - 2) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_subelement(new_facet, 0)); updateElementalConnectivity(mesh, old_node, new_node, (*element_to_subfacet_double)(facet), &(*facet_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_subelement(new_facet, 0)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserterHelper::insertCohesiveElement() { auto nb_new_facets = insertFacetsOnly(); if (nb_new_facets == 0) { return 0; } updateCohesiveData(); return nb_new_facets; } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserterHelper::insertFacetsOnly() { UInt spatial_dimension = mesh.getSpatialDimension(); if (facets_to_double_by_dim[spatial_dimension - 1]->empty()) { return 0; } if (spatial_dimension == 1) { doublePointFacet(); } else if (spatial_dimension == 2) { doubleFacets<1>(); findSubfacetToDouble<1>(); doubleSubfacet<2>(); } else if (spatial_dimension == 3) { doubleFacets<2>(); findSubfacetToDouble<2>(); doubleFacets<1>(); findSubfacetToDouble<1>(); doubleSubfacet<3>(); } return facets_to_double_by_dim[spatial_dimension - 1]->size(); } /* -------------------------------------------------------------------------- */ template void CohesiveElementInserterHelper::doubleFacets() { AKANTU_DEBUG_IN(); NewElementsEvent new_facets; auto spatial_dimension = mesh_facets.getSpatialDimension(); auto & facets_to_double = *facets_to_double_by_dim[dim]; MeshAccessor mesh_accessor(mesh_facets); for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto && old_facet = facet_to_double[0]; auto && new_facet = facet_to_double[1]; auto & facets_connectivities = mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type); facets_connectivities.resize( nb_new_facets(old_facet.type, old_facet.ghost_type)); auto facets_connectivities_begin = make_view(facets_connectivities, facets_connectivities.getNbComponent()) .begin(); // copy the connectivities Vector new_conn(facets_connectivities_begin[new_facet.element]); Vector old_conn(facets_connectivities_begin[old_facet.element]); new_conn = old_conn; // this will fail if multiple facet types exists for a given element type // \TODO handle multiple sub-facet types auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(old_facet.type); auto & subfacets_to_facets = mesh_accessor.getSubelementToElementNC( old_facet.type, old_facet.ghost_type); subfacets_to_facets.resize( nb_new_facets(old_facet.type, old_facet.ghost_type), ElementNull); auto subfacets_to_facets_begin = make_view(subfacets_to_facets, nb_subfacet_per_facet).begin(); // copy the subfacet to facets information Vector old_subfacets_to_facet( subfacets_to_facets_begin[old_facet.element]); Vector new_subfacet_to_facet( subfacets_to_facets_begin[new_facet.element]); new_subfacet_to_facet = old_subfacets_to_facet; for (auto & subfacet : old_subfacets_to_facet) { if (subfacet == ElementNull) { continue; } /// update facet_to_subfacet array mesh_accessor.getElementToSubelement(subfacet).push_back(new_facet); } new_facets.getList().push_back(new_facet); } /// update facet_to_subfacet and _segment_3 facets if any if (dim != spatial_dimension - 1) { updateSubelementToElement(dim, true); updateElementToSubelement(dim, true); updateQuadraticSegments(dim); } else { updateQuadraticSegments(dim); } mesh_facets.sendEvent(new_facets); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void CohesiveElementInserterHelper::findSubfacetToDouble() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); MeshAccessor mesh_accessor(mesh_facets); auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; auto & subfacets_to_facets = mesh_facets.getSubelementToElement(); auto & elements_to_facets = mesh_accessor.getElementToSubelement(); /// loop on every facet for (auto f : arange(2)) { for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto old_facet = facet_to_double[f]; auto new_facet = facet_to_double[1 - f]; const auto & starting_element = elements_to_facets(new_facet, 0)[0]; auto current_facet = old_facet; Vector subfacets_to_facet = subfacets_to_facets.get(old_facet); /// loop on every subfacet for (auto & subfacet : subfacets_to_facet) { if (subfacet == ElementNull) { continue; } if (dim == spatial_dimension - 2) { Vector subsubfacets_to_subfacet = subfacets_to_facets.get(subfacet); /// loop on every subsubfacet for (auto & subsubfacet : subsubfacets_to_subfacet) { if (subsubfacet == ElementNull) { continue; } Vector subsubfacet_connectivity = mesh_facets.getConnectivity(subsubfacet); std::vector element_list; std::vector facet_list; std::vector subfacet_list; /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet( starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, &subfacet_list) == false and removeElementsInVector( subfacet_list, elements_to_facets(subsubfacet)) == false) { Element new_subsubfacet{ subsubfacet.type, nb_new_facets(subsubfacet.type, subsubfacet.ghost_type)++, subsubfacet.ghost_type}; facets_to_double_by_dim[dim - 1]->push_back( Vector{subsubfacet, new_subsubfacet}); elementsOfDimToElementsOfDim(dim - 1, dim - 1) .push_back(subfacet_list); elementsOfDimToElementsOfDim(dim, dim - 1) .push_back(facet_list); elementsOfDimToElementsOfDim(dim + 1, dim - 1) .push_back(element_list); } } } else { std::vector element_list; std::vector facet_list; Vector subfacet_connectivity = mesh_facets.getConnectivity(subfacet); /// check if subfacet is to be doubled if (findElementsAroundSubfacet(starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false and removeElementsInVector(facet_list, elements_to_facets(subfacet)) == false) { Element new_subfacet{ subfacet.type, nb_new_facets(subfacet.type, subfacet.ghost_type)++, subfacet.ghost_type}; facets_to_double_by_dim[dim - 1]->push_back( Vector{subfacet, new_subfacet}); elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list); elementsOfDimToElementsOfDim(dim + 1, dim - 1) .push_back(element_list); } } } } } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::doubleNodes( const std::vector & old_nodes) { auto & position = mesh.getNodes(); auto spatial_dimension = mesh.getSpatialDimension(); auto old_nb_nodes = position.size(); position.reserve(old_nb_nodes + old_nodes.size()); doubled_nodes.reserve(doubled_nodes.size() + old_nodes.size()); auto position_begin = position.begin(spatial_dimension); for (auto && data : enumerate(old_nodes)) { auto n = std::get<0>(data); auto old_node = std::get<1>(data); decltype(old_node) new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes.push_back(Vector{old_node, new_node}); /// update position Vector coords = Vector(position_begin[old_node]); position.push_back(coords); } } /* -------------------------------------------------------------------------- */ bool CohesiveElementInserterHelper::findElementsAroundSubfacet( const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & element_list, std::vector & facet_list, std::vector * subfacet_list) { bool facet_matched = false; element_list.push_back(starting_element); std::queue elements_to_check; elements_to_check.push(starting_element); /// keep going as long as there are elements to check while (not elements_to_check.empty()) { /// check current element Element & current_element = elements_to_check.front(); const Vector facets_to_element = mesh_facets.getSubelementToElement(current_element); // for every facet of the element for (auto & current_facet : facets_to_element) { if (current_facet == ElementNull) { continue; } if (current_facet == end_facet) { facet_matched = true; } auto find_facet = std::find(facet_list.begin(), facet_list.end(), current_facet); // facet already listed or subfacet_connectivity is not in the // connectivity of current_facet; if ((find_facet != facet_list.end()) or not hasElement(mesh_facets.getConnectivity(current_facet), subfacet_connectivity)) { continue; } facet_list.push_back(current_facet); if (subfacet_list != nullptr) { const Vector subfacets_of_facet = mesh_facets.getSubelementToElement(current_facet); /// check subfacets for (const auto & current_subfacet : subfacets_of_facet) { if (current_subfacet == ElementNull) { continue; } auto find_subfacet = std::find( subfacet_list->begin(), subfacet_list->end(), current_subfacet); if ((find_subfacet != subfacet_list->end()) or not hasElement(mesh_facets.getConnectivity(current_subfacet), subfacet_connectivity)) { continue; } subfacet_list->push_back(current_subfacet); } } /// consider opposing element const auto & elements_to_facet = mesh_facets.getElementToSubelement(current_facet); UInt opposing = 0; if (elements_to_facet[0] == current_element) { opposing = 1; } auto & opposing_element = elements_to_facet[opposing]; /// skip null elements since they are on a boundary if (opposing_element == ElementNull) { continue; } /// skip this element if already added if (std::find(element_list.begin(), element_list.end(), opposing_element) != element_list.end()) { continue; } /// only regular elements have to be checked if (opposing_element.kind() == _ek_regular) { elements_to_check.push(opposing_element); } element_list.push_back(opposing_element); AKANTU_DEBUG_ASSERT(hasElement(mesh.getConnectivity(opposing_element), subfacet_connectivity), "Subfacet doesn't belong to this element"); } /// erased checked element from the list elements_to_check.pop(); } return facet_matched; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateCohesiveData() { UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; MeshAccessor mesh_accessor(mesh); MeshAccessor mesh_facets_accessor(mesh_facets); for (auto ghost_type : ghost_types) { for (auto cohesive_type : mesh.elementTypes(_element_kind = _ek_cohesive)) { auto nb_cohesive_elements = mesh.getNbElement(cohesive_type, ghost_type); nb_new_facets(cohesive_type, ghost_type) = nb_cohesive_elements; mesh_facets_accessor.getSubelementToElement(cohesive_type, ghost_type); } } auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; new_elements.reserve(new_elements.size() + facets_to_double.size()); std::array facets; auto & element_to_facet = mesh_facets_accessor.getElementToSubelement(); auto & facets_to_cohesive_elements = mesh_facets_accessor.getSubelementToElement(); for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double[0]; /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension ? 1 : 0] = old_facet; facets[third_dimension ? 0 : 1] = facet_to_double[1]; auto type_cohesive = FEEngine::getCohesiveElementType(old_facet.type); auto & facet_connectivities = mesh_facets.getConnectivity(old_facet.type, old_facet.ghost_type); auto facet_connectivity_it = facet_connectivities.begin(facet_connectivities.getNbComponent()); auto cohesive_element = Element{ type_cohesive, nb_new_facets(type_cohesive, old_facet.ghost_type)++, old_facet.ghost_type}; auto & cohesives_connectivities = mesh_accessor.getConnectivity(type_cohesive, old_facet.ghost_type); Matrix connectivity(facet_connectivities.getNbComponent(), 2); Vector facets_to_cohesive_element(2); for (auto s : arange(2)) { /// store doubled facets facets_to_cohesive_element[s] = facets[s]; // update connectivities connectivity(s) = Vector(facet_connectivity_it[facets[s].element]); /// update element_to_facet vectors element_to_facet(facets[s], 0)[1] = cohesive_element; } cohesives_connectivities.push_back(connectivity); facets_to_cohesive_elements(type_cohesive, old_facet.ghost_type) .push_back(facets_to_cohesive_element); /// add cohesive element to the element event list new_elements.push_back(cohesive_element); } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::doublePointFacet() { UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) { return; } NewElementsEvent new_facets_event; auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; const auto & element_to_facet = mesh_facets.getElementToSubelement(); auto & position = mesh.getNodes(); MeshAccessor mesh_accessor(mesh_facets); for (auto ghost_type : ghost_types) { for (auto facet_type : nb_new_facets.elementTypes( spatial_dimension - 1, ghost_type, _ek_regular)) { auto nb_new_element = nb_new_facets(facet_type, ghost_type); auto & connectivities = mesh_accessor.getConnectivity(facet_type, ghost_type); - connectivities.resize(connectivities.size() + nb_new_element); + connectivities.resize(nb_new_element); } } position.reserve(position.size() + facets_to_double.size()); for (auto facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double[0]; auto & new_facet = facet_to_double[1]; auto element = element_to_facet(new_facet)[0]; auto & facet_connectivities = mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type); auto old_node = facet_connectivities(old_facet.element); auto new_node = position.size(); /// update position position.push_back(position(old_node)); facet_connectivities(new_facet.element) = new_node; Vector segment_conectivity = mesh.getConnectivity(element); /// update facet connectivity auto it = std::find(segment_conectivity.begin(), segment_conectivity.end(), old_node); *it = new_node; doubled_nodes.push_back(Vector{old_node, new_node}); new_facets_event.getList().push_back(new_facet); } mesh_facets.sendEvent(new_facets_event); } /* -------------------------------------------------------------------------- */ template void CohesiveElementInserterHelper::doubleSubfacet() { if (spatial_dimension == 1) { return; } NewElementsEvent new_facets_event; std::vector nodes_to_double; MeshAccessor mesh_accessor(mesh_facets); auto & connectivities = mesh_accessor.getConnectivities(); auto & facets_to_double = *facets_to_double_by_dim[0]; ElementTypeMap nb_new_elements; for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto && old_element = facet_to_double[0]; nb_new_elements(old_element.type, old_element.ghost_type) = 0; } for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto && old_element = facet_to_double[0]; ++nb_new_elements(old_element.type, old_element.ghost_type); } for (auto ghost_type : ghost_types) { for (auto facet_type : nb_new_elements.elementTypes(0, ghost_type, _ek_regular)) { auto & connectivities = mesh_accessor.getConnectivity(facet_type, ghost_type); connectivities.resize(connectivities.size() + nb_new_elements(facet_type, ghost_type)); } } for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double(0); // auto & new_facet = facet_to_double(1); AKANTU_DEBUG_ASSERT( old_facet.type == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); /// list nodes double nodes_to_double.push_back(connectivities(old_facet)); } auto old_nb_doubled_nodes = doubled_nodes.size(); doubleNodes(nodes_to_double); auto double_nodes_view = make_view(doubled_nodes, 2); for (auto && data : zip(make_view(facets_to_double, 2), range(double_nodes_view.begin() + old_nb_doubled_nodes, double_nodes_view.end()), arange(facets_to_double.size()))) { // auto & old_facet = std::get<0>(data)[0]; auto & new_facet = std::get<0>(data)[1]; new_facets_event.getList().push_back(new_facet); auto & nodes = std::get<1>(data); auto old_node = nodes(0); auto new_node = nodes(1); auto f = std::get<2>(data); /// add new nodes in connectivity connectivities(new_facet) = new_node; updateElementalConnectivity(mesh, old_node, new_node, elementsOfDimToElementsOfDim(2, 0)(f), &elementsOfDimToElementsOfDim(1, 0)(f)); updateElementalConnectivity(mesh_facets, old_node, new_node, elementsOfDimToElementsOfDim(1, 0)(f)); if (spatial_dimension == 3) { updateElementalConnectivity(mesh_facets, old_node, new_node, elementsOfDimToElementsOfDim(0, 0)(f)); } } updateSubelementToElement(0, spatial_dimension == 2); updateElementToSubelement(0, spatial_dimension == 2); mesh_facets.sendEvent(new_facets_event); } } // namespace akantu