diff --git a/src/mesh/group_manager.cc b/src/mesh/group_manager.cc index bb1b89806..2955d59dd 100644 --- a/src/mesh/group_manager.cc +++ b/src/mesh/group_manager.cc @@ -1,1016 +1,1005 @@ /** * @file group_manager.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Feb 20 2018 * * @brief Stores information about ElementGroup and NodeGroup * * @section LICENSE * * Copyright (©) 2014-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 "group_manager.hh" #include "aka_csr.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ GroupManager::GroupManager(const Mesh & mesh, const ID & id, const MemoryID & mem_id) : id(id), memory_id(mem_id), mesh(mesh) { AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ GroupManager::~GroupManager() { auto eit = element_groups.begin(); auto eend = element_groups.end(); for (; eit != eend; ++eit) delete (eit->second); auto nit = node_groups.begin(); auto nend = node_groups.end(); for (; nit != nend; ++nit) delete (nit->second); } /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::createNodeGroup(const std::string & group_name, bool replace_group) { AKANTU_DEBUG_IN(); auto it = node_groups.find(group_name); if (it != node_groups.end()) { if (replace_group) { it->second->empty(); AKANTU_DEBUG_OUT(); return *(it->second); } else AKANTU_EXCEPTION( "Trying to create a node group that already exists:" << group_name); } std::stringstream sstr; sstr << this->id << ":" << group_name << "_node_group"; NodeGroup * node_group = new NodeGroup(group_name, mesh, sstr.str(), memory_id); node_groups[group_name] = node_group; AKANTU_DEBUG_OUT(); return *node_group; } /* -------------------------------------------------------------------------- */ template NodeGroup & GroupManager::createFilteredNodeGroup(const std::string & group_name, const NodeGroup & source_node_group, T & filter) { AKANTU_DEBUG_IN(); NodeGroup & node_group = this->createNodeGroup(group_name); node_group.append(source_node_group); if (T::type == FilterFunctor::_node_filter_functor) { node_group.applyNodeFilter(filter); } else { AKANTU_ERROR("ElementFilter cannot be applied to NodeGroup yet." << " Needs to be implemented."); } AKANTU_DEBUG_OUT(); return node_group; } /* -------------------------------------------------------------------------- */ void GroupManager::destroyNodeGroup(const std::string & group_name) { AKANTU_DEBUG_IN(); auto nit = node_groups.find(group_name); auto nend = node_groups.end(); if (nit != nend) { delete (nit->second); node_groups.erase(nit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, bool replace_group) { AKANTU_DEBUG_IN(); NodeGroup & new_node_group = createNodeGroup(group_name + "_nodes", replace_group); auto it = element_groups.find(group_name); if (it != element_groups.end()) { if (replace_group) { it->second->empty(); AKANTU_DEBUG_OUT(); return *(it->second); } else AKANTU_EXCEPTION("Trying to create a element group that already exists:" << group_name); } std::stringstream sstr; sstr << this->id << ":" << group_name << "_element_group"; ElementGroup * element_group = new ElementGroup( group_name, mesh, new_node_group, dimension, sstr.str(), memory_id); std::stringstream sstr_nodes; sstr_nodes << group_name << "_nodes"; node_groups[sstr_nodes.str()] = &new_node_group; element_groups[group_name] = element_group; AKANTU_DEBUG_OUT(); return *element_group; } /* -------------------------------------------------------------------------- */ void GroupManager::destroyElementGroup(const std::string & group_name, bool destroy_node_group) { AKANTU_DEBUG_IN(); auto eit = element_groups.find(group_name); auto eend = element_groups.end(); if (eit != eend) { if (destroy_node_group) destroyNodeGroup(eit->second->getNodeGroup().getName()); delete (eit->second); element_groups.erase(eit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::destroyAllElementGroups(bool destroy_node_groups) { AKANTU_DEBUG_IN(); auto eit = element_groups.begin(); auto eend = element_groups.end(); for (; eit != eend; ++eit) { if (destroy_node_groups) destroyNodeGroup(eit->second->getNodeGroup().getName()); delete (eit->second); } element_groups.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group) { AKANTU_DEBUG_IN(); if (element_groups.find(group_name) != element_groups.end()) AKANTU_EXCEPTION( "Trying to create a element group that already exists:" << group_name); ElementGroup * element_group = new ElementGroup(group_name, mesh, node_group, dimension, id + ":" + group_name + "_element_group", memory_id); element_groups[group_name] = element_group; AKANTU_DEBUG_OUT(); return *element_group; } /* -------------------------------------------------------------------------- */ template ElementGroup & GroupManager::createFilteredElementGroup( const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter) { AKANTU_DEBUG_IN(); ElementGroup * element_group = nullptr; if (T::type == FilterFunctor::_node_filter_functor) { NodeGroup & filtered_node_group = this->createFilteredNodeGroup( group_name + "_nodes", node_group, filter); element_group = &(this->createElementGroup(group_name, dimension, filtered_node_group)); } else if (T::type == FilterFunctor::_element_filter_functor) { AKANTU_ERROR( "Cannot handle an ElementFilter yet. Needs to be implemented."); } AKANTU_DEBUG_OUT(); return *element_group; } /* -------------------------------------------------------------------------- */ class ClusterSynchronizer : public DataAccessor { using DistantIDs = std::set>; public: ClusterSynchronizer(GroupManager & group_manager, UInt element_dimension, std::string cluster_name_prefix, ElementTypeMapArray & element_to_fragment, const ElementSynchronizer & element_synchronizer, UInt nb_cluster) : group_manager(group_manager), element_dimension(element_dimension), cluster_name_prefix(std::move(cluster_name_prefix)), element_to_fragment(element_to_fragment), element_synchronizer(element_synchronizer), nb_cluster(nb_cluster) {} UInt synchronize() { Communicator & comm = Communicator::getStaticCommunicator(); UInt rank = comm.whoAmI(); UInt nb_proc = comm.getNbProc(); /// find starting index to renumber local clusters Array nb_cluster_per_proc(nb_proc); nb_cluster_per_proc(rank) = nb_cluster; comm.allGather(nb_cluster_per_proc); starting_index = std::accumulate(nb_cluster_per_proc.begin(), nb_cluster_per_proc.begin() + rank, 0); UInt global_nb_fragment = std::accumulate(nb_cluster_per_proc.begin() + rank, nb_cluster_per_proc.end(), starting_index); /// create the local to distant cluster pairs with neighbors element_synchronizer.synchronizeOnce(*this, _gst_gm_clusters); /// count total number of pairs Array nb_pairs(nb_proc); // This is potentially a bug for more than // 2**31 pairs, but due to a all gatherv after // it must be int to match MPI interfaces nb_pairs(rank) = distant_ids.size(); comm.allGather(nb_pairs); UInt total_nb_pairs = std::accumulate(nb_pairs.begin(), nb_pairs.end(), 0); /// generate pairs global array UInt local_pair_index = std::accumulate(nb_pairs.storage(), nb_pairs.storage() + rank, 0); Array total_pairs(total_nb_pairs, 2); for (auto & ids : distant_ids) { total_pairs(local_pair_index, 0) = ids.first; total_pairs(local_pair_index, 1) = ids.second; ++local_pair_index; } /// communicate pairs to all processors nb_pairs *= 2; comm.allGatherV(total_pairs, nb_pairs); /// renumber clusters /// generate fragment list std::vector> global_clusters; UInt total_nb_cluster = 0; Array is_fragment_in_cluster(global_nb_fragment, 1, false); std::queue fragment_check_list; while (total_pairs.size() != 0) { /// create a new cluster ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; UInt first_fragment = total_pairs(0, 0); UInt second_fragment = total_pairs(0, 1); total_pairs.erase(0); fragment_check_list.push(first_fragment); fragment_check_list.push(second_fragment); while (!fragment_check_list.empty()) { UInt current_fragment = fragment_check_list.front(); UInt * total_pairs_end = total_pairs.storage() + total_pairs.size() * 2; UInt * fragment_found = std::find(total_pairs.storage(), total_pairs_end, current_fragment); if (fragment_found != total_pairs_end) { UInt position = fragment_found - total_pairs.storage(); UInt pair = position / 2; UInt other_index = (position + 1) % 2; fragment_check_list.push(total_pairs(pair, other_index)); total_pairs.erase(pair); } else { fragment_check_list.pop(); current_cluster.insert(current_fragment); is_fragment_in_cluster(current_fragment) = true; } } } /// add to FragmentToCluster all local fragments for (UInt c = 0; c < global_nb_fragment; ++c) { if (!is_fragment_in_cluster(c)) { ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; current_cluster.insert(c); } } /// reorganize element groups to match global clusters for (UInt c = 0; c < global_clusters.size(); ++c) { /// create new element group corresponding to current cluster std::stringstream sstr; sstr << cluster_name_prefix << "_" << c; ElementGroup & cluster = group_manager.createElementGroup(sstr.str(), element_dimension, true); auto it = global_clusters[c].begin(); auto end = global_clusters[c].end(); /// append to current element group all fragments that belong to /// the same cluster if they exist for (; it != end; ++it) { Int local_index = *it - starting_index; if (local_index < 0 || local_index >= Int(nb_cluster)) continue; std::stringstream tmp_sstr; tmp_sstr << "tmp_" << cluster_name_prefix << "_" << local_index; auto eg_it = group_manager.element_group_find(tmp_sstr.str()); AKANTU_DEBUG_ASSERT(eg_it != group_manager.element_group_end(), "Temporary fragment \"" << tmp_sstr.str() << "\" not found"); cluster.append(*(eg_it->second)); group_manager.destroyElementGroup(tmp_sstr.str(), true); } } return total_nb_cluster; } private: /// functions for parallel communications inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override { if (tag == _gst_gm_clusters) return elements.size() * sizeof(UInt); return 0; } inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag != _gst_gm_clusters) return; Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; /// for each element pack its global cluster index buffer << element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; } } inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag != _gst_gm_clusters) return; Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { UInt distant_cluster; buffer >> distant_cluster; const Element & el = *el_it; UInt local_cluster = element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; distant_ids.insert(std::make_pair(local_cluster, distant_cluster)); } } private: GroupManager & group_manager; UInt element_dimension; std::string cluster_name_prefix; ElementTypeMapArray & element_to_fragment; const ElementSynchronizer & element_synchronizer; UInt nb_cluster; DistantIDs distant_ids; UInt starting_index; }; /* -------------------------------------------------------------------------- */ /// \todo this function doesn't work in 1D UInt GroupManager::createBoundaryGroupFromGeometry() { UInt spatial_dimension = mesh.getSpatialDimension(); return createClusters(spatial_dimension - 1, "boundary"); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix, const GroupManager::ClusteringFilter & filter) { return createClusters(element_dimension, std::move(cluster_name_prefix), filter, mesh_facets); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, std::string cluster_name_prefix, const GroupManager::ClusteringFilter & filter) { std::unique_ptr mesh_facets; if (!mesh_facets && element_dimension > 0) { MeshAccessor mesh_accessor(const_cast(mesh)); mesh_facets = std::make_unique(mesh.getSpatialDimension(), mesh_accessor.getNodesSharedPtr(), "mesh_facets_for_clusters"); mesh_facets->defineMeshParent(mesh); MeshUtils::buildAllFacets(mesh, *mesh_facets, element_dimension, element_dimension - 1); } return createClusters(element_dimension, std::move(cluster_name_prefix), filter, *mesh_facets); } /* -------------------------------------------------------------------------- */ //// \todo if needed element list construction can be optimized by //// templating the filter class UInt GroupManager::createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const GroupManager::ClusteringFilter & filter, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt nb_proc = mesh.getCommunicator().getNbProc(); std::string tmp_cluster_name_prefix = cluster_name_prefix; ElementTypeMapArray * element_to_fragment = nullptr; if (nb_proc > 1 && mesh.isDistributed()) { element_to_fragment = new ElementTypeMapArray("element_to_fragment", id, memory_id); element_to_fragment->initialize( mesh, _nb_component = 1, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); // mesh.initElementTypeMapArray(*element_to_fragment, 1, element_dimension, // false, _ek_not_defined, true); tmp_cluster_name_prefix = "tmp_" + tmp_cluster_name_prefix; } ElementTypeMapArray seen_elements("seen_elements", id, memory_id); seen_elements.initialize(mesh, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); // mesh.initElementTypeMapArray(seen_elements, 1, element_dimension, false, // _ek_not_defined, true); for (auto ghost_type : ghost_types) { Element el; el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { el.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & seen_elements_array = seen_elements(type, ghost_type); for (UInt e = 0; e < nb_element; ++e) { el.element = e; if (!filter(el)) seen_elements_array(e) = true; } } } Array checked_node(mesh.getNbNodes(), 1, false); UInt nb_cluster = 0; for (auto ghost_type : ghost_types) { Element uns_el; uns_el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { uns_el.type = type; Array & seen_elements_vec = seen_elements(uns_el.type, uns_el.ghost_type); for (UInt e = 0; e < seen_elements_vec.size(); ++e) { // skip elements that have been already seen if (seen_elements_vec(e) == true) continue; // set current element uns_el.element = e; seen_elements_vec(e) = true; /// create a new cluster std::stringstream sstr; sstr << tmp_cluster_name_prefix << "_" << nb_cluster; ElementGroup & cluster = createElementGroup(sstr.str(), element_dimension, true); ++nb_cluster; // point element are cluster by themself if (element_dimension == 0) { cluster.add(uns_el); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(uns_el.type); Vector connect = mesh.getConnectivity(uns_el.type, uns_el.ghost_type) .begin(nb_nodes_per_element)[uns_el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { /// add element's nodes to the cluster UInt node = connect[n]; if (!checked_node(node)) { cluster.addNode(node); checked_node(node) = true; } } continue; } std::queue element_to_add; element_to_add.push(uns_el); /// keep looping until current cluster is complete (no more /// connected elements) while (!element_to_add.empty()) { /// take first element and erase it in the queue Element el = element_to_add.front(); element_to_add.pop(); /// if parallel, store cluster index per element if (nb_proc > 1 && mesh.isDistributed()) (*element_to_fragment)(el.type, el.ghost_type)(el.element) = nb_cluster - 1; /// add current element to the cluster cluster.add(el); const Array & element_to_facet = mesh_facets.getSubelementToElement(el.type, el.ghost_type); UInt nb_facet_per_element = element_to_facet.getNbComponent(); for (UInt f = 0; f < nb_facet_per_element; ++f) { const Element & facet = element_to_facet(el.element, f); if (facet == ElementNull) continue; const std::vector & connected_elements = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); for (UInt elem = 0; elem < connected_elements.size(); ++elem) { const Element & check_el = connected_elements[elem]; // check if this element has to be skipped if (check_el == ElementNull || check_el == el) continue; Array & seen_elements_vec_current = seen_elements(check_el.type, check_el.ghost_type); if (seen_elements_vec_current(check_el.element) == false) { seen_elements_vec_current(check_el.element) = true; element_to_add.push(check_el); } } } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = mesh.getConnectivity(el.type, el.ghost_type) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { /// add element's nodes to the cluster UInt node = connect[n]; if (!checked_node(node)) { cluster.addNode(node, false); checked_node(node) = true; } } } } } } if (nb_proc > 1 && mesh.isDistributed()) { ClusterSynchronizer cluster_synchronizer( *this, element_dimension, cluster_name_prefix, *element_to_fragment, this->mesh.getElementSynchronizer(), nb_cluster); nb_cluster = cluster_synchronizer.synchronize(); delete element_to_fragment; } if (mesh.isDistributed()) this->synchronizeGroupNames(); AKANTU_DEBUG_OUT(); return nb_cluster; } /* -------------------------------------------------------------------------- */ template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name) { std::set group_names; const auto & datas = mesh.getData(dataset_name); std::map group_dim; for (auto ghost_type : ghost_types) { for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements (" << type << ":" << ghost_type << ") in the map from MeshData (" << dataset.size() << ") " << dataset_name << " and in the mesh (" << nb_element << ")!"); for (UInt e(0); e < nb_element; ++e) { std::stringstream sstr; sstr << dataset(e); std::string gname = sstr.str(); group_names.insert(gname); auto it = group_dim.find(gname); if (it == group_dim.end()) { group_dim[gname] = mesh.getSpatialDimension(type); } else { it->second = std::max(it->second, mesh.getSpatialDimension(type)); } } } } auto git = group_names.begin(); auto gend = group_names.end(); for (; git != gend; ++git) createElementGroup(*git, group_dim[*git]); if (mesh.isDistributed()) this->synchronizeGroupNames(); Element el; for (auto ghost_type : ghost_types) { el.ghost_type = ghost_type; for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { el.type = type; const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements in the map from " "MeshData and in the mesh!"); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(el.type); Array::const_iterator> cit = mesh.getConnectivity(type, ghost_type).begin(nb_nodes_per_element); for (UInt e(0); e < nb_element; ++e, ++cit) { el.element = e; std::stringstream sstr; sstr << dataset(e); ElementGroup & group = getElementGroup(sstr.str()); group.add(el, false, false); const Vector & connect = *cit; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; group.addNode(node, false); } } } } git = group_names.begin(); for (; git != gend; ++git) { getElementGroup(*git).optimize(); } } template void GroupManager::createGroupsFromMeshData( const std::string & dataset_name); template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name); /* -------------------------------------------------------------------------- */ void GroupManager::createElementGroupFromNodeGroup( const std::string & name, const std::string & node_group_name, UInt dimension) { NodeGroup & node_group = getNodeGroup(node_group_name); ElementGroup & group = createElementGroup(name, dimension, node_group); group.fillFromNodeGroup(); } /* -------------------------------------------------------------------------- */ void GroupManager::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "GroupManager [" << std::endl; std::set node_group_seen; for (auto it(element_group_begin()); it != element_group_end(); ++it) { it->second->printself(stream, indent + 1); node_group_seen.insert(it->second->getNodeGroup().getName()); } for (auto it(node_group_begin()); it != node_group_end(); ++it) { if (node_group_seen.find(it->second->getName()) == node_group_seen.end()) it->second->printself(stream, indent + 1); } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ UInt GroupManager::getNbElementGroups(UInt dimension) const { if (dimension == _all_dimensions) return element_groups.size(); auto it = element_groups.begin(); auto end = element_groups.end(); UInt count = 0; for (; it != end; ++it) count += (it->second->getDimension() == dimension); return count; } /* -------------------------------------------------------------------------- */ -void GroupManager::checkAndAddGroups(CommunicationBuffer & buffer) { +void GroupManager::checkAndAddGroups(DynamicCommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); UInt nb_node_group; buffer >> nb_node_group; AKANTU_DEBUG_INFO("Received " << nb_node_group << " node group names"); for (UInt ng = 0; ng < nb_node_group; ++ng) { std::string node_group_name; buffer >> node_group_name; if (node_groups.find(node_group_name) == node_groups.end()) { this->createNodeGroup(node_group_name); } AKANTU_DEBUG_INFO("Received node goup name: " << node_group_name); } UInt nb_element_group; buffer >> nb_element_group; AKANTU_DEBUG_INFO("Received " << nb_element_group << " element group names"); for (UInt eg = 0; eg < nb_element_group; ++eg) { std::string element_group_name; buffer >> element_group_name; std::string node_group_name; buffer >> node_group_name; UInt dim; buffer >> dim; AKANTU_DEBUG_INFO("Received element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with node group " << node_group_name); NodeGroup & node_group = *node_groups[node_group_name]; if (element_groups.find(element_group_name) == element_groups.end()) { this->createElementGroup(element_group_name, dim, node_group); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::fillBufferWithGroupNames( DynamicCommunicationBuffer & comm_buffer) const { AKANTU_DEBUG_IN(); // packing node group names; UInt nb_groups = this->node_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " node group names"); auto nnames_it = node_groups.begin(); auto nnames_end = node_groups.end(); for (; nnames_it != nnames_end; ++nnames_it) { std::string node_group_name = nnames_it->first; comm_buffer << node_group_name; AKANTU_DEBUG_INFO("Sending node goupe name: " << node_group_name); } // packing element group names with there associated node group name nb_groups = this->element_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " element group names"); auto gnames_it = this->element_groups.begin(); auto gnames_end = this->element_groups.end(); for (; gnames_it != gnames_end; ++gnames_it) { ElementGroup & element_group = *(gnames_it->second); std::string element_group_name = gnames_it->first; std::string node_group_name = element_group.getNodeGroup().getName(); UInt dim = element_group.getDimension(); comm_buffer << element_group_name; comm_buffer << node_group_name; comm_buffer << dim; AKANTU_DEBUG_INFO("Sending element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with the node group " << node_group_name); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::synchronizeGroupNames() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); Int nb_proc = comm.getNbProc(); Int my_rank = comm.whoAmI(); if (nb_proc == 1) return; if (my_rank == 0) { for (Int p = 1; p < nb_proc; ++p) { - CommunicationStatus status; - comm.probe(p, p, status); - AKANTU_DEBUG_INFO("Got " << printMemorySize(status.size()) - << " from proc " << p); - - CommunicationBuffer recv_buffer(status.size()); + DynamicCommunicationBuffer recv_buffer; comm.receive(recv_buffer, p, p); - + AKANTU_DEBUG_INFO("Got " << printMemorySize(recv_buffer.size()) + << " from proc " << p); this->checkAndAddGroups(recv_buffer); } DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); - UInt buffer_size = comm_buffer.size(); - - comm.broadcast(buffer_size, 0); - AKANTU_DEBUG_INFO("Initiating broadcast with " << printMemorySize(comm_buffer.size())); - comm.broadcast(comm_buffer, 0); + comm.broadcast(comm_buffer); } else { DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); AKANTU_DEBUG_INFO("Sending " << printMemorySize(comm_buffer.size()) << " to proc " << 0); comm.send(comm_buffer, 0, my_rank); - UInt buffer_size = 0; - comm.broadcast(buffer_size, 0); - + DynamicCommunicationBuffer recv_buffer; + comm.broadcast(recv_buffer); AKANTU_DEBUG_INFO("Receiving broadcast with " - << printMemorySize(comm_buffer.size())); - CommunicationBuffer recv_buffer(buffer_size); - comm.broadcast(recv_buffer, 0); + << printMemorySize(recv_buffer.size())); this->checkAndAddGroups(recv_buffer); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const ElementGroup & GroupManager::getElementGroup(const std::string & name) const { auto it = element_group_find(name); if (it == element_group_end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::getElementGroup(const std::string & name) { auto it = element_group_find(name); if (it == element_group_end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const NodeGroup & GroupManager::getNodeGroup(const std::string & name) const { auto it = node_group_find(name); if (it == node_group_end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::getNodeGroup(const std::string & name) { auto it = node_group_find(name); if (it == node_group_end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } } // namespace akantu diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index f185d92e8..fd98c46b1 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,326 +1,326 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Feb 07 2018 * * @brief Stores information relevent to the notion of element and nodes * groups. * * @section LICENSE * * Copyright (©) 2014-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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GROUP_MANAGER_HH__ #define __AKANTU_GROUP_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_iterators.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template class CommunicationBufferTemplated; namespace dumper { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: #ifdef SWIGPYTHON public: using ElementGroups = std::map; using NodeGroups = std::map; private: #else using ElementGroups = std::map; using NodeGroups = std::map; #endif public: using GroupManagerTypeSet = std::set; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(const Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #ifndef SWIG #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ inline BOOST_PP_CAT(BOOST_PP_CAT(const_, group_type), _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ inline BOOST_PP_CAT(group_type, _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); #endif public: #ifndef SWIG decltype(auto) iterateNodeGroups() { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateNodeGroups() const { return make_dereference_adaptor(make_values_adaptor(node_groups)); } #endif /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* -------------------------------------------------------------------9+ ----- */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); /// destroy all element groups and the associated node groups void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) #ifndef SWIG template class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> dumper::Field * createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); template class ftype> dumper::Field * createNodalField(const ftype * field, const std::string & group_name, UInt padding_size = 0); template class ftype> dumper::Field * createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy - void checkAndAddGroups(CommunicationBufferTemplated & buffer); + void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template inline dumper::Field * createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template inline dumper::Field * createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); #endif /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// memory_id to create element and node groups MemoryID memory_id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs const Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh index 7548cd406..2f0f621dd 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,688 +1,688 @@ /** * @file mesh.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Feb 19 2018 * * @brief the class representing the 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_bbox.hh" #include "aka_event_handler_manager.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element.hh" #include "element_class.hh" #include "element_type_map.hh" #include "group_manager.hh" #include "mesh_data.hh" #include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { class Communicator; class ElementSynchronizer; class NodeSynchronizer; class PeriodicNodeSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Array, and the connectivity. The connectivity are stored in by element * types. * * In order to loop on all element you have to loop on all types like this : * @code{.cpp} for(auto & type : mesh.elementTypes()) { UInt nb_element = mesh.getNbElement(type); const Array & conn = mesh.getConnectivity(type); for(UInt e = 0; e < nb_element; ++e) { ... } } or for_each_element(mesh, [](Element & element) { std::cout << element << std::endl }); @endcode */ class Mesh : protected Memory, public EventHandlerManager, public GroupManager, public MeshData, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: /// default constructor used for chaining, the last parameter is just to /// differentiate constructors Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator); public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// mesh not distributed and not using the default communicator Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); ~Mesh() override; /// read the mesh from a file void read(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); /// write the mesh to a file void write(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); protected: void makeReady(); private: /// initialize the connectivity to NULL and other stuff void init(); /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /* ------------------------------------------------------------------------ */ /* Distributed memory methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// patitionate the mesh among the processors involved in their computation virtual void distribute(Communicator & communicator); virtual void distribute(); /// defines is the mesh is distributed or not inline bool isDistributed() const { return this->is_distributed; } /* ------------------------------------------------------------------------ */ /* Periodicity methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// set the periodicity in a given direction void makePeriodic(const SpatialDirection & direction); void makePeriodic(const SpatialDirection & direction, const ID & list_1, const ID & list_2); protected: void makePeriodic(const SpatialDirection & direction, const Array & list_1, const Array & list_2); /// Removes the face that the mesh is periodic void wipePeriodicInfo(); inline void addPeriodicSlave(UInt slave, UInt master); template void synchronizePeriodicSlaveDataWithMaster(Array & data); // update the periodic synchronizer (creates it if it does not exists) void updatePeriodicSynchronizer(); public: /// defines if the mesh is periodic or not inline bool isPeriodic() const { return (this->is_periodic != 0); } inline bool isPeriodic(const SpatialDirection & direction) const { return ((this->is_periodic & (1 << direction)) != 0); } class PeriodicSlaves; /// get the master node for a given slave nodes, except if node not a slave inline UInt getPeriodicMaster(UInt slave) const; #ifndef SWIG /// get an iterable list of slaves for a given master node inline decltype(auto) getPeriodicSlaves(UInt master) const; #endif /* ------------------------------------------------------------------------ */ /* General Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Array & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; // /// extract coordinates of nodes from a reversed element // inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, // UInt * connectivity, // UInt n_nodes); /// add a Array of connectivity for the type . inline void addConnectivityType(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { // if(event.getList().size() != 0) EventHandlerManager::sendEvent(event); } /// prepare the event to remove the elements listed void eraseElements(const Array & elements); /* ------------------------------------------------------------------------ */ template inline void removeNodesFromArray(Array & vect, const Array & new_numbering); /// initialize normals void initNormals(); /// init facets' mesh Mesh & initMeshFacets(const ID & id = "mesh_facets"); /// define parent mesh void defineMeshParent(const Mesh & mesh); /// get global connectivity array void getGlobalConnectivity(ElementTypeMapArray & global_connectivity); public: void getAssociatedElements(const Array & node_list, Array & elements); private: /// fills the nodes_to_elements structure void fillNodesToElements(); /// update the global ids, nodes type, ... std::tuple updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event); void registerGlobalDataUpdater( std::unique_ptr && global_data_updater); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the id of the mesh AKANTU_GET_MACRO(ID, Memory::id, const ID &); /// get the id of the mesh AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &); /// get the spatial dimension of the mesh = number of component of the /// coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Array aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Array &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array &); /// get the normals for the elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt); /// get the Array of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array &); - AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array &); + //AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global id of a node inline UInt getNodeLocalId(UInt global_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Array AKANTU_GET_MACRO(NodesFlags, *nodes_flags, const Array &); protected: AKANTU_GET_MACRO_NOT_CONST(NodesFlags, *nodes_flags, Array &); public: inline NodeFlag getNodeFlag(UInt local_id) const; inline Int getNodePrank(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; inline bool isPeriodicSlave(UInt n) const; inline bool isPeriodicMaster(UInt n) const; const Vector & getLowerBounds() const { return bbox.getLowerBounds(); } const Vector & getUpperBounds() const { return bbox.getUpperBounds(); } AKANTU_GET_MACRO(BBox, bbox, const BBox &); const Vector & getLocalLowerBounds() const { return bbox_local.getLowerBounds(); } const Vector & getLocalUpperBounds() const { return bbox_local.getUpperBounds(); } AKANTU_GET_MACRO(LocalBBox, bbox_local, const BBox &); /// get the connectivity Array for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ElementTypeMapArray &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the number of element for a given ghost_type and a given dimension inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_not_defined) const; /// compute the barycenter of a given element inline void getBarycenter(const Element & element, Vector & barycenter) const; void getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const; #ifndef SWIG /// get the element connected to a subelement (element of lower dimension) const auto & getElementToSubelement() const; /// get the element connected to a subelement const auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the element connected to a subelement auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the elements connected to a subelement const auto & getElementToSubelement(const Element & element) const; /// get the subelement (element of lower dimension) connected to a element const auto & getSubelementToElement() const; /// get the subelement connected to an element const auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the subelement connected to an element auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the subelement (element of lower dimension) connected to a element VectorProxy getSubelementToElement(const Element & element) const; /// get connectivity of a given element inline VectorProxy getConnectivity(const Element & element) const; inline Vector getConnectivityWithPeriodicity(const Element & element) const; protected: inline auto & getElementToSubelement(const Element & element); inline VectorProxy getSubelementToElement(const Element & element); inline VectorProxy getConnectivity(const Element & element); #endif public: /// get a name field associated to the mesh template inline const Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get a name field associated to the mesh template inline Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get a name field associated to the mesh template inline const ElementTypeMapArray & getData(const ID & data_name) const; /// get a name field associated to the mesh template inline ElementTypeMapArray & getData(const ID & data_name); template ElementTypeMap getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template dumper::Field * createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); /// templated getter returning the pointer to data in MeshData (modifiable) template inline Array & getDataPointer(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false); template inline Array & getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_); /// Facets mesh accessor inline const Mesh & getMeshFacets() const; inline Mesh & getMeshFacets(); /// Parent mesh accessor inline const Mesh & getMeshParent() const; inline bool isMeshFacets() const { return this->is_mesh_facets; } #ifndef SWIG /// return the dumper from a group and and a dumper name DumperIOHelper & getGroupDumper(const std::string & dumper_name, const std::string & group_name); #endif /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t); #ifndef SWIG /// get local connectivity of a facet for a given facet type static inline auto getFacetLocalConnectivity(const ElementType & type, UInt t = 0); /// get connectivity of facets for a given element inline auto getFacetConnectivity(const Element & element, UInt t = 0) const; #endif /// get the number of type of the surface element associated to a given /// element type static inline UInt getNbFacetTypes(const ElementType & type, UInt t = 0); #ifndef SWIG /// get the type of the surface element associated to a given element static inline constexpr auto getFacetType(const ElementType & type, UInt t = 0); /// get all the type of the surface element associated to a given element static inline constexpr auto getAllFacetTypes(const ElementType & type); #endif /// get the number of nodes in the given element list static inline UInt getNbNodesPerElementList(const Array & elements); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ #ifndef SWIG using type_iterator[[deprecated]] = ElementTypeMapArray::type_iterator; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; template ElementTypesIteratorHelper elementTypes(pack &&... _pack) const; [[deprecated("Use elementTypes instead")]] inline decltype(auto) firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).begin(); } [[deprecated("Use elementTypes instead")]] inline decltype(auto) lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).end(); } #endif AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer, const ElementSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer, ElementSynchronizer &); AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer, const NodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer, NodeSynchronizer &); AKANTU_GET_MACRO(PeriodicNodeSynchronizer, *periodic_node_synchronizer, const PeriodicNodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(PeriodicNodeSynchronizer, *periodic_node_synchronizer, PeriodicNodeSynchronizer &); // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator // &); #ifndef SWIG AKANTU_GET_MACRO(Communicator, *communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &); AKANTU_GET_MACRO(PeriodicMasterSlaves, periodic_master_slave, const auto &); #endif /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshAccessor; friend class MeshUtils; AKANTU_GET_MACRO(NodesPointer, *nodes, Array &); /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline Array & getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Array and create it if necessary inline Array & getNodesFlagsPointer(); /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline Array & getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the ghost element counter inline Array & getGhostsCounters(const ElementType & type, const GhostType & ghost_type = _ghost) { AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost, "No ghost counter for _not_ghost elements"); return ghosts_counters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline Array> & getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline Array & getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of the nodes coordinates std::shared_ptr> nodes; /// global node ids std::shared_ptr> nodes_global_ids; /// node flags (shared/periodic/...) std::shared_ptr> nodes_flags; /// processor handling the node when not local or master std::unordered_map nodes_prank; /// global number of nodes; UInt nb_global_nodes{0}; /// all class of elements present in this mesh (for heterogenous meshes) ElementTypeMapArray connectivities; /// count the references on ghost elements ElementTypeMapArray ghosts_counters; /// map to normals for all class of elements present in this mesh ElementTypeMapArray normals; /// the spatial dimension of this mesh UInt spatial_dimension{0}; /// size covered by the mesh on each direction Vector size; /// global bounding box BBox bbox; /// local bounding box BBox bbox_local; /// Extra data loaded from the mesh file // MeshData mesh_data; /// facets' mesh std::unique_ptr mesh_facets; /// parent mesh (this is set for mesh_facets meshes) const Mesh * mesh_parent{nullptr}; /// defines if current mesh is mesh_facets or not bool is_mesh_facets{false}; /// defines if the mesh is centralized or distributed bool is_distributed{false}; /// defines if the mesh is periodic bool is_periodic{false}; /// Communicator on which mesh is distributed Communicator * communicator; /// Element synchronizer std::unique_ptr element_synchronizer; /// Node synchronizer std::unique_ptr node_synchronizer; /// Node synchronizer for periodic nodes std::unique_ptr periodic_node_synchronizer; using NodesToElements = std::vector>>; /// class to update global data using external knowledge std::unique_ptr global_data_updater; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; /// periodicity local info std::unordered_map periodic_slave_master; std::unordered_multimap periodic_master_slave; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ #include "element_type_map_tmpl.hh" #include "mesh_inline_impl.cc" #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/mesh_utils/global_ids_updater.cc b/src/mesh_utils/global_ids_updater.cc index ee0c72503..ea2b7aac5 100644 --- a/src/mesh_utils/global_ids_updater.cc +++ b/src/mesh_utils/global_ids_updater.cc @@ -1,139 +1,133 @@ /** * @file global_ids_updater.cc * * @author Marco Vocialta * * @date creation: Fri Apr 13 2012 * @date last modification: Fri Dec 08 2017 * * @brief Functions of the GlobalIdsUpdater * * @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 "global_ids_updater.hh" #include "element_synchronizer.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { UInt GlobalIdsUpdater::updateGlobalIDs(UInt local_nb_new_nodes) { if (mesh.getCommunicator().getNbProc() == 1) return local_nb_new_nodes; UInt total_nb_new_nodes = this->updateGlobalIDsLocally(local_nb_new_nodes); if (mesh.isDistributed()) { this->synchronizeGlobalIDs(); } return total_nb_new_nodes; } UInt GlobalIdsUpdater::updateGlobalIDsLocally(UInt local_nb_new_nodes) { const auto & comm = mesh.getCommunicator(); - Int rank = comm.whoAmI(); Int nb_proc = comm.getNbProc(); if (nb_proc == 1) return local_nb_new_nodes; /// resize global ids array - Array & nodes_global_ids = mesh.getGlobalNodesIds(); + MeshAccessor mesh_accessor(mesh); + auto && nodes_global_ids = mesh_accessor.getNodesGlobalIds(); UInt old_nb_nodes = mesh.getNbNodes() - local_nb_new_nodes; nodes_global_ids.resize(mesh.getNbNodes(), -1); /// compute the number of global nodes based on the number of old nodes - Matrix local_master_nodes(2, nb_proc, 0); + Vector local_master_nodes(2, 0); for (UInt n = 0; n < old_nb_nodes; ++n) if (mesh.isLocalOrMasterNode(n)) - ++local_master_nodes(0, rank); + ++local_master_nodes(0); /// compute amount of local or master doubled nodes for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) if (mesh.isLocalOrMasterNode(n)) - ++local_master_nodes(1, rank); + ++local_master_nodes(1); + - comm.allGather(local_master_nodes); + auto starting_index = local_master_nodes(1); - local_master_nodes = local_master_nodes.transpose(); - UInt old_global_nodes = - std::accumulate(local_master_nodes(0).storage(), - local_master_nodes(0).storage() + nb_proc, 0); + comm.allReduce(local_master_nodes); - /// update global number of nodes - UInt total_nb_new_nodes = - std::accumulate(local_master_nodes(1).storage(), - local_master_nodes(1).storage() + nb_proc, 0); + UInt old_global_nodes = local_master_nodes(0); + UInt total_nb_new_nodes = local_master_nodes(1); if (total_nb_new_nodes == 0) return 0; /// set global ids of local and master nodes - UInt starting_index = - std::accumulate(local_master_nodes(1).storage(), - local_master_nodes(1).storage() + rank, old_global_nodes); + comm.exclusiveScan(starting_index); + starting_index += old_global_nodes; for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) { if (mesh.isLocalOrMasterNode(n)) { nodes_global_ids(n) = starting_index; ++starting_index; } else { nodes_global_ids(n) = -1; } } - MeshAccessor mesh_accessor(mesh); mesh_accessor.setNbGlobalNodes(old_global_nodes + total_nb_new_nodes); return total_nb_new_nodes; } void GlobalIdsUpdater::synchronizeGlobalIDs() { this->reduce = true; this->synchronizer.slaveReductionOnce(*this, _gst_giu_global_conn); #ifndef AKANTU_NDEBUG for (auto node : nodes_flags) { auto node_flag = mesh.getNodeFlag(node.first); if (node_flag != NodeFlag::_pure_ghost) continue; auto n = 0u; for (auto & pair : node.second) { if (std::get<1>(pair) == NodeFlag::_pure_ghost) ++n; } if (n == node.second.size()) { AKANTU_DEBUG_WARNING( "The node " << n << "is ghost on all the neighboring processors"); } } #endif this->reduce = false; this->synchronizer.synchronizeOnce(*this, _gst_giu_global_conn); } } // akantu diff --git a/src/mesh_utils/global_ids_updater_inline_impl.cc b/src/mesh_utils/global_ids_updater_inline_impl.cc index 4727529ca..fbd9a66b5 100644 --- a/src/mesh_utils/global_ids_updater_inline_impl.cc +++ b/src/mesh_utils/global_ids_updater_inline_impl.cc @@ -1,128 +1,129 @@ /** * @file global_ids_updater_inline_impl.cc * * @author Marco Vocialta * * @date creation: Fri Oct 02 2015 * @date last modification: Sun Aug 13 2017 * * @brief Implementation of the inline functions of GlobalIdsUpdater * * @section LICENSE * * Copyright (©) 2015-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 "communicator.hh" #include "global_ids_updater.hh" #include "mesh.hh" +#include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ #define __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt GlobalIdsUpdater::getNbData(const Array & elements, const SynchronizationTag & tag) const { UInt size = 0; - if (tag == _gst_giu_global_conn){ - size += Mesh::getNbNodesPerElementList(elements) * - sizeof(UInt); + if (tag == _gst_giu_global_conn) { + size += + Mesh::getNbNodesPerElementList(elements) * sizeof(UInt) + sizeof(int); #ifndef AKANTU_NDEBUG - size += sizeof(NodeFlag) * Mesh::getNbNodesPerElementList(elements) + sizeof(int); + size += sizeof(NodeFlag) * Mesh::getNbNodesPerElementList(elements); #endif } return size; } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag != _gst_giu_global_conn) return; auto & global_nodes_ids = mesh.getGlobalNodesIds(); -#ifndef AKANTU_NDEBUG buffer << int(mesh.getCommunicator().whoAmI()); -#endif + for (auto & element : elements) { /// get element connectivity const Vector current_conn = const_cast(mesh).getConnectivity(element); /// loop on all connectivity nodes for (auto node : current_conn) { UInt index = -1; if ((this->reduce and mesh.isLocalOrMasterNode(node)) or (not this->reduce and not mesh.isPureGhostNode(node))) { index = global_nodes_ids(node); } buffer << index; #ifndef AKANTU_NDEBUG buffer << mesh.getNodeFlag(node); #endif } } } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag != _gst_giu_global_conn) return; - auto & global_nodes_ids = mesh.getGlobalNodesIds(); + MeshAccessor mesh_accessor(mesh); + auto & global_nodes_ids = mesh_accessor.getNodesGlobalIds(); -#ifndef AKANTU_NDEBUG int proc; buffer >> proc; -#endif for (auto & element : elements) { /// get element connectivity Vector current_conn = const_cast(mesh).getConnectivity(element); /// loop on all connectivity nodes for (auto node : current_conn) { UInt index; buffer >> index; #ifndef AKANTU_NDEBUG NodeFlag node_flag; buffer >> node_flag; if (reduce) nodes_flags[node].push_back(std::make_pair(proc, node_flag)); #endif if (index == UInt(-1)) continue; - if (mesh.isSlaveNode(node)) + if (mesh.isSlaveNode(node)) { global_nodes_ids(node) = index; + mesh_accessor.setNodePrank(node, proc); + } } } } -} // akantu +} // namespace akantu #endif /* __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index 85cc447b1..b92532bbf 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,1821 +1,1810 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Marco Vocialta * * @date creation: Fri Aug 20 2010 * @date last modification: Wed Feb 21 2018 * * @brief All mesh utils necessary for various tasks * * @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 "mesh_utils.hh" #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" #include "mesh_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); for_each_element(mesh, [&](auto && element) { Vector conn = mesh.getConnectivity(element); for (auto && node : conn) { ++node_to_elem.rowOffset(node); } }, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list // Element e; node_to_elem.beginInsertions(); for_each_element(mesh, [&](auto && element) { Vector conn = mesh.getConnectivity(element); for (auto && node : conn) { node_to_elem.insertInRow(node, element); } }, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).size(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(spatial_dimension - 1, ghost_type)) { mesh.getConnectivity(type, ghost_type).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT( mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension); /// sort facets and generate sub-facets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension) { AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the // sub-elements // example: if the function is called with mesh = mesh_facets const Mesh * mesh_facets_parent = nullptr; try { mesh_facets_parent = &mesh_facets.getMeshParent(); } catch (...) { } mesh_facets.defineMeshParent(mesh); MeshAccessor mesh_accessor(mesh_facets); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const auto mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; std::vector connected_elements; // init the SubelementToElement data to improve performance for (auto && ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(dimension, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type); auto facet_types = mesh.getAllFacetTypes(type); for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); mesh_accessor.getElementToSubelement(facet_type, ghost_type); mesh_accessor.getConnectivity(facet_type, ghost_type); } } } const ElementSynchronizer * synchronizer = nullptr; if (mesh.isDistributed()) { synchronizer = &(mesh.getElementSynchronizer()); } Element current_element; for (auto && ghost_type : ghost_types) { GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; for (auto && type : mesh.elementTypes(dimension, ghost_type)) { auto facet_types = mesh.getAllFacetTypes(type); current_element.type = type; for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); auto nb_element = mesh.getNbElement(type, ghost_type); auto element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, ghost_type); auto connectivity_facets = &mesh_facets.getConnectivity(facet_type, ghost_type); auto nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft); const auto & element_connectivity = mesh.getConnectivity(type, ghost_type); Matrix facet_local_connectivity( mesh.getFacetLocalConnectivity(type, ft)); auto nb_nodes_per_facet = connectivity_facets->getNbComponent(); Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; for (UInt f = 0; f < nb_facet_per_element; ++f) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = element_connectivity(el, facet_local_connectivity(f, n)); UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.clear(); // loop over the other nodes to search intersecting elements, // which are the elements that share another node with the // starting element after first_node UInt local_el = 0; auto first_node_elements = node_to_elem.begin(facet(0)); auto first_node_elements_end = node_to_elem.end(facet(0)); for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { auto node_elements_begin = node_to_elem.begin(facet(n)); auto node_elements_end = node_to_elem.end(facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; connected_elements.clear(); for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (not(counter(el_f) == nb_nodes_per_facet - 1)) continue; ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } if (minimum_el != current_element) continue; bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; if (full_ghost_facet) continue; if (boundary_only and nb_element_connected_to_facet != 1) continue; std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); if (nb_element_connected_to_facet == 1) { /// boundary facet elements.push_back(ElementNull); } else if (nb_element_connected_to_facet == 2) { /// internal facet elements.push_back(connected_elements[1]); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = {_not_ghost, _not_ghost}; for (UInt el = 0; el < connected_elements.size(); ++el) gt[el] = connected_elements[el].ghost_type; if ((gt[0] == _not_ghost) xor (gt[1] == _not_ghost)) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { prank[el] = synchronizer->getRank(connected_elements[el]); } // ugly trick from Marco detected :P bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement( facet_type, facet_ghost_type); } } else { /// facet of facet for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { elements.push_back(connected_elements[i]); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->size() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type == _not_defined) continue; Array & subelement_to_element = mesh_facets.getSubelementToElement(loc_el.type, loc_el.ghost_type); UInt nb_facet_per_loc_element = subelement_to_element.getNbComponent(); for (UInt f_in = 0; f_in < nb_facet_per_loc_element; ++f_in) { auto & el = subelement_to_element(loc_el.element, f_in); if (el.type != _not_defined) continue; el.type = facet_type; el.element = current_facet; el.ghost_type = facet_ghost_type; break; } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = &mesh_accessor.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_accessor.getElementToSubelement( facet_type, facet_ghost_type); } } } } } } // restore the parent of mesh_facet if (mesh_facets_parent) mesh_facets.defineMeshParent(*mesh_facets_parent); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, Array & local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.size(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element) * nb_nodes_per_element, renumbering_map); old_nodes_numbers.resize(renumbering_map.size()); for (auto & renumber_pair : renumbering_map) { old_nodes_numbers(renumber_pair.second) = renumber_pair.first; } renumbering_map.clear(); MeshAccessor mesh_accessor(mesh); /// copy the renumbered connectivity to the right place auto & local_conn = mesh_accessor.getConnectivity(type); local_conn.resize(nb_local_element); memcpy(local_conn.storage(), local_connectivities.storage(), nb_local_element * nb_nodes_per_element * sizeof(UInt)); auto & ghost_conn = mesh_accessor.getConnectivity(type, _ghost); ghost_conn.resize(nb_ghost_element); std::memcpy(ghost_conn.storage(), local_connectivities.storage() + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost); ghost_counter.resize(nb_ghost_element, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity( Array & list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes.storage(); UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; auto it = renumbering_map.find(node); if (it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(_all_dimensions, ghost_type, _ek_not_defined)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); Array & connectivity = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity.size()); renumberNodesInConnectivity( connectivity, nb_element * nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); for (auto && pair : renumbering_map) { new_numbering(std::get<0>(pair)) = std::get<1>(pair); } for (UInt i = 0; i < new_numbering.size(); ++i) { if (new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_COHESIVE_ELEMENT) /* -------------------------------------------------------------------------- */ UInt MeshUtils::insertCohesiveElements( Mesh & mesh, Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion, Array & doubled_nodes, Array & new_elements, bool only_double_facets) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion); if (elements_to_insert > 0) { if (spatial_dimension == 1) { doublePointFacet(mesh, mesh_facets, doubled_nodes); } else { doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes, true); findSubfacetToDouble(mesh_facets); if (spatial_dimension == 2) { doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes); } else if (spatial_dimension == 3) { doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false); findSubfacetToDouble(mesh_facets); doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes); } } if (!only_double_facets) updateCohesiveData(mesh, mesh_facets, new_elements); } return elements_to_insert; } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doubleNodes(Mesh & mesh, const std::vector & old_nodes, Array & doubled_nodes) { AKANTU_DEBUG_IN(); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt old_nb_nodes = position.size(); UInt new_nb_nodes = old_nb_nodes + old_nodes.size(); UInt old_nb_doubled_nodes = doubled_nodes.size(); UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size(); position.resize(new_nb_nodes); doubled_nodes.resize(new_nb_doubled_nodes); Array::iterator> position_begin = position.begin(spatial_dimension); for (UInt n = 0; n < old_nodes.size(); ++n) { UInt new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n]; doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node; /// update position std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1, position_begin + new_node); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, UInt facet_dimension, Array & doubled_nodes, bool facet_mode) { AKANTU_DEBUG_IN(); for (auto gt_facet : ghost_types) { for (auto && type_facet : mesh_facets.elementTypes(facet_dimension, gt_facet)) { auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); if (nb_facet_to_double == 0) continue; // this while fail if multiple facet types // \TODO handle multiple sub-facet types auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto nb_nodes_per_facet = conn_facet.getNbComponent(); auto old_nb_facet = conn_facet.size(); auto new_nb_facet = old_nb_facet + nb_facet_to_double; #ifndef AKANTU_NDEBUG // memory initialization are slow but help debug conn_facet.resize(new_nb_facet, UInt(-1)); #else conn_facet.resize(new_nb_facet); #endif auto conn_facet_begin = conn_facet.begin(nb_nodes_per_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); #ifndef AKANTU_NDEBUG subfacet_to_facet.resize(new_nb_facet, ElementNull); #else subfacet_to_facet.resize(new_nb_facet); #endif auto subfacet_to_facet_begin = subfacet_to_facet.begin(nb_subfacet_per_facet); Element new_facet{type_facet, old_nb_facet, gt_facet}; auto conn_facet_new_it = conn_facet_begin + new_facet.element; auto subfacet_to_facet_new_it = subfacet_to_facet_begin + new_facet.element; for (UInt facet = 0; facet < nb_facet_to_double; ++facet, ++new_facet.element, ++conn_facet_new_it, ++subfacet_to_facet_new_it) { UInt old_facet = facets_to_double(facet); /// adding a new facet by copying original one /// copy connectivity in new facet *conn_facet_new_it = conn_facet_begin[old_facet]; /// update subfacet_to_facet *subfacet_to_facet_new_it = subfacet_to_facet_begin[old_facet]; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; /// update facet_to_subfacet array mesh_facets.getElementToSubelement(subfacet).push_back(new_facet); } } /// update facet_to_subfacet and _segment_3 facets if any if (not facet_mode) { updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true); updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true); updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } else updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt MeshUtils::updateFacetToDouble( Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); UInt nb_facets_to_double = 0.; 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 & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); - ElementType el_type = _not_defined; - GhostType el_gt = _casper; - UInt nb_facet_per_element = 0; Element old_facet_el{type_facet, 0, gt_facet}; - - Array * facet_to_element = nullptr; + UInt nb_facets = mesh_facets.getNbElement(type_facet, gt_facet); for (UInt f = 0; f < f_insertion.size(); ++f) { if (f_insertion(f) == false) continue; ++nb_facets_to_double; if (element_to_facet(f)[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || element_to_facet(f)[1].kind() == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } f_to_double.push_back(f); - UInt new_facet = mesh_facets.getNbElement(type_facet, gt_facet) + - f_to_double.size() - 1; + UInt new_facet = nb_facets + f_to_double.size() - 1; old_facet_el.element = f; /// update facet_to_element vector - Element & elem_to_update = element_to_facet(f)[1]; + auto & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.element; - if (elem_to_update.ghost_type != el_gt || - elem_to_update.type != el_type) { - el_type = elem_to_update.type; - el_gt = elem_to_update.ghost_type; - facet_to_element = - &mesh_facets.getSubelementToElement(el_type, el_gt); - nb_facet_per_element = facet_to_element->getNbComponent(); - } - - auto begin = facet_to_element->storage() + el * nb_facet_per_element; + auto & facet_to_element = mesh_facets.getSubelementToElement( + elem_to_update.type, elem_to_update.ghost_type); + auto el_facets = Vector( + make_view(facet_to_element, facet_to_element.getNbComponent()) + .begin()[el]); auto f_update = - std::find(begin, begin + nb_facet_per_element, old_facet_el); + std::find(el_facets.begin(), el_facets.end(), old_facet_el); - AKANTU_DEBUG_ASSERT(f_update != begin + nb_facet_per_element, - "Facet not found"); + AKANTU_DEBUG_ASSERT(f_update != el_facets.end(), "Facet not found"); f_update->element = new_facet; /// update elements connected to facet - std::vector first_facet_list = element_to_facet(f); + const auto & first_facet_list = element_to_facet(f); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1]; + element_to_facet(new_facet)[1] = ElementNull; element_to_facet(f)[1] = ElementNull; - element_to_facet(new_facet)[1] = ElementNull; } } } AKANTU_DEBUG_OUT(); return nb_facets_to_double; } /* -------------------------------------------------------------------------- */ void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { for (auto type : mesh_facets.elementTypes(_all_dimensions, gt)) { mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false); mesh_facets.getDataPointer>( "facets_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "elements_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "subfacets_to_subsubfacet_double", type, gt, 1, false); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::findSubfacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 1) return; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); GhostType gt_subfacet = _casper; ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet); GhostType gt_subsubfacet = _casper; Array * conn_subfacet = nullptr; Array * sf_to_double = nullptr; Array> * sf_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; Array> * el_to_subfacet_double = nullptr; UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet); UInt nb_subsubfacet; UInt nb_nodes_per_sf_el; if (subsubfacet_mode) { nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet); nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet); } else nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array * subsubfacet_to_subfacet = nullptr; UInt old_nb_facet = subfacet_to_facet.size() - nb_facet_to_double; Element current_facet{type_facet, 0, gt_facet}; std::vector element_list; std::vector facet_list; std::vector * subfacet_list; if (subsubfacet_mode) subfacet_list = new std::vector; /// map to filter subfacets Array> * facet_to_subfacet = nullptr; /// this is used to make sure that both new and old facets are /// checked UInt facets[2]; /// loop on every facet for (UInt f_index = 0; f_index < 2; ++f_index) { for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { facets[bool(f_index)] = facets_to_double(facet); facets[!bool(f_index)] = old_nb_facet + facet; UInt old_facet = facets[0]; UInt new_facet = facets[1]; Element & starting_element = element_to_facet(new_facet)[0]; current_facet.element = old_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; if (subsubfacet_mode) { subsubfacet_to_subfacet = &mesh_facets.getSubelementToElement( type_subfacet, gt_subfacet); } else { conn_subfacet = &mesh_facets.getConnectivity(type_subfacet, gt_subfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subfacet, gt_subfacet); } } if (subsubfacet_mode) { /// loop on every subsubfacet for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) { Element & subsubfacet = (*subsubfacet_to_subfacet)(subfacet.element, ssf); if (subsubfacet == ElementNull) continue; if (gt_subsubfacet != subsubfacet.ghost_type) { gt_subsubfacet = subsubfacet.ghost_type; conn_subfacet = &mesh_facets.getConnectivity(type_subsubfacet, gt_subsubfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subsubfacet, gt_subsubfacet); sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subsubfacet, gt_subsubfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subsubfacet, gt_subsubfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subsubfacet, gt_subsubfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subsubfacet, gt_subsubfacet); } UInt global_ssf = subsubfacet.element; Vector subsubfacet_connectivity( conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet( mesh_facets, starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, subfacet_list) == false && removeElementsInVector(*subfacet_list, (*facet_to_subfacet)(global_ssf)) == false) { sf_to_double->push_back(global_ssf); sf_to_subfacet_double->push_back(*subfacet_list); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } else { const UInt global_sf = subfacet.element; Vector subfacet_connectivity( conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subfacet is to be doubled if (findElementsAroundSubfacet( mesh_facets, starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false && removeElementsInVector( facet_list, (*facet_to_subfacet)(global_sf)) == false) { sf_to_double->push_back(global_sf); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } } } if (subsubfacet_mode) delete subfacet_list; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; MeshAccessor mesh_facets_accessor(mesh_facets); for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); auto & facet_to_coh_element = mesh_facets_accessor.getSubelementToElement(type_cohesive, gt_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet); UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive.size(); UInt new_nb_cohesive_elements = conn_cohesive.size() + nb_facet_to_double; UInt old_nb_facet = element_to_facet.size() - nb_facet_to_double; facet_to_coh_element.resize(new_nb_cohesive_elements); conn_cohesive.resize(new_nb_cohesive_elements); UInt new_elements_old_size = new_elements.size(); new_elements.resize(new_elements_old_size + nb_facet_to_double); Element c_element{type_cohesive, 0, gt_facet}; Element f_element{type_facet, 0, gt_facet}; UInt facets[2]; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension ? 1 : 0] = f_to_double(facet); facets[third_dimension ? 0 : 1] = old_nb_facet + facet; UInt cohesive_element = old_nb_cohesive_elements + facet; /// store doubled facets f_element.element = facets[0]; facet_to_coh_element(cohesive_element, 0) = f_element; f_element.element = facets[1]; facet_to_coh_element(cohesive_element, 1) = f_element; /// modify cohesive elements' connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { conn_cohesive(cohesive_element, n) = conn_facet(facets[0], n); conn_cohesive(cohesive_element, n + nb_nodes_per_facet) = conn_facet(facets[1], n); } /// update element_to_facet vectors c_element.element = cohesive_element; element_to_facet(facets[0])[1] = c_element; element_to_facet(facets[1])[1] = c_element; /// add cohesive element to the element event list new_elements(new_elements_old_size + facet) = c_element; } } } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) return; auto & position = mesh.getNodes(); for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); const auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); auto new_nb_facet = element_to_facet.size(); auto old_nb_facet = element_to_facet.size() - nb_facet_to_double; auto old_nb_nodes = position.size(); auto new_nb_nodes = old_nb_nodes + nb_facet_to_double; position.resize(new_nb_nodes); conn_facet.resize(new_nb_facet); auto old_nb_doubled_nodes = doubled_nodes.size(); doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double); for (auto && data_facet : enumerate(facets_to_double)) { const auto & old_facet = std::get<1>(data_facet); auto facet = std::get<0>(data_facet); auto new_facet = old_nb_facet + facet; auto el = element_to_facet(new_facet)[0]; auto old_node = conn_facet(old_facet); auto new_node = old_nb_nodes + facet; /// update position position(new_node) = position(old_node); conn_facet(new_facet) = new_node; Vector conn_segment = mesh.getConnectivity(el); /// update facet connectivity auto it = std::find(conn_segment.begin(), conn_segment.end(), old_node); *it = new_node; doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node; doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, GhostType gt_facet, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (type_facet != _segment_3) return; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); UInt old_nb_facet = mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); /// this ones matter only for segments in 3D Array> * el_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; if (third_dim_segments) { el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_facet, gt_facet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_facet, gt_facet); } std::vector middle_nodes; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt node = conn_facet(old_facet, 2); if (!mesh.isPureGhostNode(node)) middle_nodes.push_back(node); } UInt n = doubled_nodes.size(); doubleNodes(mesh, middle_nodes, doubled_nodes); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt old_node = conn_facet(old_facet, 2); if (mesh.isPureGhostNode(old_node)) continue; UInt new_node = doubled_nodes(n, 1); UInt new_facet = old_nb_facet + facet; conn_facet(new_facet, 2) = new_node; if (third_dim_segments) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_facet(new_facet)); updateElementalConnectivity(mesh, old_node, new_node, (*el_to_subfacet_double)(facet), &(*f_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_facet(new_facet)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); /// update subfacet_to_facet vector ElementType type_facet = _not_defined; GhostType gt_facet = _casper; Array * subfacet_to_facet = nullptr; UInt nb_subfacet_per_facet = 0; UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) - nb_subfacet_to_double; Array> * facet_list = nullptr; if (facet_mode) facet_list = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); else facet_list = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); Element old_subfacet_el{type_subfacet, 0, gt_subfacet}; Element new_subfacet_el{type_subfacet, 0, gt_subfacet}; for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { old_subfacet_el.element = sf_to_double(sf); new_subfacet_el.element = old_nb_subfacet + sf; for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) { Element & facet = (*facet_list)(sf)[f]; if (facet.type != type_facet || facet.ghost_type != gt_facet) { type_facet = facet.type; gt_facet = facet.ghost_type; subfacet_to_facet = &mesh_facets.getSubelementToElement(type_facet, gt_facet); nb_subfacet_per_facet = subfacet_to_facet->getNbComponent(); } Element * sf_update = std::find( subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet, subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, old_subfacet_el); AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet != subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, "Subfacet not found"); *sf_update = new_subfacet_el; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); Array> & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); Array> * facet_to_subfacet_double = nullptr; if (facet_mode) { facet_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); } else { facet_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); } UInt old_nb_subfacet = facet_to_subfacet.size(); facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (spatial_dimension == 1) return; for (auto gt_subfacet : ghost_types) { for (auto type_subfacet : mesh_facets.elementTypes(0, gt_subfacet)) { auto & sf_to_double = mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); if (nb_subfacet_to_double == 0) continue; AKANTU_DEBUG_ASSERT( type_subfacet == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); auto & conn_subfacet = mesh_facets.getConnectivity(type_subfacet, gt_subfacet); UInt old_nb_subfacet = conn_subfacet.size(); UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double; conn_subfacet.resize(new_nb_subfacet); std::vector nodes_to_double; UInt old_nb_doubled_nodes = doubled_nodes.size(); /// double nodes for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_subfacet = sf_to_double(sf); nodes_to_double.push_back(conn_subfacet(old_subfacet)); } doubleNodes(mesh, nodes_to_double, doubled_nodes); /// add new nodes in connectivity for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt new_subfacet = old_nb_subfacet + sf; UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); conn_subfacet(new_subfacet) = new_node; } /// update facet and element connectivity Array> & f_to_subfacet_double = mesh_facets.getData>("facets_to_subfacet_double", type_subfacet, gt_subfacet); Array> & el_to_subfacet_double = mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); Array> * sf_to_subfacet_double = nullptr; if (spatial_dimension == 3) sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0); UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); updateElementalConnectivity(mesh, old_node, new_node, el_to_subfacet_double(sf), &f_to_subfacet_double(sf)); updateElementalConnectivity(mesh_facets, old_node, new_node, f_to_subfacet_double(sf)); if (spatial_dimension == 3) updateElementalConnectivity(mesh_facets, old_node, new_node, (*sf_to_subfacet_double)(sf)); } if (spatial_dimension == 2) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true); } else if (spatial_dimension == 3) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets( Mesh & mesh_facets, const ElementTypeMapArray & remote_global_connectivities, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray local_global_connectivities( "local_global_connectivity", mesh_facets.getID(), mesh_facets.getMemoryID()); local_global_connectivities.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _ghost_type = gt_facet, _with_nb_nodes_per_element = true, _with_nb_element = true); mesh_facets.getGlobalConnectivity(local_global_connectivities); /// loop on every facet for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet); auto & local_global_connectivity = local_global_connectivities(type_facet, gt_facet); const auto & remote_global_connectivity = remote_global_connectivities(type_facet, gt_facet); auto & element_per_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); auto nb_nodes_per_facet = connectivity.getNbComponent(); auto nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); for (auto && data : zip(make_view(connectivity, nb_nodes_per_facet), make_view(local_global_connectivity, nb_nodes_per_facet), make_view(remote_global_connectivity, nb_nodes_per_facet), make_view(subfacet_to_facet, subfacet_to_facet.getNbComponent()), make_view(element_per_facet))) { auto & conn = std::get<0>(data); auto & local_gconn = std::get<1>(data); const auto & remote_gconn = std::get<2>(data); /// skip facet if connectivities are the same if (local_gconn == remote_gconn) continue; /// re-arrange connectivity auto conn_tmp = conn; auto begin = local_gconn.begin(); auto end = local_gconn.end(); std::transform(remote_gconn.begin(), remote_gconn.end(), conn.begin(), [&](auto && gnode) { auto it = std::find(begin, end, gnode); AKANTU_DEBUG_ASSERT(it != end, "Node not found"); return conn_tmp(it - begin); }); /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { auto begin = remote_gconn.storage(); /// find first node auto it = std::find(begin, begin + remote_gconn.size(), local_gconn(0)); UInt n, start = it - begin; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (n = 1; n < nb_nodes_per_P1_facet && local_gconn(n) == remote_gconn((start + n) % nb_nodes_per_P1_facet); ++n) ; /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) continue; } /// update data to invert facet auto & element_per_facet = std::get<4>(data); std::swap(element_per_facet[0], element_per_facet[1]); auto & subfacets_of_facet = std::get<3>(data); std::swap(subfacets_of_facet(0), subfacets_of_facet(1)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call " "the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID(), mesh.getMemoryID()); barycenters.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = _all_dimensions); // mesh.initElementTypeMapArray(barycenters, spatial_dimension, // _all_dimensions); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(_all_dimensions, ghost_type)) { element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & barycenters_arr = barycenters(type, ghost_type); barycenters_arr.resize(nb_element); auto bary = barycenters_arr.begin(spatial_dimension); auto bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { element.element = el; mesh.getBarycenter(element, *bary); } } } MeshAccessor mesh_accessor(mesh); for (Int sp(spatial_dimension); sp >= 1; --sp) { if (mesh.getNbElement(sp) == 0) continue; for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(sp, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh_accessor.getSubelementToElement(type, ghost_type).set(ElementNull); } for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { mesh_accessor.getElementToSubelement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh.getElementToSubelement(type, ghost_type).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (auto ghost_type : ghost_types) { facet_element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { facet_element.type = type; auto & element_to_subelement = mesh.getElementToSubelement(type, ghost_type); const auto & connectivity = mesh.getConnectivity(type, ghost_type); for (auto && data : enumerate( make_view(connectivity, mesh.getNbNodesPerElement(type)))) { const auto & facet = std::get<1>(data); facet_element.element = std::get<0>(data); std::map element_seen_counter; auto nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(type)); // count the number of node in common between the facet and the other // element connected to the nodes of the facet for (auto node : arange(nb_nodes_per_facet)) { for (auto & elem : nodes_to_elements.getRow(facet(node))) { auto cit = element_seen_counter.find(elem); if (cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } // check which are the connected elements std::vector connected_elements; for (auto && cit : element_seen_counter) { if (cit.second == nb_nodes_per_facet) connected_elements.push_back(cit.first); } // add the connected elements as sub-elements for (auto & connected_element : connected_elements) { element_to_subelement(facet_element.element) .push_back(connected_element); } // add the element as sub-element to the connected elements for (auto & connected_element : connected_elements) { Vector subelements_to_element = mesh.getSubelementToElement(connected_element); // find the position where to insert the element auto it = std::find(subelements_to_element.begin(), subelements_to_element.end(), ElementNull); AKANTU_DEBUG_ASSERT( it != subelements_to_element.end(), "The element " << connected_element << " seems to have too many facets!! (" << (it - subelements_to_element.begin()) << " < " << mesh.getNbFacetsPerElement(connected_element.type) << ")"); *it = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MeshUtils::findElementsAroundSubfacet( const Mesh & mesh_facets, const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & element_list, std::vector & facet_list, std::vector * subfacet_list) { AKANTU_DEBUG_IN(); bool facet_matched = false; element_list.clear(); facet_list.clear(); if (third_dim_points) { subfacet_list->clear(); } 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; // facet already listed if (std::find(facet_list.begin(), facet_list.end(), current_facet) != facet_list.end()) continue; // subfacet_connectivity is not in the connectivity of current_facet; if ((std::find(facet_list.begin(), facet_list.end(), current_facet) != facet_list.end()) or not hasElement(mesh_facets.getConnectivity(current_facet), subfacet_connectivity)) continue; facet_list.push_back(current_facet); if (third_dim_points) { 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; if ((std::find(subfacet_list->begin(), subfacet_list->end(), current_subfacet) == subfacet_list->end()) and hasElement(mesh_facets.getConnectivity(current_subfacet), subfacet_connectivity)) 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_facets.getMeshParent().getConnectivity(opposing_element), subfacet_connectivity), "Subfacet doesn't belong to this element"); } /// erased checked element from the list elements_to_check.pop(); } AKANTU_DEBUG_OUT(); return facet_matched; } /* -------------------------------------------------------------------------- */ void MeshUtils::updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, const std::vector * #if defined(AKANTU_COHESIVE_ELEMENT) facet_list #endif ) { AKANTU_DEBUG_IN(); for (auto & element : element_list) { if (element.type == _not_defined) continue; Vector connectivity = mesh.getConnectivity(element); #if defined(AKANTU_COHESIVE_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.getMeshFacets().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() + 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 #endif { 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(); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index c434ffa6d..3fd2238de 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,707 +1,705 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Fabian Barras * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Solid mechanics model for cohesive elements * * @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 "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), *model.cohesive_synchronizer) {} /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto cohesive_nodes_event = dynamic_cast(&nodes_event); if (not cohesive_nodes_event) return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); auto & nodes_flags = mesh_accessor.getNodesFlags(); auto nb_old_nodes = nodes_flags.size(); nodes_flags.resize(nb_old_nodes + local_nb_new_nodes); for (auto && data : zip(old_nodes, new_nodes)) { UInt old_node, new_node; std::tie(old_node, new_node) = data; nodes_flags(new_node) = nodes_flags(old_node); } model.updateCohesiveSynchronizers(); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } Vector nb_new_stuff = {nb_new_nodes, elements_event.getList().size()}; const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_stuff, SynchronizerOperation::_sum); if (nb_new_stuff(1) > 0) { mesh.sendEvent(elements_event); MeshUtils::resetFacetToDouble(mesh.getMeshFacets()); } - if (nb_new_nodes > 0) { + if (nb_new_stuff(0) > 0) { mesh.sendEvent(nodes_event); // mesh.sendEvent(global_ids_updater.getChangedNodeEvent()); } - return std::make_tuple(nb_new_nodes, nb_new_stuff(1)); + return std::make_tuple(nb_new_stuff(0), nb_new_stuff(1)); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, _gst_material_id); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.is_extrinsic; - std::cout << "Extrinsic " << is_extrinsic << std::endl; - inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = dynamic_cast(mesh_facets.getElementSynchronizer()); synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, _gst_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt cohesive_index = UInt(-1); for (auto && material : enumerate(materials)) { if (dynamic_cast(std::get<1>(material).get())) { cohesive_index = std::get<0>(material); break; } } if (cohesive_index == UInt(-1)) AKANTU_EXCEPTION("No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = dynamic_cast(*materials[mat_index]); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initMaterials(); if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.size() == 0) continue; type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) continue; Matrix quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; // for (UInt q = 0; q < nb_quad_per_facet; ++q) { // for (UInt s = 0; s < Model::spatial_dimension; ++s) { // el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + // f * nb_quad_per_facet + q, // s) = quad_f(global_facet * nb_quad_per_facet + q, // s); // } // } //} } } } /// loop over non cohesive materials for (auto && material : materials) { if (dynamic_cast(material.get())) continue; /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = dynamic_cast(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { auto * mat = dynamic_cast(material.get()); if (mat == nullptr) /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } this->synchronize(_gst_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); for (auto & mat : materials) { auto * mat_cohesive = dynamic_cast(mat.get()); if (mat_cohesive) { /// check which not ghost cohesive elements are to be created mat_cohesive->checkInsertion(); } } /// communicate data among processors // this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onElementsAdded(element_list, event); if (is_extrinsic) resizeFacetStress(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onNodesAdded(new_nodes, event); const CohesiveNewNodesEvent * cohesive_event; if ((cohesive_event = dynamic_cast(&event)) == nullptr) return; const auto & old_nodes = cohesive_event->getOldNodesList(); auto copy = [this, &new_nodes, &old_nodes](auto & arr) { UInt new_node, old_node; auto view = make_view(arr, spatial_dimension); auto begin = view.begin(); for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; auto old_ = begin + old_node; auto new_ = begin + new_node; *new_ = *old_; } }; copy(*displacement); copy(*blocked_dofs); if (velocity) copy(*velocity); if (acceleration) copy(*acceleration); if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); // if (external_force) // copy(*external_force); // if (internal_force) // copy(*internal_force); if (displacement_increment) copy(*displacement_increment); copy(getDOFManager().getSolution("displacement")); // this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::afterSolveStep() { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ for (auto & mat : materials) { if (mat->isFiniteDeformation()) mat->computeAllCauchyStresses(_not_ghost); } SolidMechanicsModel::afterSolveStep(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/master_element_info_per_processor.cc b/src/synchronizer/master_element_info_per_processor.cc index 27514c1de..823297ceb 100644 --- a/src/synchronizer/master_element_info_per_processor.cc +++ b/src/synchronizer/master_element_info_per_processor.cc @@ -1,457 +1,451 @@ /** * @file master_element_info_per_processor.cc * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Tue Feb 20 2018 * * @brief Helper class to distribute a mesh * * @section LICENSE * * Copyright (©) 2016-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_iterators.hh" #include "communicator.hh" #include "element_group.hh" #include "element_info_per_processor.hh" #include "element_synchronizer.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MasterElementInfoPerProc::MasterElementInfoPerProc( ElementSynchronizer & synchronizer, UInt message_cnt, UInt root, ElementType type, const MeshPartition & partition) : ElementInfoPerProc(synchronizer, message_cnt, root, type), partition(partition), all_nb_local_element(nb_proc, 0), all_nb_ghost_element(nb_proc, 0), all_nb_element_to_send(nb_proc, 0) { Vector size(5); size(0) = (UInt)type; if (type != _not_defined) { nb_nodes_per_element = Mesh::getNbNodesPerElement(type); nb_element = mesh.getNbElement(type); const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { this->all_nb_local_element[partition_num(el)]++; for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { this->all_nb_ghost_element[*part]++; } this->all_nb_element_to_send[partition_num(el)] += ghost_partition.getNbCols(el) + 1; } /// tag info auto && tag_names = this->mesh.getTagNames(type); this->nb_tags = tag_names.size(); size(4) = nb_tags; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { size(1) = this->all_nb_local_element[p]; size(2) = this->all_nb_ghost_element[p]; size(3) = this->all_nb_element_to_send[p]; AKANTU_DEBUG_INFO( "Sending connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_SIZES) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_SIZES)); } else { this->nb_local_element = this->all_nb_local_element[p]; this->nb_ghost_element = this->all_nb_ghost_element[p]; } } } else { for (UInt p = 0; p < this->nb_proc; ++p) { if (p != this->root) { AKANTU_DEBUG_INFO( "Sending empty connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_SIZES) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_SIZES)); } } } } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizeConnectivities() { const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector> buffers(this->nb_proc); const auto & connectivities = this->mesh.getConnectivity(this->type, _not_ghost); /// copying the local connectivity for (auto && part_conn : zip(partition_num, make_view(connectivities, this->nb_nodes_per_element))) { auto && part = std::get<0>(part_conn); auto && conn = std::get<1>(part_conn); for (UInt i = 0; i < conn.size(); ++i) { buffers[part].push_back(conn[i]); } } /// copying the connectivity of ghost element for (auto && tuple : enumerate(make_view(connectivities, this->nb_nodes_per_element))) { auto && el = std::get<0>(tuple); auto && conn = std::get<1>(tuple); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { UInt proc = *part; for (UInt i = 0; i < conn.size(); ++i) { buffers[proc].push_back(conn[i]); } } } #ifndef AKANTU_NDEBUG for (auto p : arange(this->nb_proc)) { UInt size = this->nb_nodes_per_element * (this->all_nb_local_element[p] + this->all_nb_ghost_element[p]); AKANTU_DEBUG_ASSERT( buffers[p].size() == size, "The connectivity data packed in the buffer are not correct"); } #endif /// send all connectivity and ghost information to all processors std::vector requests; for (auto p : arange(this->nb_proc)) { if (p == this->root) continue; auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_CONNECTIVITY); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } Array & old_nodes = this->getNodesGlobalIds(); /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, buffers[root], all_nb_local_element[root], all_nb_ghost_element[root], type, old_nodes); comm.waitAll(requests); comm.freeCommunicationRequest(requests); } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizePartitions() { const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector> buffers(this->partition.getNbPartition()); /// splitting the partition information to send them to processors Vector count_by_proc(nb_proc, 0); for (UInt el = 0; el < nb_element; ++el) { UInt proc = partition_num(el); buffers[proc].push_back(ghost_partition.getNbCols(el)); UInt i(0); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[proc].push_back(*part); } } for (UInt el = 0; el < nb_element; ++el) { UInt i(0); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[*part].push_back(partition_num(el)); } } #ifndef AKANTU_NDEBUG for (UInt p = 0; p < this->nb_proc; ++p) { AKANTU_DEBUG_ASSERT(buffers[p].size() == (this->all_nb_ghost_element[p] + this->all_nb_element_to_send[p]), "Data stored in the buffer are most probably wrong"); } #endif std::vector requests; /// last data to compute the communication scheme for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->root) continue; auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_PARTITIONS); AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } if (Mesh::getSpatialDimension(this->type) == this->mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); this->fillCommunicationScheme(buffers[this->rank]); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeTags() { AKANTU_DEBUG_IN(); if (this->nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } - UInt mesh_data_sizes_buffer_length; - /// tag info auto tag_names = mesh.getTagNames(type); // Make sure the tags are sorted (or at least not in random order), // because they come from a map !! std::sort(tag_names.begin(), tag_names.end()); // Sending information about the tags in mesh_data: name, data type and // number of components of the underlying array associated to the current // type DynamicCommunicationBuffer mesh_data_sizes_buffer; for (auto && tag_name : tag_names) { mesh_data_sizes_buffer << tag_name; mesh_data_sizes_buffer << mesh.getTypeCode(tag_name); mesh_data_sizes_buffer << mesh.getNbComponent(tag_name, type); } - mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.size(); AKANTU_DEBUG_INFO( "Broadcasting the size of the information about the mesh data tags: (" - << mesh_data_sizes_buffer_length << ")."); - comm.broadcast(mesh_data_sizes_buffer_length, root); + << mesh_data_sizes_buffer.size() << ")."); AKANTU_DEBUG_INFO( "Broadcasting the information about the mesh data tags, addr " << (void *)mesh_data_sizes_buffer.storage()); - if (mesh_data_sizes_buffer_length != 0) - comm.broadcast(mesh_data_sizes_buffer, root); + comm.broadcast(mesh_data_sizes_buffer, root); - if (mesh_data_sizes_buffer_length != 0) { - // Sending the actual data to each processor - std::vector buffers(nb_proc); - // Loop over each tag for the current type - for (auto && tag_name : tag_names) { - // Type code of the current tag (i.e. the tag named *names_it) - this->fillTagBuffer(buffers, tag_name); - } + if (mesh_data_sizes_buffer.size() == 0) + return; - std::vector requests; - for (UInt p = 0; p < nb_proc; ++p) { - if (p == root) - continue; + // Sending the actual data to each processor + std::vector buffers(nb_proc); + // Loop over each tag for the current type + for (auto && tag_name : tag_names) { + // Type code of the current tag (i.e. the tag named *names_it) + this->fillTagBuffer(buffers, tag_name); + } - auto && tag = - Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA); - AKANTU_DEBUG_INFO("Sending " << buffers[p].size() - << " bytes of mesh data to proc " << p - << " TAG(" << tag << ")"); + std::vector requests; + for (UInt p = 0; p < nb_proc; ++p) { + if (p == root) + continue; - requests.push_back(comm.asyncSend(buffers[p], p, tag)); - } + auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA); + AKANTU_DEBUG_INFO("Sending " << buffers[p].size() + << " bytes of mesh data to proc " << p + << " TAG(" << tag << ")"); - // Loop over each tag for the current type - for (auto && tag_name : tag_names) { - // Reinitializing the mesh data on the master - this->fillMeshData(buffers[root], tag_name, - mesh.getTypeCode(tag_name), - mesh.getNbComponent(tag_name, type)); - } + requests.push_back(comm.asyncSend(buffers[p], p, tag)); + } - comm.waitAll(requests); - comm.freeCommunicationRequest(requests); - requests.clear(); + // Loop over each tag for the current type + for (auto && tag_name : tag_names) { + // Reinitializing the mesh data on the master + this->fillMeshData(buffers[root], tag_name, mesh.getTypeCode(tag_name), + mesh.getNbComponent(tag_name, type)); } + comm.waitAll(requests); + comm.freeCommunicationRequest(requests); + requests.clear(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MasterElementInfoPerProc::fillTagBufferTemplated( std::vector & buffers, const std::string & tag_name) { const auto & data = mesh.getElementalDataArray(tag_name, type); const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); // Not possible to use the iterator because it potentially triggers the // creation of complex // type templates (such as akantu::Vector< std::vector > which don't // implement the right interface // (e.g. operator<< in that case). // typename Array::template const_iterator< Vector > data_it = // data.begin(data.getNbComponent()); // typename Array::template const_iterator< Vector > data_end = // data.end(data.getNbComponent()); const T * data_it = data.storage(); const T * data_end = data.storage() + data.size() * data.getNbComponent(); const UInt * part = partition_num.storage(); /// copying the data, element by element for (; data_it != data_end; ++part) { for (UInt j(0); j < data.getNbComponent(); ++j, ++data_it) { buffers[*part] << *data_it; } } data_it = data.storage(); /// copying the data for the ghost element for (UInt el(0); data_it != data_end; data_it += data.getNbComponent(), ++el) { auto it = ghost_partition.begin(el); auto end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; for (UInt j(0); j < data.getNbComponent(); ++j) { buffers[proc] << data_it[j]; } } } } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::fillTagBuffer( std::vector & buffers, const std::string & tag_name) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ this->fillTagBufferTemplated(buffers, \ tag_name); \ break; \ } MeshDataTypeCode data_type_code = mesh.getTypeCode(tag_name); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); std::vector buffers(nb_proc); using ElementToGroup = std::vector>; ElementToGroup element_to_group(nb_element); for (auto & eg : ElementGroupsIterable(mesh)) { const auto & name = eg.getName(); for (const auto & element : eg.getElements(type, _not_ghost)) { element_to_group[element].push_back(name); } auto eit = eg.begin(type, _not_ghost); if (eit != eg.end(type, _not_ghost)) const_cast &>(eg.getElements(type)).empty(); } const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); /// copying the data, element by element for (auto && pair : zip(partition_num, element_to_group)) { buffers[std::get<0>(pair)] << std::get<1>(pair); } /// copying the data for the ghost element for (auto && pair : enumerate(element_to_group)) { auto && el = std::get<0>(pair); auto it = ghost_partition.begin(el); auto end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; buffers[proc] << std::get<1>(pair); } } std::vector requests; for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->rank) continue; auto && tag = Tag::genTag(this->rank, p, Tag::_ELEMENT_GROUP); AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } this->fillElementGroupsFromBuffer(buffers[this->rank]); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/node_info_per_processor.cc b/src/synchronizer/node_info_per_processor.cc index 71fbdbf55..e8f6ed990 100644 --- a/src/synchronizer/node_info_per_processor.cc +++ b/src/synchronizer/node_info_per_processor.cc @@ -1,866 +1,851 @@ /** * @file node_info_per_processor.cc * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Wed Nov 08 2017 * * @brief Please type the brief for file: Helper classes to create the * distributed synchronizer and distribute a mesh * * @section LICENSE * * Copyright (©) 2016-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 "node_info_per_processor.hh" #include "communicator.hh" #include "node_group.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NodeInfoPerProc::NodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root) : MeshAccessor(synchronizer.getMesh()), synchronizer(synchronizer), comm(synchronizer.getCommunicator()), rank(comm.whoAmI()), nb_proc(comm.getNbProc()), root(root), mesh(synchronizer.getMesh()), spatial_dimension(synchronizer.getMesh().getSpatialDimension()), message_count(message_cnt) {} /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::synchronize() { synchronizeNodes(); synchronizeTypes(); synchronizeGroups(); synchronizePeriodicity(); synchronizeTags(); } /* -------------------------------------------------------------------------- */ template void NodeInfoPerProc::fillNodeGroupsFromBuffer(CommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); std::vector> node_to_group; buffer >> node_to_group; AKANTU_DEBUG_ASSERT(node_to_group.size() == mesh.getNbGlobalNodes(), "Not the good amount of nodes where transmitted"); const auto & global_nodes = mesh.getGlobalNodesIds(); for (auto && data : enumerate(global_nodes)) { for (const auto & node : node_to_group[std::get<1>(data)]) { mesh.getNodeGroup(node).add(std::get<0>(data), false); } } for (auto && ng_data : mesh.iterateNodeGroups()) { ng_data.optimize(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillNodesType() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); auto & nodes_flags = this->getNodesFlags(); Array nodes_set(nb_nodes); nodes_set.set(0); enum NodeSet { NORMAL_SET = 1, GHOST_SET = 2, }; Array already_seen(nb_nodes, 1, false); for (auto gt : ghost_types) { UInt set = NORMAL_SET; if (gt == _ghost) set = GHOST_SET; already_seen.set(false); for (auto && type : mesh.elementTypes(_all_dimensions, gt, _ek_not_defined)) { const auto & connectivity = mesh.getConnectivity(type, gt); for (auto & conn : make_view(connectivity, connectivity.getNbComponent())) { for (UInt n = 0; n < conn.size(); ++n) { AKANTU_DEBUG_ASSERT(conn(n) < nb_nodes, "Node " << conn(n) << " bigger than number of nodes " << nb_nodes); if (!already_seen(conn(n))) { nodes_set(conn(n)) += set; already_seen(conn(n)) = true; } } } } } nodes_flags.resize(nb_nodes); for (UInt i = 0; i < nb_nodes; ++i) { if (nodes_set(i) == NORMAL_SET) nodes_flags(i) = NodeFlag::_normal; else if (nodes_set(i) == GHOST_SET) nodes_flags(i) = NodeFlag::_pure_ghost; else if (nodes_set(i) == (GHOST_SET + NORMAL_SET)) nodes_flags(i) = NodeFlag::_master; else { AKANTU_EXCEPTION("Gni ?"); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillCommunicationScheme(const Array & master_info) { AKANTU_DEBUG_IN(); Communications & communications = this->synchronizer.getCommunications(); { // send schemes auto it = master_info.begin_reinterpret(2, master_info.size() / 2); auto end = master_info.end_reinterpret(2, master_info.size() / 2); std::map> send_array_per_proc; for (; it != end; ++it) { const Vector & send_info = *it; send_array_per_proc[send_info(0)].push_back(send_info(1)); } for (auto & send_schemes : send_array_per_proc) { auto & scheme = communications.createSendScheme(send_schemes.first); auto & sends = send_schemes.second; std::sort(sends.begin(), sends.end()); std::transform(sends.begin(), sends.end(), sends.begin(), [this](UInt g) -> UInt { return mesh.getNodeLocalId(g); }); scheme.copy(sends); } } { // receive schemes std::map> recv_array_per_proc; for (auto node : arange(mesh.getNbNodes())) { if (mesh.isSlaveNode(node)) { recv_array_per_proc[mesh.getNodePrank(node)].push_back( mesh.getNodeGlobalId(node)); } } for (auto & recv_schemes : recv_array_per_proc) { auto & scheme = communications.createRecvScheme(recv_schemes.first); auto & recvs = recv_schemes.second; std::sort(recvs.begin(), recvs.end()); std::transform(recvs.begin(), recvs.end(), recvs.begin(), [this](UInt g) -> UInt { return mesh.getNodeLocalId(g); }); scheme.copy(recvs); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillPeriodicPairs(const Array & global_pairs, std::vector & missing_nodes) { this->wipePeriodicInfo(); auto & nodes_flags = this->getNodesFlags(); auto checkIsLocal = [&](auto && global_node) { auto && node = mesh.getNodeLocalId(global_node); if (node == UInt(-1)) { auto & global_nodes = this->getNodesGlobalIds(); node = global_nodes.size(); global_nodes.push_back(global_node); nodes_flags.push_back(NodeFlag::_pure_ghost); missing_nodes.push_back(global_node); std::cout << "Missing node " << node << std::endl; } return node; }; for (auto && pairs : make_view(global_pairs, 2)) { UInt slave = checkIsLocal(pairs(0)); UInt master = checkIsLocal(pairs(1)); this->addPeriodicSlave(slave, master); } this->markMeshPeriodic(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::receiveMissingPeriodic( DynamicCommunicationBuffer & buffer) { auto & nodes = this->getNodes(); Communications & communications = this->synchronizer.getCommunications(); std::size_t nb_nodes; buffer >> nb_nodes; for (auto _[[gnu::unused]] : arange(nb_nodes)) { Vector pos(spatial_dimension); Int prank; buffer >> pos; buffer >> prank; UInt node = nodes.size(); this->setNodePrank(node, prank); nodes.push_back(pos); auto & scheme = communications.createRecvScheme(prank); scheme.push_back(node); } while (buffer.getLeftToUnpack() != 0) { Int prank; UInt gnode; buffer >> gnode; buffer >> prank; auto node = mesh.getNodeLocalId(gnode); AKANTU_DEBUG_ASSERT(node != UInt(-1), "I cannot send the node " << gnode << " to proc " << prank << " because it is note a local node"); auto & scheme = communications.createSendScheme(prank); scheme.push_back(node); } } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillNodalData(DynamicCommunicationBuffer & buffer, std::string tag_name) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, _, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ auto & nodal_data = \ mesh.getNodalData(tag_name); \ nodal_data.resize(mesh.getNbNodes()); \ for (auto && data : make_view(nodal_data)) { \ buffer >> data; \ } \ break; \ } MeshDataTypeCode data_type_code = mesh.getTypeCode(tag_name, MeshDataType::_nodal); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ MasterNodeInfoPerProc::MasterNodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root) : NodeInfoPerProc(synchronizer, message_cnt, root), all_nodes(0, synchronizer.getMesh().getSpatialDimension()) { UInt nb_global_nodes = this->mesh.getNbGlobalNodes(); this->comm.broadcast(nb_global_nodes, this->root); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeNodes() { this->nodes_per_proc.resize(nb_proc); this->nb_nodes_per_proc.resize(nb_proc); Array local_nodes(0, spatial_dimension); Array & nodes = this->getNodes(); all_nodes.copy(nodes); nodes_pranks.resize(nodes.size(), UInt(-1)); for (UInt p = 0; p < nb_proc; ++p) { UInt nb_nodes = 0; // UInt * buffer; Array * nodes_to_send; Array & nodespp = nodes_per_proc[p]; if (p != root) { nodes_to_send = new Array(0, spatial_dimension); AKANTU_DEBUG_INFO("Receiving number of nodes from proc " << p << " " << Tag::genTag(p, 0, Tag::_NB_NODES)); comm.receive(nb_nodes, p, Tag::genTag(p, 0, Tag::_NB_NODES)); nodespp.resize(nb_nodes); this->nb_nodes_per_proc(p) = nb_nodes; AKANTU_DEBUG_INFO("Receiving list of nodes from proc " << p << " " << Tag::genTag(p, 0, Tag::_NODES)); comm.receive(nodespp, p, Tag::genTag(p, 0, Tag::_NODES)); } else { Array & local_ids = this->getNodesGlobalIds(); this->nb_nodes_per_proc(p) = local_ids.size(); nodespp.copy(local_ids); nodes_to_send = &local_nodes; } Array::const_scalar_iterator it = nodespp.begin(); Array::const_scalar_iterator end = nodespp.end(); /// get the coordinates for the selected nodes for (; it != end; ++it) { Vector coord(nodes.storage() + spatial_dimension * *it, spatial_dimension); nodes_to_send->push_back(coord); } if (p != root) { /// send them for distant processors AKANTU_DEBUG_INFO("Sending coordinates to proc " << p << " " << Tag::genTag(this->rank, 0, Tag::_COORDINATES)); comm.send(*nodes_to_send, p, Tag::genTag(this->rank, 0, Tag::_COORDINATES)); delete nodes_to_send; } } /// construct the local nodes coordinates nodes.copy(local_nodes); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeTypes() { // > std::multimap> nodes_to_proc; std::vector> nodes_flags_per_proc(nb_proc); std::vector> nodes_prank_per_proc(nb_proc); if (mesh.isPeriodic()) all_periodic_flags.copy(this->getNodesFlags()); // arrays containing pairs of (proc, node) std::vector> nodes_to_send_per_proc(nb_proc); for (UInt p = 0; p < nb_proc; ++p) { nodes_flags_per_proc[p].resize(nb_nodes_per_proc(p), NodeFlag(0xFF)); nodes_prank_per_proc[p].resize(nb_nodes_per_proc(p), -1); } this->fillNodesType(); auto is_master = [](auto && flag) { return (flag & NodeFlag::_shared_mask) == NodeFlag::_master; }; auto is_local = [](auto && flag) { return (flag & NodeFlag::_shared_mask) == NodeFlag::_normal; }; for (auto p : arange(nb_proc)) { auto & nodes_flags = nodes_flags_per_proc[p]; if (p != root) { AKANTU_DEBUG_INFO( "Receiving first nodes types from proc " << p << " " << Tag::genTag(this->rank, this->message_count, Tag::_NODES_TYPE)); comm.receive(nodes_flags, p, Tag::genTag(p, 0, Tag::_NODES_TYPE)); } else { nodes_flags.copy(this->getNodesFlags()); } // stack all processors claiming to be master for a node for (auto local_node : arange(nb_nodes_per_proc(p))) { auto global_node = nodes_per_proc[p](local_node); if (is_master(nodes_flags(local_node))) { nodes_to_proc.insert( std::make_pair(global_node, std::make_pair(p, local_node))); } else if (is_local(nodes_flags(local_node))) { nodes_pranks[global_node] = p; } } } for (auto i : arange(mesh.getNbGlobalNodes())) { auto it_range = nodes_to_proc.equal_range(i); if (it_range.first == nodes_to_proc.end() || it_range.first->first != i) continue; // pick the first processor out of the multi-map as the actual master UInt master_proc = (it_range.first)->second.first; nodes_pranks[i] = master_proc; for (auto && data : range(it_range.first, it_range.second)) { auto proc = data.second.first; auto node = data.second.second; if (proc != master_proc) { // store the info on all the slaves for a given master nodes_flags_per_proc[proc](node) = NodeFlag::_slave; nodes_to_send_per_proc[master_proc].push_back(proc); nodes_to_send_per_proc[master_proc].push_back(i); } } } /// Fills the nodes prank per proc for (auto && data : zip(arange(nb_proc), nodes_per_proc, nodes_prank_per_proc, nodes_flags_per_proc)) { for (auto && node_data : zip(std::get<1>(data), std::get<2>(data), std::get<3>(data))) { if (std::get<2>(node_data) == NodeFlag::_normal) { std::get<1>(node_data) = std::get<0>(data); } else { std::get<1>(node_data) = nodes_pranks(std::get<0>(node_data)); } } } std::vector requests_send_type; std::vector requests_send_master_info; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Sending nodes types to proc " << p << " " << Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)); requests_send_type.push_back( comm.asyncSend(nodes_flags_per_proc[p], p, Tag::genTag(this->rank, 0, Tag::_NODES_TYPE))); requests_send_type.push_back( comm.asyncSend(nodes_prank_per_proc[p], p, Tag::genTag(this->rank, 2, Tag::_NODES_TYPE))); auto & nodes_to_send = nodes_to_send_per_proc[p]; AKANTU_DEBUG_INFO("Sending nodes master info to proc " << p << " " << Tag::genTag(this->rank, 1, Tag::_NODES_TYPE)); requests_send_master_info.push_back(comm.asyncSend( nodes_to_send, p, Tag::genTag(this->rank, 1, Tag::_NODES_TYPE))); } else { this->getNodesFlags().copy(nodes_flags_per_proc[p]); for (auto && data : enumerate(nodes_prank_per_proc[p])) { auto node = std::get<0>(data); if (not(mesh.isMasterNode(node) or mesh.isLocalNode(node))) { this->setNodePrank(node, std::get<1>(data)); } } this->fillCommunicationScheme(nodes_to_send_per_proc[root]); } } comm.waitAll(requests_send_type); comm.freeCommunicationRequest(requests_send_type); requests_send_type.clear(); comm.waitAll(requests_send_master_info); comm.freeCommunicationRequest(requests_send_master_info); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); UInt nb_total_nodes = mesh.getNbGlobalNodes(); DynamicCommunicationBuffer buffer; using NodeToGroup = std::vector>; NodeToGroup node_to_group; node_to_group.resize(nb_total_nodes); GroupManager::const_node_group_iterator ngi = mesh.node_group_begin(); GroupManager::const_node_group_iterator nge = mesh.node_group_end(); for (; ngi != nge; ++ngi) { NodeGroup & ng = *(ngi->second); std::string name = ngi->first; NodeGroup::const_node_iterator nit = ng.begin(); NodeGroup::const_node_iterator nend = ng.end(); for (; nit != nend; ++nit) { node_to_group[*nit].push_back(name); } nit = ng.begin(); if (nit != nend) ng.empty(); } buffer << node_to_group; std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == this->rank) continue; AKANTU_DEBUG_INFO("Sending node groups to proc " << p << " " << Tag::genTag(this->rank, p, Tag::_NODE_GROUP)); requests.push_back(comm.asyncSend( buffer, p, Tag::genTag(this->rank, p, Tag::_NODE_GROUP))); } this->fillNodeGroupsFromBuffer(buffer); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizePeriodicity() { bool is_periodic = mesh.isPeriodic(); comm.broadcast(is_periodic, root); if (not is_periodic) return; std::vector requests; std::vector> periodic_info_to_send_per_proc; for (auto p : arange(nb_proc)) { periodic_info_to_send_per_proc.emplace_back(0, 2); auto && periodic_info = periodic_info_to_send_per_proc.back(); for (UInt proc_local_node : arange(nb_nodes_per_proc(p))) { UInt global_node = nodes_per_proc[p](proc_local_node); if ((all_periodic_flags[global_node] & NodeFlag::_periodic_mask) == NodeFlag::_periodic_slave) { periodic_info.push_back( Vector{global_node, mesh.getPeriodicMaster(global_node)}); } } if (p == root) continue; auto && tag = Tag::genTag(this->rank, p, Tag::_PERIODIC_SLAVES); AKANTU_DEBUG_INFO("Sending periodic info to proc " << p << " " << tag); requests.push_back(comm.asyncSend(periodic_info, p, tag)); } CommunicationStatus status; std::vector buffers(nb_proc); std::vector> proc_missings(nb_proc); auto nodes_it = all_nodes.begin(spatial_dimension); for (UInt p = 0; p < nb_proc; ++p) { auto & proc_missing = proc_missings[p]; if (p != root) { auto && tag = Tag::genTag(p, 0, Tag::_PERIODIC_NODES); comm.probe(p, tag, status); proc_missing.resize(status.size()); comm.receive(proc_missing, p, tag); } else { fillPeriodicPairs(periodic_info_to_send_per_proc[root], proc_missing); } auto & buffer = buffers[p]; buffer.reserve((spatial_dimension * sizeof(Real) + sizeof(Int)) * proc_missing.size()); buffer << proc_missing.size(); for (auto && node : proc_missing) { buffer << *(nodes_it + node); buffer << nodes_pranks(node); } } for (UInt p = 0; p < nb_proc; ++p) { for (auto && node : proc_missings[p]) { auto & buffer = buffers[nodes_pranks(node)]; buffer << node; buffer << p; } } for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { auto && tag_send = Tag::genTag(p, 1, Tag::_PERIODIC_NODES); requests.push_back(comm.asyncSend(buffers[p], p, tag_send)); } else { receiveMissingPeriodic(buffers[p]); } } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::fillTagBuffers( std::vector & buffers, const std::string & tag_name) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, _, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ auto & nodal_data = \ mesh.getNodalData(tag_name); \ for (auto && data : enumerate(nodes_per_proc)) { \ auto proc = std::get<0>(data); \ auto & nodes = std::get<1>(data); \ auto & buffer = buffers[proc]; \ for (auto & node : nodes) { \ for (auto i : arange(nodal_data.getNbComponent())) { \ buffer << nodal_data(node, i); \ } \ } \ } \ break; \ } MeshDataTypeCode data_type_code = mesh.getTypeCode(tag_name, MeshDataType::_nodal); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } // namespace akantu /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeTags() { /// tag info auto tag_names = mesh.getTagNames(); DynamicCommunicationBuffer tags_buffer; for (auto && tag_name : tag_names) { tags_buffer << tag_name; tags_buffer << mesh.getTypeCode(tag_name, MeshDataType::_nodal); tags_buffer << mesh.getNbComponent(tag_name); } - UInt tags_buffer_length = tags_buffer.size(); AKANTU_DEBUG_INFO( - "Broadcasting the size of the information about the mesh data tags: (" - << tags_buffer_length << ")."); - comm.broadcast(tags_buffer_length, root); - - if (tags_buffer_length != 0) - comm.broadcast(tags_buffer, root); + "Broadcasting the information about the nodes mesh data tags: (" + << tags_buffer.size() << ")."); + comm.broadcast(tags_buffer, root); for (auto && tag_data : enumerate(tag_names)) { auto tag_count = std::get<0>(tag_data); auto & tag_name = std::get<1>(tag_data); std::vector buffers; std::vector requests; buffers.resize(nb_proc); fillTagBuffers(buffers, tag_name); for (auto && data : enumerate(buffers)) { auto && proc = std::get<0>(data); auto & buffer = std::get<1>(data); if (proc == root) { fillNodalData(buffer, tag_name); } else { auto && tag = Tag::genTag(this->rank, tag_count, Tag::_MESH_DATA); requests.push_back(comm.asyncSend(buffer, proc, tag)); } } comm.waitAll(requests); comm.freeCommunicationRequest(requests); } } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ SlaveNodeInfoPerProc::SlaveNodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root) : NodeInfoPerProc(synchronizer, message_cnt, root) { UInt nb_global_nodes = 0; comm.broadcast(nb_global_nodes, root); this->setNbGlobalNodes(nb_global_nodes); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeNodes() { AKANTU_DEBUG_INFO("Sending list of nodes to proc " << root << " " << Tag::genTag(this->rank, 0, Tag::_NB_NODES) << " " << Tag::genTag(this->rank, 0, Tag::_NODES)); Array & local_ids = this->getNodesGlobalIds(); Array & nodes = this->getNodes(); UInt nb_nodes = local_ids.size(); comm.send(nb_nodes, root, Tag::genTag(this->rank, 0, Tag::_NB_NODES)); comm.send(local_ids, root, Tag::genTag(this->rank, 0, Tag::_NODES)); /* --------<<<<-COORDINATES---------------------------------------------- */ nodes.resize(nb_nodes); AKANTU_DEBUG_INFO("Receiving coordinates from proc " << root << " " << Tag::genTag(root, 0, Tag::_COORDINATES)); comm.receive(nodes, root, Tag::genTag(root, 0, Tag::_COORDINATES)); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeTypes() { this->fillNodesType(); auto & nodes_types = this->getNodesFlags(); AKANTU_DEBUG_INFO("Sending first nodes types to proc " << root << "" << Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)); comm.send(nodes_types, root, Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)); AKANTU_DEBUG_INFO("Receiving nodes types from proc " << root << " " << Tag::genTag(root, 0, Tag::_NODES_TYPE)); comm.receive(nodes_types, root, Tag::genTag(root, 0, Tag::_NODES_TYPE)); Array nodes_prank(nodes_types.size()); comm.receive(nodes_prank, root, Tag::genTag(root, 2, Tag::_NODES_TYPE)); for (auto && data : enumerate(nodes_prank)) { auto node = std::get<0>(data); if (not(mesh.isMasterNode(node) or mesh.isLocalNode(node))) { this->setNodePrank(node, std::get<1>(data)); } } AKANTU_DEBUG_INFO("Receiving nodes master info from proc " << root << " " << Tag::genTag(root, 1, Tag::_NODES_TYPE)); CommunicationStatus status; comm.probe(root, Tag::genTag(root, 1, Tag::_NODES_TYPE), status); Array nodes_master_info(status.size()); if (nodes_master_info.size() > 0) comm.receive(nodes_master_info, root, Tag::genTag(root, 1, Tag::_NODES_TYPE)); this->fillCommunicationScheme(nodes_master_info); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Receiving node groups from proc " << root << " " << Tag::genTag(root, this->rank, Tag::_NODE_GROUP)); - CommunicationStatus status; - comm.probe(root, Tag::genTag(root, this->rank, Tag::_NODE_GROUP), - status); - - CommunicationBuffer buffer(status.size()); + DynamicCommunicationBuffer buffer; comm.receive(buffer, root, Tag::genTag(root, this->rank, Tag::_NODE_GROUP)); this->fillNodeGroupsFromBuffer(buffer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizePeriodicity() { bool is_periodic; comm.broadcast(is_periodic, root); if (not is_periodic) return; CommunicationStatus status; auto && tag = Tag::genTag(root, this->rank, Tag::_PERIODIC_SLAVES); comm.probe(root, tag, status); Array periodic_info(status.size() / 2, 2); comm.receive(periodic_info, root, tag); std::vector proc_missing; fillPeriodicPairs(periodic_info, proc_missing); auto && tag_missing_request = Tag::genTag(this->rank, 0, Tag::_PERIODIC_NODES); comm.send(proc_missing, root, tag_missing_request); DynamicCommunicationBuffer buffer; auto && tag_missing = Tag::genTag(this->rank, 1, Tag::_PERIODIC_NODES); comm.receive(buffer, root, tag_missing); receiveMissingPeriodic(buffer); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeTags() { - UInt tags_buffer_length{0}; - comm.broadcast(tags_buffer_length, root); - if (tags_buffer_length == 0) { - return; - } - - CommunicationBuffer tags_buffer; - tags_buffer.resize(tags_buffer_length); + DynamicCommunicationBuffer tags_buffer; comm.broadcast(tags_buffer, root); std::vector tag_names; while (tags_buffer.getLeftToUnpack() > 0) { std::string name; MeshDataTypeCode code; UInt nb_components; tags_buffer >> name; tags_buffer >> code; tags_buffer >> nb_components; mesh.registerNodalData(name, nb_components, code); tag_names.push_back(name); } for (auto && tag_data : enumerate(tag_names)) { auto tag_count = std::get<0>(tag_data); auto & tag_name = std::get<1>(tag_data); DynamicCommunicationBuffer buffer; auto && tag = Tag::genTag(this->root, tag_count, Tag::_MESH_DATA); comm.receive(buffer, this->root, tag); fillNodalData(buffer, tag_name); } } } // namespace akantu diff --git a/src/synchronizer/node_synchronizer.cc b/src/synchronizer/node_synchronizer.cc index fcb7010ba..43fe4a2b3 100644 --- a/src/synchronizer/node_synchronizer.cc +++ b/src/synchronizer/node_synchronizer.cc @@ -1,221 +1,226 @@ /** -* @file node_synchronizer.cc -* -* @author Nicolas Richart -* -* @date creation: Fri Jun 18 2010 -* @date last modification: Wed Nov 15 2017 -* -* @brief Implementation of the node synchronizer -* -* @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 . -* -*/ + * @file node_synchronizer.cc + * + * @author Nicolas Richart + * + * @date creation: Fri Jun 18 2010 + * @date last modification: Wed Nov 15 2017 + * + * @brief Implementation of the node synchronizer + * + * @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 "node_synchronizer.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NodeSynchronizer::NodeSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id, const bool register_to_event_manager, EventHandlerPriority event_priority) : SynchronizerImpl(mesh.getCommunicator(), id, memory_id), mesh(mesh) { AKANTU_DEBUG_IN(); if (register_to_event_manager) { this->mesh.registerEventHandler(*this, event_priority); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NodeSynchronizer::~NodeSynchronizer() = default; /* -------------------------------------------------------------------------- */ Int NodeSynchronizer::getRank(const UInt & node) const { return this->mesh.getNodePrank(node); } /* -------------------------------------------------------------------------- */ void NodeSynchronizer::onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent &) { std::map> nodes_per_proc; // recreates fully the schemes due to changes of global ids // \TODO add an event to handle global id changes for (auto && data : communications.iterateSchemes(_recv)) { auto & scheme = data.second; scheme.resize(0); } for (auto && local_id : arange(mesh.getNbNodes())) { if (not mesh.isSlaveNode(local_id)) continue; // local, master or pure ghost auto global_id = mesh.getNodeGlobalId(local_id); auto proc = mesh.getNodePrank(local_id); + AKANTU_DEBUG_ASSERT( + proc != -1, + "The node " << local_id << " does not have a valid associated prank"); nodes_per_proc[proc].push_back(global_id); auto & scheme = communications.createScheme(proc, _recv); scheme.push_back(local_id); } std::vector send_requests; for (auto && pair : communications.iterateSchemes(_recv)) { auto proc = pair.first; + AKANTU_DEBUG_ASSERT(proc != UInt(-1), + "For real I should send something to proc -1"); // if proc not in nodes_per_proc this should insert an empty array to send send_requests.push_back(communicator.asyncSend( nodes_per_proc[proc], proc, Tag::genTag(rank, proc, 0xcafe))); } for (auto && data : communications.iterateSchemes(_send)) { auto proc = data.first; auto & scheme = data.second; CommunicationStatus status; auto tag = Tag::genTag(proc, rank, 0xcafe); communicator.probe(proc, tag, status); scheme.resize(status.size()); communicator.receive(scheme, proc, tag); std::transform(scheme.begin(), scheme.end(), scheme.begin(), [&](auto & gnode) { return mesh.getNodeLocalId(gnode); }); } // communicator.receiveAnyNumber( // send_requests, // [&](auto && proc, auto && nodes) { // auto & scheme = communications.createScheme(proc, _send); // scheme.resize(nodes.size()); // for (auto && data : enumerate(nodes)) { // auto global_id = std::get<1>(data); // auto local_id = mesh.getNodeLocalId(global_id); // AKANTU_DEBUG_ASSERT(local_id != UInt(-1), // "The global node " << global_id // << "is not known on rank " // << rank); // scheme[std::get<0>(data)] = local_id; // } // }, // Tag::genTag(rank, count, 0xcafe)); // ++count; communicator.waitAll(send_requests); communicator.freeCommunicationRequest(send_requests); } /* -------------------------------------------------------------------------- */ UInt NodeSynchronizer::sanityCheckDataSize(const Array & nodes, const SynchronizationTag & tag, bool from_comm_desc) const { UInt size = SynchronizerImpl::sanityCheckDataSize(nodes, tag, from_comm_desc); // global id if (tag != _gst_giu_global_conn) { size += sizeof(UInt) * nodes.size(); } // flag size += sizeof(NodeFlag) * nodes.size(); // positions size += mesh.getSpatialDimension() * sizeof(Real) * nodes.size(); return size; } /* -------------------------------------------------------------------------- */ void NodeSynchronizer::packSanityCheckData( CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const { auto dim = mesh.getSpatialDimension(); for (auto && node : nodes) { if (tag != _gst_giu_global_conn) { buffer << mesh.getNodeGlobalId(node); } buffer << mesh.getNodeFlag(node); buffer << Vector(mesh.getNodes().begin(dim)[node]); } } /* -------------------------------------------------------------------------- */ void NodeSynchronizer::unpackSanityCheckData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag, UInt proc, UInt rank) const { auto dim = mesh.getSpatialDimension(); // std::set skip_conn_tags{_gst_smmc_facets_conn, // _gst_giu_global_conn}; // bool is_skip_tag_conn = skip_conn_tags.find(tag) != skip_conn_tags.end(); auto periodic = [&](auto && flag) { return flag & NodeFlag::_periodic_mask; }; auto distrib = [&](auto && flag) { return flag & NodeFlag::_shared_mask; }; for (auto && node : nodes) { if (tag != _gst_giu_global_conn) { UInt global_id; buffer >> global_id; AKANTU_DEBUG_ASSERT(global_id == mesh.getNodeGlobalId(node), "The nodes global ids do not match: " << global_id << " != " << mesh.getNodeGlobalId(node)); } NodeFlag flag; buffer >> flag; AKANTU_DEBUG_ASSERT( (periodic(flag) == periodic(mesh.getNodeFlag(node))) and (((distrib(flag) == NodeFlag::_master) and (distrib(mesh.getNodeFlag(node)) == NodeFlag::_slave)) or // master to slave ((distrib(flag) == NodeFlag::_slave) and (distrib(mesh.getNodeFlag(node)) == NodeFlag::_master)) or // reverse comm slave to master (distrib(mesh.getNodeFlag(node)) == NodeFlag::_pure_ghost or // pure ghost nodes distrib(flag) == NodeFlag::_pure_ghost)), - "The node flags do not make sense: " - << std::hex << "0x" << flag << " and 0x" << mesh.getNodeFlag(node)); + "The node flags: " + << flag << " and " << mesh.getNodeFlag(node)); Vector pos_remote(dim); buffer >> pos_remote; Vector pos(mesh.getNodes().begin(dim)[node]); auto dist = pos_remote.distance(pos); if (not Math::are_float_equal(dist, 0.)) { AKANTU_EXCEPTION("Unpacking an unknown value for the node " << node << "(position " << pos << " != buffer " << pos_remote << ") [" << dist << "] - tag: " << tag << " comm from " << proc << " to " << rank); } } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/slave_element_info_per_processor.cc b/src/synchronizer/slave_element_info_per_processor.cc index e38400012..7756b1363 100644 --- a/src/synchronizer/slave_element_info_per_processor.cc +++ b/src/synchronizer/slave_element_info_per_processor.cc @@ -1,198 +1,192 @@ /** * @file slave_element_info_per_processor.cc * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Tue Nov 07 2017 * * @brief Helper class to distribute a mesh * * @section LICENSE * * Copyright (©) 2016-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 "communicator.hh" #include "element_info_per_processor.hh" #include "element_synchronizer.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SlaveElementInfoPerProc::SlaveElementInfoPerProc( ElementSynchronizer & synchronizer, UInt message_cnt, UInt root) : ElementInfoPerProc(synchronizer, message_cnt, root, _not_defined) { Vector size(5); comm.receive(size, this->root, Tag::genTag(this->root, this->message_count, Tag::_SIZES)); this->type = (ElementType)size[0]; this->nb_local_element = size[1]; this->nb_ghost_element = size[2]; this->nb_element_to_receive = size[3]; this->nb_tags = size[4]; if (this->type != _not_defined) this->nb_nodes_per_element = Mesh::getNbNodesPerElement(type); } /* -------------------------------------------------------------------------- */ bool SlaveElementInfoPerProc::needSynchronize() { return this->type != _not_defined; } /* -------------------------------------------------------------------------- */ void SlaveElementInfoPerProc::synchronizeConnectivities() { Array local_connectivity( (this->nb_local_element + this->nb_ghost_element) * this->nb_nodes_per_element); AKANTU_DEBUG_INFO("Receiving connectivities from proc " << root); comm.receive( local_connectivity, this->root, Tag::genTag(this->root, this->message_count, Tag::_CONNECTIVITY)); auto & old_nodes = this->getNodesGlobalIds(); AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(this->mesh, local_connectivity, this->nb_local_element, this->nb_ghost_element, this->type, old_nodes); } /* -------------------------------------------------------------------------- */ void SlaveElementInfoPerProc::synchronizePartitions() { Array local_partitions(this->nb_element_to_receive + this->nb_ghost_element * 2); AKANTU_DEBUG_INFO("Receiving partition informations from proc " << root); this->comm.receive(local_partitions, this->root, Tag::genTag(root, this->message_count, Tag::_PARTITIONS)); if (Mesh::getSpatialDimension(this->type) == this->mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); this->fillCommunicationScheme(local_partitions); } } /* -------------------------------------------------------------------------- */ void SlaveElementInfoPerProc::synchronizeTags() { AKANTU_DEBUG_IN(); if (this->nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } /* --------<<<<-TAGS------------------------------------------------- */ - UInt mesh_data_sizes_buffer_length = 0; - CommunicationBuffer mesh_data_sizes_buffer; - - AKANTU_DEBUG_INFO( - "Receiving the size of the information about the mesh data tags."); - comm.broadcast(mesh_data_sizes_buffer_length, root); - - if (mesh_data_sizes_buffer_length != 0) { - mesh_data_sizes_buffer.resize(mesh_data_sizes_buffer_length); - AKANTU_DEBUG_INFO( - "Receiving the information about the mesh data tags, addr " - << (void *)mesh_data_sizes_buffer.storage()); - comm.broadcast(mesh_data_sizes_buffer, root); - AKANTU_DEBUG_INFO("Size of the information about the mesh data: " - << mesh_data_sizes_buffer_length); - - std::vector tag_names; - std::vector tag_type_codes; - std::vector tag_nb_component; - tag_names.resize(nb_tags); - tag_type_codes.resize(nb_tags); - tag_nb_component.resize(nb_tags); - CommunicationBuffer mesh_data_buffer; - UInt type_code_int; - for (UInt i(0); i < nb_tags; ++i) { - mesh_data_sizes_buffer >> tag_names[i]; - mesh_data_sizes_buffer >> type_code_int; - tag_type_codes[i] = static_cast(type_code_int); - mesh_data_sizes_buffer >> tag_nb_component[i]; - } - - std::vector::const_iterator names_it = tag_names.begin(); - std::vector::const_iterator names_end = tag_names.end(); - - CommunicationStatus mesh_data_comm_status; - AKANTU_DEBUG_INFO("Checking size of data to receive for mesh data TAG(" - << Tag::genTag(root, this->message_count, Tag::_MESH_DATA) - << ")"); - comm.probe(root, - Tag::genTag(root, this->message_count, Tag::_MESH_DATA), - mesh_data_comm_status); - UInt mesh_data_buffer_size(mesh_data_comm_status.size()); - AKANTU_DEBUG_INFO("Receiving " - << mesh_data_buffer_size << " bytes of mesh data TAG(" - << Tag::genTag(root, this->message_count, Tag::_MESH_DATA) - << ")"); - mesh_data_buffer.resize(mesh_data_buffer_size); - comm.receive(mesh_data_buffer, root, - Tag::genTag(root, this->message_count, Tag::_MESH_DATA)); - - // Loop over each tag for the current type - UInt k(0); - for (; names_it != names_end; ++names_it, ++k) { - this->fillMeshData(mesh_data_buffer, *names_it, tag_type_codes[k], - tag_nb_component[k]); - } + DynamicCommunicationBuffer mesh_data_sizes_buffer; + comm.broadcast(mesh_data_sizes_buffer, root); + AKANTU_DEBUG_INFO("Size of the information about the mesh data: " + << mesh_data_sizes_buffer.size()); + + if (mesh_data_sizes_buffer.size() == 0) + return; + + AKANTU_DEBUG_INFO("Receiving the information about the mesh data tags, addr " + << (void *)mesh_data_sizes_buffer.storage()); + + std::vector tag_names; + std::vector tag_type_codes; + std::vector tag_nb_component; + tag_names.resize(nb_tags); + tag_type_codes.resize(nb_tags); + tag_nb_component.resize(nb_tags); + CommunicationBuffer mesh_data_buffer; + UInt type_code_int; + for (UInt i(0); i < nb_tags; ++i) { + mesh_data_sizes_buffer >> tag_names[i]; + mesh_data_sizes_buffer >> type_code_int; + tag_type_codes[i] = static_cast(type_code_int); + mesh_data_sizes_buffer >> tag_nb_component[i]; + } + + std::vector::const_iterator names_it = tag_names.begin(); + std::vector::const_iterator names_end = tag_names.end(); + + CommunicationStatus mesh_data_comm_status; + AKANTU_DEBUG_INFO("Checking size of data to receive for mesh data TAG(" + << Tag::genTag(root, this->message_count, Tag::_MESH_DATA) + << ")"); + comm.probe(root, + Tag::genTag(root, this->message_count, Tag::_MESH_DATA), + mesh_data_comm_status); + UInt mesh_data_buffer_size(mesh_data_comm_status.size()); + AKANTU_DEBUG_INFO("Receiving " + << mesh_data_buffer_size << " bytes of mesh data TAG(" + << Tag::genTag(root, this->message_count, Tag::_MESH_DATA) + << ")"); + mesh_data_buffer.resize(mesh_data_buffer_size); + comm.receive(mesh_data_buffer, root, + Tag::genTag(root, this->message_count, Tag::_MESH_DATA)); + + // Loop over each tag for the current type + UInt k(0); + for (; names_it != names_end; ++names_it, ++k) { + this->fillMeshData(mesh_data_buffer, *names_it, tag_type_codes[k], + tag_nb_component[k]); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SlaveElementInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); UInt my_rank = comm.whoAmI(); AKANTU_DEBUG_INFO("Receiving element groups from proc " << root << " TAG(" << Tag::genTag(root, my_rank, Tag::_ELEMENT_GROUP) << ")"); CommunicationStatus status; comm.probe(root, Tag::genTag(root, my_rank, Tag::_ELEMENT_GROUP), status); CommunicationBuffer buffer(status.size()); comm.receive(buffer, root, Tag::genTag(root, my_rank, Tag::_ELEMENT_GROUP)); this->fillElementGroupsFromBuffer(buffer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -} // akantu +} // namespace akantu