diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc index 58664e486..b9dc383a0 100644 --- a/src/mesh/element_group.cc +++ b/src/mesh/element_group.cc @@ -1,204 +1,213 @@ /** * @file element_group.cc * * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Mon Jan 22 2018 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_csr.hh" #include "dumpable.hh" #include "dumpable_inline_impl.hh" #include "group_manager.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_utils.hh" #include #include #include #include "element_group.hh" #if defined(AKANTU_USE_IOHELPER) #include "dumper_iohelper_paraview.hh" #endif namespace akantu { /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const std::string & group_name, const Mesh & mesh, NodeGroup & node_group, UInt dimension, const std::string & id, const MemoryID & mem_id) : Memory(id, mem_id), mesh(mesh), name(group_name), elements("elements", id, mem_id), node_group(node_group), dimension(dimension) { AKANTU_DEBUG_IN(); #if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_" + group_name, group_name, true); this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), _all_dimensions); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const ElementGroup & /*other*/) = default; /* -------------------------------------------------------------------------- */ -void ElementGroup::clear() { elements.free(); } +void ElementGroup::clear() { + elements.free(); +} + +/* -------------------------------------------------------------------------- */ +void ElementGroup::clear(ElementType type, GhostType ghost_type) { + if (elements.exists(type, ghost_type)) { + elements(type, ghost_type).clear(); + } +} /* -------------------------------------------------------------------------- */ bool ElementGroup::empty() const { return elements.empty(); } /* -------------------------------------------------------------------------- */ void ElementGroup::append(const ElementGroup & other_group) { AKANTU_DEBUG_IN(); node_group.append(other_group.node_group); /// loop on all element types in all dimensions for (auto ghost_type : ghost_types) { for (auto type : other_group.elementTypes(_ghost_type = ghost_type, _element_kind = _ek_not_defined)) { const Array & other_elem_list = other_group.elements(type, ghost_type); UInt nb_other_elem = other_elem_list.size(); Array * elem_list; UInt nb_elem = 0; /// create current type if doesn't exists, otherwise get information if (elements.exists(type, ghost_type)) { elem_list = &elements(type, ghost_type); nb_elem = elem_list->size(); } else { elem_list = &(elements.alloc(0, 1, type, ghost_type)); } /// append new elements to current list elem_list->resize(nb_elem + nb_other_elem); std::copy(other_elem_list.begin(), other_elem_list.end(), elem_list->begin() + nb_elem); /// remove duplicates std::sort(elem_list->begin(), elem_list->end()); Array::iterator<> end = std::unique(elem_list->begin(), elem_list->end()); elem_list->resize(end - elem_list->begin()); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) { ; } stream << space << "ElementGroup [" << std::endl; stream << space << " + name: " << name << std::endl; stream << space << " + dimension: " << dimension << std::endl; elements.printself(stream, indent + 1); node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void ElementGroup::optimize() { // increasing the locality of data when iterating on the element of a group for (auto ghost_type : ghost_types) { for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { Array & els = elements(type, ghost_type); std::sort(els.begin(), els.end()); Array::iterator<> end = std::unique(els.begin(), els.end()); els.resize(end - els.begin()); } } node_group.optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::fillFromNodeGroup() { CSR node_to_elem; MeshUtils::buildNode2Elements(this->mesh, node_to_elem, this->dimension); std::set seen; Array::const_iterator<> itn = this->node_group.begin(); Array::const_iterator<> endn = this->node_group.end(); for (; itn != endn; ++itn) { CSR::iterator ite = node_to_elem.begin(*itn); CSR::iterator ende = node_to_elem.end(*itn); for (; ite != ende; ++ite) { const Element & elem = *ite; if (this->dimension != _all_dimensions && this->dimension != Mesh::getSpatialDimension(elem.type)) { continue; } if (seen.find(elem) != seen.end()) { continue; } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(elem.type); Array::const_iterator> conn_it = this->mesh.getConnectivity(elem.type, elem.ghost_type) .begin(nb_nodes_per_element); const Vector & conn = conn_it[elem.element]; UInt count = 0; for (UInt n = 0; n < conn.size(); ++n) { count += (this->node_group.getNodes().find(conn(n)) != UInt(-1) ? 1 : 0); } if (count == nb_nodes_per_element) { this->add(elem); } seen.insert(elem); } } this->optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::addDimension(UInt dimension) { this->dimension = std::max(dimension, this->dimension); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 1eb483f1c..b1a8eaa48 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,201 +1,202 @@ /** * @file element_group.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Wed Nov 08 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element_type_map.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_GROUP_HH_ #define AKANTU_ELEMENT_GROUP_HH_ namespace akantu { class Mesh; class Element; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ElementGroup : private Memory, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ElementGroup(const std::string & name, const Mesh & mesh, NodeGroup & node_group, UInt dimension = _all_dimensions, const std::string & id = "element_group", const MemoryID & mem_id = 0); ElementGroup(const ElementGroup & /*unused*/); /* ------------------------------------------------------------------------ */ /* Type definitions */ /* ------------------------------------------------------------------------ */ public: using ElementList = ElementTypeMapArray; using NodeList = Array; /* ------------------------------------------------------------------------ */ /* Element iterator */ /* ------------------------------------------------------------------------ */ using type_iterator = ElementList::type_iterator; [[deprecated("Use elementTypes instead")]] inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; template inline decltype(auto) elementTypes(pack &&... _pack) const { return elements.elementTypes(_pack...); } using const_element_iterator = Array::const_iterator; inline const_element_iterator begin(ElementType type, GhostType ghost_type = _not_ghost) const; inline const_element_iterator end(ElementType type, GhostType ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// empty the element group void clear(); + void clear(ElementType type, GhostType ghost_type = _not_ghost); bool empty() const __attribute__((warn_unused_result)); /// append another group to this group /// BE CAREFUL: it doesn't conserve the element order void append(const ElementGroup & other_group); /// add an element to the group. By default the it does not add the nodes to /// the group inline void add(const Element & el, bool add_nodes = false, bool check_for_duplicate = true); /// \todo fix the default for add_nodes : make it coherent with the other /// method inline void add(ElementType type, UInt element, GhostType ghost_type = _not_ghost, bool add_nodes = true, bool check_for_duplicate = true); inline void addNode(UInt node_id, bool check_for_duplicate = true); inline void removeNode(UInt node_id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the elements based on the underlying node group. virtual void fillFromNodeGroup(); // sort and remove duplicated values void optimize(); /// change the dimension if needed void addDimension(UInt dimension); private: inline void addElement(ElementType elem_type, UInt elem_id, GhostType ghost_type); friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const Array & getElements(ElementType type, GhostType ghost_type = _not_ghost) const; AKANTU_GET_MACRO(Elements, elements, const ElementTypeMapArray &); AKANTU_GET_MACRO_NOT_CONST(Elements, elements, ElementTypeMapArray &); template auto size(Args &&... pack) const { return elements.size(std::forward(pack)...); } // AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); AKANTU_GET_MACRO(NodeGroup, node_group, const NodeGroup &); AKANTU_GET_MACRO_NOT_CONST(NodeGroup, node_group, NodeGroup &); AKANTU_GET_MACRO(Dimension, dimension, UInt); AKANTU_GET_MACRO(Name, name, std::string); inline UInt getNbNodes() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Mesh to which this group belongs const Mesh & mesh; /// name of the group std::string name; /// list of elements composing the group ElementList elements; /// sub list of nodes which are composing the elements NodeGroup & node_group; /// group dimension UInt dimension{_all_dimensions}; /// empty arry for the iterator to work when an element type not present Array empty_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ElementGroup & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "element.hh" #include "element_group_inline_impl.hh" #endif /* AKANTU_ELEMENT_GROUP_HH_ */ diff --git a/src/mesh/group_manager.cc b/src/mesh/group_manager.cc index 105b1f491..70dba4f1b 100644 --- a/src/mesh/group_manager.cc +++ b/src/mesh/group_manager.cc @@ -1,1036 +1,1033 @@ /** * @file group_manager.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Feb 20 2018 * * @brief Stores information about ElementGroup and NodeGroup * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "group_manager.hh" #include "aka_csr.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ GroupManager::GroupManager(Mesh & mesh, const ID & id, const MemoryID & mem_id) : id(id), memory_id(mem_id), mesh(mesh) { AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ GroupManager::~GroupManager() = default; /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::createNodeGroup(const std::string & group_name, bool replace_group) { AKANTU_DEBUG_IN(); auto it = node_groups.find(group_name); if (it != node_groups.end()) { if (replace_group) { it->second.reset(); } else { AKANTU_EXCEPTION( "Trying to create a node group that already exists:" << group_name); } } std::stringstream sstr; sstr << this->id << ":" << group_name << "_node_group"; auto && ptr = std::make_unique(group_name, mesh, sstr.str(), memory_id); auto & node_group = *ptr; // \todo insert_or_assign in c++17 if (it != node_groups.end()) { it->second = std::move(ptr); } else { node_groups[group_name] = std::move(ptr); } AKANTU_DEBUG_OUT(); return node_group; } /* -------------------------------------------------------------------------- */ template NodeGroup & GroupManager::createFilteredNodeGroup(const std::string & group_name, const NodeGroup & source_node_group, T & filter) { AKANTU_DEBUG_IN(); NodeGroup & node_group = this->createNodeGroup(group_name); node_group.append(source_node_group); if (T::type == FilterFunctor::_node_filter_functor) { node_group.applyNodeFilter(filter); } else { AKANTU_ERROR("ElementFilter cannot be applied to NodeGroup yet." << " Needs to be implemented."); } AKANTU_DEBUG_OUT(); return node_group; } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, bool replace_group) { AKANTU_DEBUG_IN(); auto it = element_groups.find(group_name); if (it != element_groups.end()) { if (replace_group) { it->second.reset(); } else { AKANTU_EXCEPTION("Trying to create a element group that already exists:" << group_name); } } NodeGroup & new_node_group = createNodeGroup(group_name + "_nodes", replace_group); auto && ptr = std::make_unique( group_name, mesh, new_node_group, dimension, this->id + ":" + group_name + "_element_group", memory_id); auto & element_group = *ptr; if (it != element_groups.end()) { it->second = std::move(ptr); } else { element_groups[group_name] = std::move(ptr); } AKANTU_DEBUG_OUT(); return element_group; } /* -------------------------------------------------------------------------- */ void GroupManager::destroyElementGroup(const std::string & group_name, bool destroy_node_group) { AKANTU_DEBUG_IN(); auto eit = element_groups.find(group_name); if (eit != element_groups.end()) { if (destroy_node_group) { destroyNodeGroup(eit->second->getNodeGroup().getName()); } element_groups.erase(eit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::destroyNodeGroup(const std::string & group_name) { AKANTU_DEBUG_IN(); auto nit = node_groups.find(group_name); if (nit != node_groups.end()) { node_groups.erase(nit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group) { AKANTU_DEBUG_IN(); if (element_groups.find(group_name) != element_groups.end()) { AKANTU_EXCEPTION( "Trying to create a element group that already exists:" << group_name); } auto && ptr = std::make_unique( group_name, mesh, node_group, dimension, id + ":" + group_name + "_element_group", memory_id); auto & element_group = *ptr; element_groups[group_name] = std::move(ptr); AKANTU_DEBUG_OUT(); return element_group; } /* -------------------------------------------------------------------------- */ template ElementGroup & GroupManager::createFilteredElementGroup( const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter) { AKANTU_DEBUG_IN(); if (T::type == FilterFunctor::_node_filter_functor) { auto & filtered_node_group = this->createFilteredNodeGroup( group_name + "_nodes", node_group, filter); AKANTU_DEBUG_OUT(); return this->createElementGroup(group_name, dimension, filtered_node_group); } if (T::type == FilterFunctor::_element_filter_functor) { AKANTU_ERROR( "Cannot handle an ElementFilter yet. Needs to be implemented."); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ class ClusterSynchronizer : public DataAccessor { using DistantIDs = std::set>; public: ClusterSynchronizer(GroupManager & group_manager, UInt element_dimension, std::string cluster_name_prefix, ElementTypeMapArray & element_to_fragment, const ElementSynchronizer & element_synchronizer, UInt nb_cluster) : group_manager(group_manager), element_dimension(element_dimension), cluster_name_prefix(std::move(cluster_name_prefix)), element_to_fragment(element_to_fragment), element_synchronizer(element_synchronizer), nb_cluster(nb_cluster) {} UInt synchronize() { Communicator & comm = Communicator::getStaticCommunicator(); UInt rank = comm.whoAmI(); UInt nb_proc = comm.getNbProc(); /// find starting index to renumber local clusters Array nb_cluster_per_proc(nb_proc); nb_cluster_per_proc(rank) = nb_cluster; comm.allGather(nb_cluster_per_proc); starting_index = std::accumulate(nb_cluster_per_proc.begin(), nb_cluster_per_proc.begin() + rank, 0U); UInt global_nb_fragment = std::accumulate(nb_cluster_per_proc.begin() + rank, nb_cluster_per_proc.end(), starting_index); /// create the local to distant cluster pairs with neighbors element_synchronizer.synchronizeOnce(*this, SynchronizationTag::_gm_clusters); /// count total number of pairs Array nb_pairs(nb_proc); // This is potentially a bug for more than // 2**31 pairs, but due to a all gatherv after // it must be int to match MPI interfaces nb_pairs(rank) = distant_ids.size(); comm.allGather(nb_pairs); UInt total_nb_pairs = std::accumulate(nb_pairs.begin(), nb_pairs.end(), 0); /// generate pairs global array UInt local_pair_index = std::accumulate(nb_pairs.storage(), nb_pairs.storage() + rank, 0); Array total_pairs(total_nb_pairs, 2); for (const auto & ids : distant_ids) { total_pairs(local_pair_index, 0) = ids.first; total_pairs(local_pair_index, 1) = ids.second; ++local_pair_index; } /// communicate pairs to all processors nb_pairs *= 2; comm.allGatherV(total_pairs, nb_pairs); /// renumber clusters /// generate fragment list std::vector> global_clusters; UInt total_nb_cluster = 0; Array is_fragment_in_cluster(global_nb_fragment, 1, false); std::queue fragment_check_list; while (not total_pairs.empty()) { /// create a new cluster ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; UInt first_fragment = total_pairs(0, 0); UInt second_fragment = total_pairs(0, 1); total_pairs.erase(0); fragment_check_list.push(first_fragment); fragment_check_list.push(second_fragment); while (!fragment_check_list.empty()) { UInt current_fragment = fragment_check_list.front(); UInt * total_pairs_end = total_pairs.storage() + total_pairs.size() * 2; UInt * fragment_found = std::find(total_pairs.storage(), total_pairs_end, current_fragment); if (fragment_found != total_pairs_end) { UInt position = fragment_found - total_pairs.storage(); UInt pair = position / 2; UInt other_index = (position + 1) % 2; fragment_check_list.push(total_pairs(pair, other_index)); total_pairs.erase(pair); } else { fragment_check_list.pop(); current_cluster.insert(current_fragment); is_fragment_in_cluster(current_fragment) = true; } } } /// add to FragmentToCluster all local fragments for (UInt c = 0; c < global_nb_fragment; ++c) { if (!is_fragment_in_cluster(c)) { ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; current_cluster.insert(c); } } /// reorganize element groups to match global clusters for (UInt c = 0; c < global_clusters.size(); ++c) { /// create new element group corresponding to current cluster std::stringstream sstr; sstr << cluster_name_prefix << "_" << c; ElementGroup & cluster = group_manager.createElementGroup(sstr.str(), element_dimension, true); auto it = global_clusters[c].begin(); auto end = global_clusters[c].end(); /// append to current element group all fragments that belong to /// the same cluster if they exist for (; it != end; ++it) { Int local_index = *it - starting_index; if (local_index < 0 || local_index >= Int(nb_cluster)) { continue; } std::stringstream tmp_sstr; tmp_sstr << "tmp_" << cluster_name_prefix << "_" << local_index; AKANTU_DEBUG_ASSERT(group_manager.elementGroupExists(tmp_sstr.str()), "Temporary fragment \"" << tmp_sstr.str() << "\" not found"); cluster.append(group_manager.getElementGroup(tmp_sstr.str())); group_manager.destroyElementGroup(tmp_sstr.str(), true); } } return total_nb_cluster; } private: /// functions for parallel communications inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_gm_clusters) { return elements.size() * sizeof(UInt); } return 0; } inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag != SynchronizationTag::_gm_clusters) { return; } Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; /// for each element pack its global cluster index buffer << element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; } } inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag != SynchronizationTag::_gm_clusters) { return; } Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { UInt distant_cluster; buffer >> distant_cluster; const Element & el = *el_it; UInt local_cluster = element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; distant_ids.insert(std::make_pair(local_cluster, distant_cluster)); } } private: GroupManager & group_manager; UInt element_dimension; std::string cluster_name_prefix; ElementTypeMapArray & element_to_fragment; const ElementSynchronizer & element_synchronizer; UInt nb_cluster; DistantIDs distant_ids; UInt starting_index; }; /* -------------------------------------------------------------------------- */ /// \todo this function doesn't work in 1D UInt GroupManager::createBoundaryGroupFromGeometry() { UInt spatial_dimension = mesh.getSpatialDimension(); return createClusters(spatial_dimension - 1, "boundary"); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, Mesh & mesh_facets, const std::string & cluster_name_prefix, const GroupManager::ClusteringFilter & filter) { return createClusters(element_dimension, cluster_name_prefix, filter, mesh_facets); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, const std::string & cluster_name_prefix, const GroupManager::ClusteringFilter & filter) { MeshAccessor mesh_accessor(mesh); auto mesh_facets = std::make_unique(mesh.getSpatialDimension(), mesh_accessor.getNodesSharedPtr(), "mesh_facets_for_clusters"); mesh_facets->defineMeshParent(mesh); MeshUtils::buildAllFacets(mesh, *mesh_facets, element_dimension, element_dimension - 1); return createClusters(element_dimension, cluster_name_prefix, filter, *mesh_facets); } /* -------------------------------------------------------------------------- */ //// \todo if needed element list construction can be optimized by //// templating the filter class UInt GroupManager::createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const GroupManager::ClusteringFilter & filter, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt nb_proc = mesh.getCommunicator().getNbProc(); std::string tmp_cluster_name_prefix = cluster_name_prefix; ElementTypeMapArray * element_to_fragment = nullptr; if (nb_proc > 1 && mesh.isDistributed()) { element_to_fragment = new ElementTypeMapArray("element_to_fragment", id, memory_id); element_to_fragment->initialize( mesh, _nb_component = 1, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); // mesh.initElementTypeMapArray(*element_to_fragment, 1, element_dimension, // false, _ek_not_defined, true); tmp_cluster_name_prefix = "tmp_" + tmp_cluster_name_prefix; } ElementTypeMapArray seen_elements("seen_elements", id, memory_id); seen_elements.initialize(mesh, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); // mesh.initElementTypeMapArray(seen_elements, 1, element_dimension, false, // _ek_not_defined, true); for (auto ghost_type : ghost_types) { Element el; el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { el.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & seen_elements_array = seen_elements(type, ghost_type); for (UInt e = 0; e < nb_element; ++e) { el.element = e; if (!filter(el)) { seen_elements_array(e) = true; } } } } Array checked_node(mesh.getNbNodes(), 1, false); UInt nb_cluster = 0; for (auto ghost_type : ghost_types) { Element uns_el; uns_el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { uns_el.type = type; Array & seen_elements_vec = seen_elements(uns_el.type, uns_el.ghost_type); for (UInt e = 0; e < seen_elements_vec.size(); ++e) { // skip elements that have been already seen if (seen_elements_vec(e)) { continue; } // set current element uns_el.element = e; seen_elements_vec(e) = true; /// create a new cluster std::stringstream sstr; sstr << tmp_cluster_name_prefix << "_" << nb_cluster; ElementGroup & cluster = createElementGroup(sstr.str(), element_dimension, true); ++nb_cluster; // point element are cluster by themself if (element_dimension == 0) { cluster.add(uns_el); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(uns_el.type); Vector connect = mesh.getConnectivity(uns_el.type, uns_el.ghost_type) .begin(nb_nodes_per_element)[uns_el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { /// add element's nodes to the cluster UInt node = connect[n]; if (!checked_node(node)) { cluster.addNode(node); checked_node(node) = true; } } continue; } std::queue element_to_add; element_to_add.push(uns_el); /// keep looping until current cluster is complete (no more /// connected elements) while (!element_to_add.empty()) { /// take first element and erase it in the queue Element el = element_to_add.front(); element_to_add.pop(); /// if parallel, store cluster index per element if (nb_proc > 1 && mesh.isDistributed()) { (*element_to_fragment)(el.type, el.ghost_type)(el.element) = nb_cluster - 1; } /// add current element to the cluster cluster.add(el); const Array & element_to_facet = mesh_facets.getSubelementToElement(el.type, el.ghost_type); UInt nb_facet_per_element = element_to_facet.getNbComponent(); for (UInt f = 0; f < nb_facet_per_element; ++f) { const Element & facet = element_to_facet(el.element, f); if (facet == ElementNull) { continue; } const std::vector & connected_elements = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); for (UInt elem = 0; elem < connected_elements.size(); ++elem) { const Element & check_el = connected_elements[elem]; // check if this element has to be skipped if (check_el == ElementNull || check_el == el) { continue; } Array & seen_elements_vec_current = seen_elements(check_el.type, check_el.ghost_type); if (not seen_elements_vec_current(check_el.element)) { seen_elements_vec_current(check_el.element) = true; element_to_add.push(check_el); } } } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = mesh.getConnectivity(el.type, el.ghost_type) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { /// add element's nodes to the cluster UInt node = connect[n]; if (!checked_node(node)) { cluster.addNode(node, false); checked_node(node) = true; } } } } } } if (nb_proc > 1 && mesh.isDistributed()) { ClusterSynchronizer cluster_synchronizer( *this, element_dimension, cluster_name_prefix, *element_to_fragment, this->mesh.getElementSynchronizer(), nb_cluster); nb_cluster = cluster_synchronizer.synchronize(); delete element_to_fragment; } if (mesh.isDistributed()) { this->synchronizeGroupNames(); } AKANTU_DEBUG_OUT(); return nb_cluster; } /* -------------------------------------------------------------------------- */ template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name) { std::set group_names; const auto & datas = mesh.getData(dataset_name); std::map group_dim; for (auto ghost_type : ghost_types) { for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements (" << type << ":" << ghost_type << ") in the map from MeshData (" << dataset.size() << ") " << dataset_name << " and in the mesh (" << nb_element << ")!"); for (UInt e(0); e < nb_element; ++e) { std::stringstream sstr; sstr << dataset(e); std::string gname = sstr.str(); group_names.insert(gname); auto it = group_dim.find(gname); if (it == group_dim.end()) { group_dim[gname] = mesh.getSpatialDimension(type); } else { it->second = std::max(it->second, mesh.getSpatialDimension(type)); } } } } - auto git = group_names.begin(); - auto gend = group_names.end(); - for (; git != gend; ++git) { - createElementGroup(*git, group_dim[*git]); + for (auto && name : group_names) { + createElementGroup(name, group_dim[name]); } if (mesh.isDistributed()) { this->synchronizeGroupNames(); } Element el; for (auto ghost_type : ghost_types) { el.ghost_type = ghost_type; for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { el.type = type; const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements in the map from " "MeshData and in the mesh!"); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(el.type); Array::const_iterator> cit = mesh.getConnectivity(type, ghost_type).begin(nb_nodes_per_element); for (UInt e(0); e < nb_element; ++e, ++cit) { el.element = e; std::stringstream sstr; sstr << dataset(e); ElementGroup & group = getElementGroup(sstr.str()); group.add(el, false, false); const Vector & connect = *cit; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; group.addNode(node, false); } } } } - git = group_names.begin(); - for (; git != gend; ++git) { - getElementGroup(*git).optimize(); + for (auto && name : group_names) { + getElementGroup(name).optimize(); } } template void GroupManager::createGroupsFromMeshData( const std::string & dataset_name); template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name); /* -------------------------------------------------------------------------- */ void GroupManager::createElementGroupFromNodeGroup( const std::string & name, const std::string & node_group_name, UInt dimension) { NodeGroup & node_group = getNodeGroup(node_group_name); ElementGroup & group = createElementGroup(name, dimension, node_group); group.fillFromNodeGroup(); } /* -------------------------------------------------------------------------- */ void GroupManager::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "GroupManager [" << std::endl; std::set node_group_seen; for (auto & group : iterateElementGroups()) { group.printself(stream, indent + 1); node_group_seen.insert(group.getNodeGroup().getName()); } for (auto & group : iterateNodeGroups()) { if (node_group_seen.find(group.getName()) == node_group_seen.end()) { group.printself(stream, indent + 1); } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ UInt GroupManager::getNbElementGroups(UInt dimension) const { if (dimension == _all_dimensions) { return element_groups.size(); } return std::count_if( element_groups.begin(), element_groups.end(), [dimension](auto && eg) { return eg.second->getDimension() == dimension; }); } /* -------------------------------------------------------------------------- */ void GroupManager::checkAndAddGroups(DynamicCommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); UInt nb_node_group; buffer >> nb_node_group; AKANTU_DEBUG_INFO("Received " << nb_node_group << " node group names"); for (UInt ng = 0; ng < nb_node_group; ++ng) { std::string node_group_name; buffer >> node_group_name; if (node_groups.find(node_group_name) == node_groups.end()) { this->createNodeGroup(node_group_name); } AKANTU_DEBUG_INFO("Received node goup name: " << node_group_name); } UInt nb_element_group; buffer >> nb_element_group; AKANTU_DEBUG_INFO("Received " << nb_element_group << " element group names"); for (UInt eg = 0; eg < nb_element_group; ++eg) { std::string element_group_name; buffer >> element_group_name; std::string node_group_name; buffer >> node_group_name; UInt dim; buffer >> dim; AKANTU_DEBUG_INFO("Received element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with node group " << node_group_name); NodeGroup & node_group = *node_groups[node_group_name]; if (element_groups.find(element_group_name) == element_groups.end()) { this->createElementGroup(element_group_name, dim, node_group); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::fillBufferWithGroupNames( DynamicCommunicationBuffer & comm_buffer) const { AKANTU_DEBUG_IN(); // packing node group names; UInt nb_groups = this->node_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " node group names"); auto nnames_it = node_groups.begin(); auto nnames_end = node_groups.end(); for (; nnames_it != nnames_end; ++nnames_it) { std::string node_group_name = nnames_it->first; comm_buffer << node_group_name; AKANTU_DEBUG_INFO("Sending node goupe name: " << node_group_name); } // packing element group names with there associated node group name nb_groups = this->element_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " element group names"); auto gnames_it = this->element_groups.begin(); auto gnames_end = this->element_groups.end(); for (; gnames_it != gnames_end; ++gnames_it) { ElementGroup & element_group = *(gnames_it->second); std::string element_group_name = gnames_it->first; std::string node_group_name = element_group.getNodeGroup().getName(); UInt dim = element_group.getDimension(); comm_buffer << element_group_name; comm_buffer << node_group_name; comm_buffer << dim; AKANTU_DEBUG_INFO("Sending element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with the node group " << node_group_name); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::synchronizeGroupNames() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); Int nb_proc = comm.getNbProc(); Int my_rank = comm.whoAmI(); if (nb_proc == 1) { return; } if (my_rank == 0) { for (Int p = 1; p < nb_proc; ++p) { DynamicCommunicationBuffer recv_buffer; auto tag = Tag::genTag(p, 0, Tag::_element_group); comm.receive(recv_buffer, p, tag); AKANTU_DEBUG_INFO("Got " << printMemorySize(recv_buffer.size()) << " from proc " << p << " " << tag); this->checkAndAddGroups(recv_buffer); } DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); AKANTU_DEBUG_INFO("Initiating broadcast with " << printMemorySize(comm_buffer.size())); comm.broadcast(comm_buffer); } else { DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); auto tag = Tag::genTag(my_rank, 0, Tag::_element_group); AKANTU_DEBUG_INFO("Sending " << printMemorySize(comm_buffer.size()) << " to proc " << 0 << " " << tag); comm.send(comm_buffer, 0, tag); DynamicCommunicationBuffer recv_buffer; comm.broadcast(recv_buffer); AKANTU_DEBUG_INFO("Receiving broadcast with " << printMemorySize(recv_buffer.size())); this->checkAndAddGroups(recv_buffer); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const ElementGroup & GroupManager::getElementGroup(const std::string & name) const { auto it = element_groups.find(name); if (it == element_groups.end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::getElementGroup(const std::string & name) { auto it = element_groups.find(name); if (it == element_groups.end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const NodeGroup & GroupManager::getNodeGroup(const std::string & name) const { auto it = node_groups.find(name); if (it == node_groups.end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::getNodeGroup(const std::string & name) { auto it = node_groups.find(name); if (it == node_groups.end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ template void GroupManager::renameGroup(GroupsType & groups, const std::string & name, const std::string & new_name) { auto it = groups.find(name); if (it == groups.end()) { AKANTU_EXCEPTION("There are no group named " << name << " associated to the group manager: " << id); } auto && group_ptr = std::move(it->second); group_ptr->name = new_name; groups.erase(it); groups[new_name] = std::move(group_ptr); } /* -------------------------------------------------------------------------- */ void GroupManager::renameElementGroup(const std::string & name, const std::string & new_name) { renameGroup(element_groups, name, new_name); } /* -------------------------------------------------------------------------- */ void GroupManager::renameNodeGroup(const std::string & name, const std::string & new_name) { renameGroup(node_groups, name, new_name); } /* -------------------------------------------------------------------------- */ void GroupManager::copyElementGroup(const std::string & name, const std::string & new_name) { const auto & grp = getElementGroup(name); auto & new_grp = createElementGroup(new_name, grp.getDimension()); new_grp.getElements().copy(grp.getElements()); } /* -------------------------------------------------------------------------- */ void GroupManager::copyNodeGroup(const std::string & name, const std::string & new_name) { const auto & grp = getNodeGroup(name); auto & new_grp = createNodeGroup(new_name); new_grp.getNodes().copy(grp.getNodes()); } } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 4461283aa..3405d129a 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,464 +1,464 @@ /** * @file structural_mechanics_model.cc * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Wed Feb 21 2018 * * @brief Model implementation for Structural Mechanics elements * * * 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 "structural_mechanics_model.hh" #include "dof_manager.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "shape_structural.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_iohelper_paraview.hh" #include "group_manager_inline_impl.hh" #endif /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model_inline_impl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), time_step(NAN), f_m2a(1.0), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) { nb_degree_of_freedom = 3; } else if (spatial_dimension == 3) { nb_degree_of_freedom = 6; } else { AKANTU_TO_IMPLEMENT(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("structural_mechanics_model", id, true); #endif this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { // <<<< This is the SolidMechanicsModel implementation for future ref >>>> // material_index.initialize(mesh, _element_kind = _ek_not_defined, // _default_value = UInt(-1), _with_nb_element = // true); // material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, // _with_nb_element = true); // Model::initFullImpl(options); // // initialize pbc // if (this->pbc_pair.size() != 0) // this->initPBC(); // // initialize the materials // if (this->parser.getLastParsedFile() != "") { // this->instantiateMaterials(); // } // this->initMaterials(); // this->initBC(*this, *displacement, *displacement_increment, // *external_force); // <<<< END >>>> Model::initFullImpl(options); // Initializing stresses ElementTypeMap stress_components; /// TODO this is ugly af, maybe add a function to FEEngine for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { UInt nb_components = 0; // Getting number of components for each element type #define GET_(type) nb_components = ElementClass::getNbStressComponents() AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_); #undef GET_ stress_components(nb_components, type); } stress.initialize( getFEEngine(), _spatial_dimension = _all_dimensions, - _element_kind = _ek_structural, _all_ghost_types = true, + _element_kind = _ek_structural, _nb_component = [&stress_components](ElementType type, GhostType /*unused*/) -> UInt { return stress_components(type); }); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFEEngineBoundary() { /// TODO: this function should not be reimplemented /// we're just avoiding a call to Model::initFEEngineBoundary() } /* -------------------------------------------------------------------------- */ // void StructuralMechanicsModel::setTimeStep(Real time_step) { // this->time_step = time_step; // #if defined(AKANTU_USE_IOHELPER) // this->mesh.getDumper().setTimeStep(time_step); // #endif // } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(velocity, spatial_dimension, "velocity"); this->allocNodalField(acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { for (auto && type : mesh.elementTypes(_element_kind = _ek_structural)) { // computeRotationMatrix(type); element_material.alloc(mesh.getNbElement(type), 1, type); } getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); getDOFManager().getMatrix("K").zero(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(ElementType type) { Mesh & mesh = getFEEngine().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if (!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom * nb_nodes_per_element * nb_degree_of_freedom * nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).zero(); Array rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.zero(); #define COMPUTE_ROTATION_MATRIX(type) computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); auto T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { auto & T = *T_it; auto & R = *R_it; for (UInt k = 0; k < nb_nodes_per_element; ++k) { for (UInt i = 0; i < nb_degree_of_freedom; ++i) { for (UInt j = 0; j < nb_degree_of_freedom; ++j) { T(k * nb_degree_of_freedom + i, k * nb_degree_of_freedom + j) = R(i, j); } } } } } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; return mesh.createNodalField(uint_nodal_fields[field_name], group_name); } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } UInt padding_size = 0; if (padding_flag) { padding_size = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0, padding_size); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation, group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "force") { return mesh.createStridedNodalField(external_force, group_name, n, 0, padding_size); } if (field_name == "momentum") { return mesh.createStridedNodalField( external_force, group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force, group_name, n, 0, padding_size); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField( internal_force, group_name, nb_degree_of_freedom - n, n, padding_size); } return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*unused*/, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; if (field_name == "element_index_by_material") { field = mesh.createElementalField( field_name, group_name, spatial_dimension, kind); } return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") { assembleStiffnessMatrix(); } } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); auto & dof_manager = getDOFManager(); internal_force->zero(); computeStresses(); assembleInternalForce(); dof_manager.assembleToResidual("displacement", *internal_force, -1); dof_manager.assembleToResidual("displacement", *external_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes std::tuple StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce() { for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { assembleInternalForce(type, _not_ghost); // assembleInternalForce(type, _ghost); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce(ElementType type, GhostType gt) { auto & fem = getFEEngine(); auto & sigma = stress(type, gt); auto ndof = getNbDegreeOfFreedom(type); auto nb_nodes = mesh.getNbNodesPerElement(type); auto ndof_per_elem = ndof * nb_nodes; Array BtSigma(fem.getNbIntegrationPoints(type) * mesh.getNbElement(type), ndof_per_elem, "BtSigma"); fem.computeBtD(sigma, BtSigma, type, gt); Array intBtSigma(0, ndof_per_elem, "intBtSigma"); fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt); BtSigma.resize(0); getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force, type, gt, 1); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/element_info_per_processor_tmpl.hh b/src/synchronizer/element_info_per_processor_tmpl.hh index 080a71fb7..febb3e794 100644 --- a/src/synchronizer/element_info_per_processor_tmpl.hh +++ b/src/synchronizer/element_info_per_processor_tmpl.hh @@ -1,146 +1,146 @@ /** * @file element_info_per_processor_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Tue Feb 20 2018 * * @brief Helper classes to create the distributed synchronizer and distribute * a mesh * * * 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 "element_group.hh" #include "element_info_per_processor.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_INFO_PER_PROCESSOR_TMPL_HH_ #define AKANTU_ELEMENT_INFO_PER_PROCESSOR_TMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template void ElementInfoPerProc::fillMeshDataTemplated(BufferType & buffer, const std::string & tag_name, UInt nb_component) { AKANTU_DEBUG_ASSERT(this->mesh.getNbElement(this->type) == nb_local_element, "Did not got enought informations for the tag " << tag_name << " and the element type " << this->type << ":" << "_not_ghost." << " Got " << nb_local_element << " values, expected " << mesh.getNbElement(this->type)); mesh.getElementalData(tag_name); Array & data = mesh.getElementalDataArrayAlloc( tag_name, this->type, _not_ghost, nb_component); data.resize(nb_local_element); /// unpacking the data, element by element for (UInt i(0); i < nb_local_element; ++i) { for (UInt j(0); j < nb_component; ++j) { buffer >> data(i, j); } } AKANTU_DEBUG_ASSERT(mesh.getNbElement(this->type, _ghost) == nb_ghost_element, "Did not got enought informations for the tag " << tag_name << " and the element type " << this->type << ":" << "_ghost." << " Got " << nb_ghost_element << " values, expected " << mesh.getNbElement(this->type, _ghost)); Array & data_ghost = mesh.getElementalDataArrayAlloc( tag_name, this->type, _ghost, nb_component); data_ghost.resize(nb_ghost_element); /// unpacking the ghost data, element by element for (UInt j(0); j < nb_ghost_element; ++j) { for (UInt k(0); k < nb_component; ++k) { buffer >> data_ghost(j, k); } } } /* -------------------------------------------------------------------------- */ template void ElementInfoPerProc::fillMeshData(BufferType & buffer, const std::string & tag_name, const MeshDataTypeCode & type_code, UInt nb_component) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case MeshDataTypeCode::BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ fillMeshDataTemplated(buffer, tag_name, \ nb_component); \ break; \ } switch (type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not determine the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ template void ElementInfoPerProc::fillElementGroupsFromBuffer( CommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); Element el; el.type = type; for (auto ghost_type : ghost_types) { UInt nb_element = mesh.getNbElement(type, ghost_type); el.ghost_type = ghost_type; for (UInt e = 0; e < nb_element; ++e) { el.element = e; std::vector element_to_group; buffer >> element_to_group; AKANTU_DEBUG_ASSERT(e < mesh.getNbElement(type, ghost_type), "The mesh does not have the element " << e); - for (auto && element : element_to_group) { - mesh.getElementGroup(element).add(el, false, false); + for (auto && group : element_to_group) { + mesh.getElementGroup(group).add(el, false, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* AKANTU_ELEMENT_INFO_PER_PROCESSOR_TMPL_HH_ */ diff --git a/src/synchronizer/master_element_info_per_processor.cc b/src/synchronizer/master_element_info_per_processor.cc index 1421dbd4f..9df05f725 100644 --- a/src/synchronizer/master_element_info_per_processor.cc +++ b/src/synchronizer/master_element_info_per_processor.cc @@ -1,456 +1,453 @@ /** * @file master_element_info_per_processor.cc * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Tue Feb 20 2018 * * @brief Helper class to distribute a mesh * * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "communicator.hh" #include "element_group.hh" #include "element_info_per_processor.hh" #include "element_synchronizer.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MasterElementInfoPerProc::MasterElementInfoPerProc( ElementSynchronizer & synchronizer, UInt message_cnt, UInt root, ElementType type, const MeshPartition & partition) : ElementInfoPerProc(synchronizer, message_cnt, root, type), partition(partition), all_nb_local_element(nb_proc, 0), all_nb_ghost_element(nb_proc, 0), all_nb_element_to_send(nb_proc, 0) { Vector size(5); size(0) = (UInt)type; if (type != _not_defined) { nb_nodes_per_element = Mesh::getNbNodesPerElement(type); nb_element = mesh.getNbElement(type); const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { this->all_nb_local_element[partition_num(el)]++; for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { this->all_nb_ghost_element[*part]++; } this->all_nb_element_to_send[partition_num(el)] += ghost_partition.getNbCols(el) + 1; } /// tag info auto && tag_names = this->mesh.getTagNames(type); this->nb_tags = tag_names.size(); size(4) = nb_tags; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { size(1) = this->all_nb_local_element[p]; size(2) = this->all_nb_ghost_element[p]; size(3) = this->all_nb_element_to_send[p]; AKANTU_DEBUG_INFO( "Sending connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_sizes) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_sizes)); } else { this->nb_local_element = this->all_nb_local_element[p]; this->nb_ghost_element = this->all_nb_ghost_element[p]; } } } else { for (UInt p = 0; p < this->nb_proc; ++p) { if (p != this->root) { AKANTU_DEBUG_INFO( "Sending empty connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_sizes) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_sizes)); } } } } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizeConnectivities() { const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector> buffers(this->nb_proc); const auto & connectivities = this->mesh.getConnectivity(this->type, _not_ghost); /// copying the local connectivity for (auto && part_conn : zip(partition_num, make_view(connectivities, this->nb_nodes_per_element))) { auto && part = std::get<0>(part_conn); auto && conn = std::get<1>(part_conn); for (UInt i = 0; i < conn.size(); ++i) { buffers[part].push_back(conn[i]); } } /// copying the connectivity of ghost element for (auto && tuple : enumerate(make_view(connectivities, this->nb_nodes_per_element))) { auto && el = std::get<0>(tuple); auto && conn = std::get<1>(tuple); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { UInt proc = *part; for (UInt i = 0; i < conn.size(); ++i) { buffers[proc].push_back(conn[i]); } } } #ifndef AKANTU_NDEBUG for (auto p : arange(this->nb_proc)) { UInt size = this->nb_nodes_per_element * (this->all_nb_local_element[p] + this->all_nb_ghost_element[p]); AKANTU_DEBUG_ASSERT( buffers[p].size() == size, "The connectivity data packed in the buffer are not correct"); } #endif /// send all connectivity and ghost information to all processors std::vector requests; for (auto p : arange(this->nb_proc)) { if (p == this->root) { continue; } auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_connectivity); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } Array & old_nodes = this->getNodesGlobalIds(); /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, buffers[root], all_nb_local_element[root], all_nb_ghost_element[root], type, old_nodes); Communicator::waitAll(requests); Communicator::freeCommunicationRequest(requests); } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizePartitions() { const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector> buffers(this->partition.getNbPartition()); /// splitting the partition information to send them to processors Vector count_by_proc(nb_proc, 0); for (UInt el = 0; el < nb_element; ++el) { UInt proc = partition_num(el); buffers[proc].push_back(ghost_partition.getNbCols(el)); UInt i(0); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[proc].push_back(*part); } } for (UInt el = 0; el < nb_element; ++el) { UInt i(0); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[*part].push_back(partition_num(el)); } } #ifndef AKANTU_NDEBUG for (UInt p = 0; p < this->nb_proc; ++p) { AKANTU_DEBUG_ASSERT(buffers[p].size() == (this->all_nb_ghost_element[p] + this->all_nb_element_to_send[p]), "Data stored in the buffer are most probably wrong"); } #endif std::vector requests; /// last data to compute the communication scheme for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->root) { continue; } auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_partitions); AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } if (Mesh::getSpatialDimension(this->type) == this->mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); this->fillCommunicationScheme(buffers[this->rank]); } Communicator::waitAll(requests); Communicator::freeCommunicationRequest(requests); } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeTags() { AKANTU_DEBUG_IN(); if (this->nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } /// tag info auto tag_names = mesh.getTagNames(type); // Make sure the tags are sorted (or at least not in random order), // because they come from a map !! std::sort(tag_names.begin(), tag_names.end()); // Sending information about the tags in mesh_data: name, data type and // number of components of the underlying array associated to the current // type DynamicCommunicationBuffer mesh_data_sizes_buffer; for (auto && tag_name : tag_names) { mesh_data_sizes_buffer << tag_name; mesh_data_sizes_buffer << mesh.getTypeCode(tag_name); mesh_data_sizes_buffer << mesh.getNbComponent(tag_name, type); } AKANTU_DEBUG_INFO( "Broadcasting the size of the information about the mesh data tags: (" << mesh_data_sizes_buffer.size() << ")."); AKANTU_DEBUG_INFO( "Broadcasting the information about the mesh data tags, addr " << (void *)mesh_data_sizes_buffer.storage()); comm.broadcast(mesh_data_sizes_buffer, root); if (mesh_data_sizes_buffer.empty()) { return; } // Sending the actual data to each processor std::vector buffers(nb_proc); // Loop over each tag for the current type for (auto && tag_name : tag_names) { // Type code of the current tag (i.e. the tag named *names_it) this->fillTagBuffer(buffers, tag_name); } std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == root) { continue; } auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_mesh_data); AKANTU_DEBUG_INFO("Sending " << buffers[p].size() << " bytes of mesh data to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } // Loop over each tag for the current type for (auto && tag_name : tag_names) { // Reinitializing the mesh data on the master this->fillMeshData(buffers[root], tag_name, mesh.getTypeCode(tag_name), mesh.getNbComponent(tag_name, type)); } Communicator::waitAll(requests); Communicator::freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MasterElementInfoPerProc::fillTagBufferTemplated( std::vector & buffers, const std::string & tag_name) { const auto & data = mesh.getElementalDataArray(tag_name, type); const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); // Not possible to use the iterator because it potentially triggers the // creation of complex // type templates (such as akantu::Vector< std::vector > which don't // implement the right interface // (e.g. operator<< in that case). // typename Array::template const_iterator< Vector > data_it = // data.begin(data.getNbComponent()); // typename Array::template const_iterator< Vector > data_end = // data.end(data.getNbComponent()); const T * data_it = data.storage(); const T * data_end = data.storage() + data.size() * data.getNbComponent(); const UInt * part = partition_num.storage(); /// copying the data, element by element for (; data_it != data_end; ++part) { for (UInt j(0); j < data.getNbComponent(); ++j, ++data_it) { buffers[*part] << *data_it; } } data_it = data.storage(); /// copying the data for the ghost element for (UInt el(0); data_it != data_end; data_it += data.getNbComponent(), ++el) { auto it = ghost_partition.begin(el); auto end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; for (UInt j(0); j < data.getNbComponent(); ++j) { buffers[proc] << data_it[j]; } } } } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::fillTagBuffer( std::vector & buffers, const std::string & tag_name) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case MeshDataTypeCode::BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ this->fillTagBufferTemplated(buffers, \ tag_name); \ break; \ } MeshDataTypeCode data_type_code = mesh.getTypeCode(tag_name); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); std::vector buffers(nb_proc); using ElementToGroup = std::vector>; ElementToGroup element_to_group(nb_element); for (auto & eg : mesh.iterateElementGroups()) { const auto & name = eg.getName(); for (const auto & element : eg.getElements(type, _not_ghost)) { element_to_group[element].push_back(name); } - auto eit = eg.begin(type, _not_ghost); - if (eit != eg.end(type, _not_ghost)) { - const_cast &>(eg.getElements(type)).zero(); - } + eg.clear(type, _not_ghost); } const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); /// copying the data, element by element for (auto && pair : zip(partition_num, element_to_group)) { buffers[std::get<0>(pair)] << std::get<1>(pair); } /// copying the data for the ghost element for (auto && pair : enumerate(element_to_group)) { auto && el = std::get<0>(pair); auto it = ghost_partition.begin(el); auto end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; buffers[proc] << std::get<1>(pair); } } std::vector requests; for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->rank) { continue; } auto && tag = Tag::genTag(this->rank, p, Tag::_element_group); AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } this->fillElementGroupsFromBuffer(buffers[this->rank]); Communicator::waitAll(requests); Communicator::freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu