diff --git a/src/common/aka_enum_macros.hh b/src/common/aka_enum_macros.hh index 5f38371fd..abfa27b5b 100644 --- a/src/common/aka_enum_macros.hh +++ b/src/common/aka_enum_macros.hh @@ -1,116 +1,133 @@ /** * @file aka_enum_macros.hh * * @author Nicolas Richart * * @date creation Wed Oct 31 2018 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_ENUM_MACROS_HH__ #define __AKANTU_AKA_ENUM_MACROS_HH__ #define AKANTU_PP_ENUM(s, data, i, elem) \ BOOST_PP_TUPLE_REM() \ elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) #if (defined(__GNUC__) || defined(__GNUG__)) #define AKA_GCC_VERSION \ (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) #if AKA_GCC_VERSION < 60000 #define AKANTU_ENUM_HASH(type_name) \ namespace std { \ template <> struct hash<::akantu::type_name> { \ using argument_type = ::akantu::type_name; \ size_t operator()(const argument_type & e) const noexcept { \ auto ue = underlying_type_t(e); \ return uh(ue); \ } \ \ private: \ const hash> uh{}; \ }; \ } #else #define AKANTU_ENUM_HASH(type_name) #endif // AKA_GCC_VERSION #endif // GNU #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ ({BOOST_PP_CAT(data, elem), BOOST_PP_STRINGIZE(elem)}) #define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data, elem)}) #define AKANTU_CLASS_ENUM_DECLARE(type_name, list) \ enum class type_name { \ BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ }; #define AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, prefix) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace std { \ inline string to_string(const ::akantu::type_name & type) { \ using namespace akantu; \ static unordered_map<::akantu::type_name, string> convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, prefix, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << std::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM_(type_name, list, prefix) \ inline std::istream & operator>>(std::istream & stream, type_name & type) { \ std::string str; \ stream >> str; \ static std::unordered_map convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, prefix, list))}; \ - type = convert.at(str); \ + try { \ + type = convert.at(str); \ + } catch (std::out_of_range &) { \ + std::ostringstream values; \ + std::for_each(convert.begin(), convert.end(), [&values](auto && pair) { \ + static bool first = true; \ + if (not first) \ + values << ", "; \ + values << "\"" << pair.first << "\""; \ + first = false; \ + }); \ + AKANTU_EXCEPTION("The value " << str << " is not a valid " \ + << BOOST_PP_STRINGIZE(type_name) \ + << " valid values are " << values.str()); \ + } \ return stream; \ } #define AKANTU_CLASS_ENUM_OUTPUT_STREAM(type_name, list) \ AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, type_name::_) #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, ) #define AKANTU_CLASS_ENUM_INPUT_STREAM(type_name, list) \ AKANTU_ENUM_INPUT_STREAM_(type_name, list, type_name::_) #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ AKANTU_ENUM_INPUT_STREAM_(type_name, list, ) #endif /* __AKANTU_AKA_ENUM_MACROS_HH__ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index b92532bbf..ed51811d6 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,1810 +1,1817 @@ /** * @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)); + + if(nb_local_element > 0) { + 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)); - + + if(nb_ghost_element > 0) { + 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); Element old_facet_el{type_facet, 0, gt_facet}; 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 = nb_facets + f_to_double.size() - 1; old_facet_el.element = f; /// update facet_to_element vector auto & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.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(el_facets.begin(), el_facets.end(), old_facet_el); AKANTU_DEBUG_ASSERT(f_update != el_facets.end(), "Facet not found"); f_update->element = new_facet; /// update elements connected to facet 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; } } } 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/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc index 9a1289a88..0961702e6 100644 --- a/src/solver/sparse_matrix_aij.cc +++ b/src/solver/sparse_matrix_aij.cc @@ -1,291 +1,291 @@ /** * @file sparse_matrix_aij.cc * * @author Nicolas Richart * * @date creation: Fri Aug 21 2015 * @date last modification: Mon Dec 04 2017 * * @brief Implementation of the AIJ sparse matrix * * @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 "sparse_matrix_aij.hh" #include "aka_iterators.hh" #include "dof_manager_default.hh" #include "dof_synchronizer.hh" #include "solver_vector_default.hh" #include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager), irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a") {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager), irn(matrix.irn, id + ":irn"), jcn(matrix.jcn, id + ":jcn"), a(matrix.a, id + ":a") {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::~SparseMatrixAIJ() = default; /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); auto begin = blocked_dofs.begin(); auto end = blocked_dofs.end(); auto is_blocked = [&](auto && i) -> bool { auto il = this->dof_manager.globalToLocalEquationNumber(i); return std::binary_search(begin, end, il); }; for (auto && ij_a : zip(irn, jcn, a)) { UInt ni = std::get<0>(ij_a) - 1; UInt nj = std::get<1>(ij_a) - 1; if (is_blocked(ni) or is_blocked(nj)) { // clang-format off std::get<2>(ij_a) = std::get<0>(ij_a) != std::get<1>(ij_a) ? 0. : this->dof_manager.isLocalOrMasterDOF(ni) ? block_val : 0.; // clang-format on } } this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); UInt m = this->size_; auto & comm = dof_manager.getCommunicator(); // write header if (comm.whoAmI() == 0) { outfile << "%%MatrixMarket matrix coordinate pattern"; if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << m << " " << m << " " << this->nb_non_zero << std::endl; } for (auto p : arange(comm.getNbProc())) { // write content if (comm.whoAmI() == p) { for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn.storage()[i] << " " << this->jcn.storage()[i] << " 1" << std::endl; } } comm.barrier(); } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); auto & comm = dof_manager.getCommunicator(); // open and set the properties of the stream std::ofstream outfile; if (0 == comm.whoAmI()) { outfile.open(filename.c_str()); } else { outfile.open(filename.c_str(), std::ios_base::app); } outfile.precision(std::numeric_limits::digits10); // write header decltype(nb_non_zero) nnz = this->nb_non_zero; comm.allReduce(nnz); if (comm.whoAmI() == 0) { outfile << "%%MatrixMarket matrix coordinate real"; if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << this->size_ << " " << this->size_ << " " << nnz << std::endl; } for (auto p : arange(comm.getNbProc())) { // write content if (comm.whoAmI() == p) { for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i) << std::endl; } } comm.barrier(); } // time to end outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const Array & x, Array & y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); y *= beta; auto i_it = this->irn.begin(); auto j_it = this->jcn.begin(); auto a_it = this->a.begin(); auto a_end = this->a.end(); auto x_it = x.begin_reinterpret(x.size() * x.getNbComponent()); auto y_it = y.begin_reinterpret(x.size() * x.getNbComponent()); for (; a_it != a_end; ++i_it, ++j_it, ++a_it) { Int i = this->dof_manager.globalToLocalEquationNumber(*i_it - 1); Int j = this->dof_manager.globalToLocalEquationNumber(*j_it - 1); const Real & A = *a_it; y_it[i] += alpha * A * x_it[j]; if ((this->matrix_type == _symmetric) && (i != j)) y_it[j] += alpha * A * x_it[i]; } if (this->dof_manager.hasSynchronizer()) this->dof_manager.getSynchronizer().reduceSynchronize(y); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const SolverVector & _x, SolverVector & _y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); auto && x = aka::as_type(_x).getVector(); auto && y = aka::as_type(_y).getVector(); this->matVecMul(x, y, alpha, beta); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) { AKANTU_DEBUG_IN(); const auto & mat = aka::as_type(matrix); AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(), "The to matrix don't have the same profiles"); memcpy(a.storage(), mat.getA().storage(), nb_non_zero * sizeof(Real)); this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyProfile(const SparseMatrix & other) { auto & A = aka::as_type(other); SparseMatrix::clearProfile(); - this->irn = A.irn; - this->jcn = A.jcn; + this->irn.copy(A.irn); + this->jcn.copy(A.jcn); this->irn_jcn_k.clear(); UInt i, j, k; for(auto && data : enumerate(irn, jcn)) { std::tie(k, i, j) = data; this->irn_jcn_k[this->key(i - 1, j - 1)] = k; } this->nb_non_zero = this->irn.size(); this->a.resize(this->nb_non_zero); this->a.set(0.); this->size_ = A.size_; - this->profile_release++; + this->profile_release = A.profile_release; this->value_release++; } /* -------------------------------------------------------------------------- */ template void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const { UInt i, j; Real A_ij; for (auto && tuple : zip(irn, jcn, a)) { std::tie(i, j, A_ij) = tuple; B.add(i - 1, j - 1, alpha * A_ij); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const { if (aka::is_of_type(B)) { this->addMeToTemplated(aka::as_type(B), alpha); } else { // this->addMeToTemplated(*this, alpha); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::mul(Real alpha) { this->a *= alpha; this->value_release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::clear() { a.set(0.); this->value_release++; } } // namespace akantu diff --git a/src/synchronizer/node_info_per_processor.cc b/src/synchronizer/node_info_per_processor.cc index ddae37c59..f81bcae63 100644 --- a/src/synchronizer/node_info_per_processor.cc +++ b/src/synchronizer/node_info_per_processor.cc @@ -1,851 +1,849 @@ /** * @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; - + for (const auto & send_info : make_view(master_info, 2)) { 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); + std::cout << "Proc " << rank << " sends " << sends.size() + << " nodes to proc " << send_schemes.first << std::endl; } } { // 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); + std::cout << "Proc " << rank << " receives " << recvs.size() + << " nodes to proc " << recv_schemes.first << std::endl; + } } 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)) { + 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 MeshDataTypeCode::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, + for (const auto & node : nodespp) { + Vector coord(nodes.storage() + spatial_dimension * node, 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 MeshDataTypeCode::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); } AKANTU_DEBUG_INFO( "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)); 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() { 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/test/test_gtest_utils.hh b/test/test_gtest_utils.hh index 98551f437..f136abf9b 100644 --- a/test/test_gtest_utils.hh +++ b/test/test_gtest_utils.hh @@ -1,232 +1,242 @@ /** * @file test_gtest_utils.hh * * @author Nicolas Richart * * @date creation: Tue Nov 14 2017 * @date last modification: Wed Feb 21 2018 * * @brief Utils to help write tests * * @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_common.hh" #include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_GTEST_UTILS_HH__ #define __AKANTU_TEST_GTEST_UTILS_HH__ namespace { /* -------------------------------------------------------------------------- */ template <::akantu::ElementType t> using element_type_t = std::integral_constant<::akantu::ElementType, t>; /* -------------------------------------------------------------------------- */ template struct gtest_list {}; template struct gtest_list> { using type = ::testing::Types; }; template using gtest_list_t = typename gtest_list::type; /* -------------------------------------------------------------------------- */ template struct tuple_concat {}; template struct tuple_concat, std::tuple> { using type = std::tuple; }; template using tuple_concat_t = typename tuple_concat::type; /* -------------------------------------------------------------------------- */ template