diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc index 849ddeb67..73cb90e31 100644 --- a/src/mesh/element_group.cc +++ b/src/mesh/element_group.cc @@ -1,205 +1,227 @@ /** * @file element_group.cc * * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Dec 09 2020 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ +#include "element_group.hh" #include "aka_csr.hh" #include "dumpable.hh" #include "dumpable_inline_impl.hh" #include "group_manager.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_utils.hh" +#if defined(AKANTU_COHESIVE_ELEMENT) +#include "cohesive_element_inserter.hh" +#endif #include #include #include - -#include "element_group.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" namespace akantu { /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const std::string & group_name, const Mesh & mesh, NodeGroup & node_group, UInt dimension, const std::string & id) : mesh(mesh), name(group_name), elements("elements", id), node_group(node_group), dimension(dimension) { AKANTU_DEBUG_IN(); this->registerDumper("paraview_" + group_name, group_name, true); this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), _all_dimensions); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const ElementGroup & /*other*/) = default; /* -------------------------------------------------------------------------- */ void ElementGroup::clear() { elements.free(); } /* -------------------------------------------------------------------------- */ void ElementGroup::clear(ElementType type, GhostType ghost_type) { if (elements.exists(type, ghost_type)) { elements(type, ghost_type).clear(); } } /* -------------------------------------------------------------------------- */ bool ElementGroup::empty() const { return elements.empty(); } /* -------------------------------------------------------------------------- */ void ElementGroup::append(const ElementGroup & other_group) { AKANTU_DEBUG_IN(); node_group.append(other_group.node_group); /// loop on all element types in all dimensions for (auto ghost_type : ghost_types) { for (auto type : other_group.elementTypes(_ghost_type = ghost_type, _element_kind = _ek_not_defined)) { const Array & other_elem_list = other_group.elements(type, ghost_type); UInt nb_other_elem = other_elem_list.size(); Array * elem_list; UInt nb_elem = 0; /// create current type if doesn't exists, otherwise get information if (elements.exists(type, ghost_type)) { elem_list = &elements(type, ghost_type); nb_elem = elem_list->size(); } else { elem_list = &(elements.alloc(0, 1, type, ghost_type)); } /// append new elements to current list elem_list->resize(nb_elem + nb_other_elem); std::copy(other_elem_list.begin(), other_elem_list.end(), elem_list->begin() + nb_elem); } } this->optimize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) { ; } stream << space << "ElementGroup [" << std::endl; stream << space << " + name: " << name << std::endl; stream << space << " + dimension: " << dimension << std::endl; elements.printself(stream, indent + 1); node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void ElementGroup::optimize() { // increasing the locality of data when iterating on the element of a group for (auto ghost_type : ghost_types) { for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { - Array & els = elements(type, ghost_type); + auto & els = elements(type, ghost_type); std::sort(els.begin(), els.end()); - Array::iterator<> end = std::unique(els.begin(), els.end()); + auto end = std::unique(els.begin(), els.end()); els.resize(end - els.begin()); } } node_group.optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::fillFromNodeGroup() { CSR node_to_elem; MeshUtils::buildNode2Elements(this->mesh, node_to_elem, this->dimension); std::set seen; - 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; + for (auto node : node_group) { + auto ite = node_to_elem.begin(node); + auto ende = node_to_elem.end(node); + for (auto && element : node_to_elem.getRow(node)) { if (this->dimension != _all_dimensions && - this->dimension != Mesh::getSpatialDimension(elem.type)) { + this->dimension != Mesh::getSpatialDimension(element.type)) { continue; } - if (seen.find(elem) != seen.end()) { + + if (seen.find(element) != 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]; + auto nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); + auto conn = mesh.getConnectivity(element); UInt count = 0; - for (UInt n = 0; n < conn.size(); ++n) { - count += - (this->node_group.getNodes().find(conn(n)) != UInt(-1) ? 1 : 0); + for (auto && n : conn) { + count += (this->node_group.getNodes().find(n) != UInt(-1) ? 1 : 0); } if (count == nb_nodes_per_element) { - this->add(elem); + this->add(element); } - seen.insert(elem); + seen.insert(element); } } this->optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::addDimension(UInt dimension) { this->dimension = std::max(dimension, this->dimension); } +/* -------------------------------------------------------------------------- */ +void ElementGroup::onNodesAdded(const Array & new_nodes, + const NewNodesEvent & event) { +#if defined(AKANTU_COHESIVE_ELEMENT) + if (aka::is_of_type(event)) { + // nodes might have changed in the connectivity + node_group.clear(); + for (auto ghost_type : ghost_types) { + for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { + auto & els = elements(type, ghost_type); + auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); + auto && conn_it = make_view(mesh.getConnectivity(type, ghost_type), + nb_nodes_per_element) + .begin(); + for (auto element : els) { + auto && conn = conn_it[element]; + for (auto && n : conn) { + node_group.add(n, false); + } + } + } + } + node_group.optimize(); + } +#endif +} /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 6b5cfd2a7..9824916b6 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,205 +1,208 @@ /** * @file element_group.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Mon Mar 08 2021 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.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; +class NewNodesEvent; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ElementGroup : 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"); 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); + void onNodesAdded(const Array & new_nodes, const NewNodesEvent & event); + 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)...); } decltype(auto) getElementsIterable(ElementType type, GhostType ghost_type = _not_ghost) const; // 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 e2a9598df..1c127862c 100644 --- a/src/mesh/group_manager.cc +++ b/src/mesh/group_manager.cc @@ -1,1034 +1,1043 @@ /** * @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: Thu Nov 12 2020 * * @brief Stores information about ElementGroup and NodeGroup * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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) : id(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()); 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"); 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"); 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); 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); 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)); } } } } 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); } } } } 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()); } +/* -------------------------------------------------------------------------- */ +void GroupManager::onNodesAdded(const Array & new_nodes, + const NewNodesEvent & event) { + for (auto && group : make_values_adaptor(element_groups)) { + group->onNodesAdded(new_nodes, event); + } +} +/* -------------------------------------------------------------------------- */ + } // namespace akantu diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index 40c30bde6..f19a7d966 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,349 +1,352 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Fri Jul 24 2020 * * @brief Stores information relevent to the notion of element and nodes * groups. * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_GROUP_MANAGER_HH_ #define AKANTU_GROUP_MANAGER_HH_ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_iterators.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { -class ElementGroup; -class NodeGroup; -class Mesh; class Element; +class ElementGroup; class ElementSynchronizer; +class Mesh; +class NewNodesEvent; +class NodeGroup; template class CommunicationBufferTemplated; namespace dumpers { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: using ElementGroups = std::map>; using NodeGroups = std::map>; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(Mesh & mesh, const ID & id = "group_manager"); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ [[deprecated( \ "use iterate(Element|Node)Groups or " \ "(element|node)GroupExists")]] inline BOOST_PP_CAT(BOOST_PP_CAT(const_, \ group_type), \ _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ [[deprecated("use iterate(Element|Node)Groups or " \ "(element|node)GroupExists")]] inline BOOST_PP_CAT(group_type, \ _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); public: decltype(auto) iterateNodeGroups() { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateNodeGroups() const { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateElementGroups() { return make_dereference_adaptor(make_values_adaptor(element_groups)); } decltype(auto) iterateElementGroups() const { return make_dereference_adaptor(make_values_adaptor(element_groups)); } /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* ------------------------------------------------------------------------ */ public: class ClusteringFilter { public: virtual bool operator()(const Element & /*unused*/) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /* ------------------------------------------------------------------------ */ /// renames an element group void renameElementGroup(const std::string & name, const std::string & new_name); /// renames a node group void renameNodeGroup(const std::string & name, const std::string & new_name); /// copy an existing element group void copyElementGroup(const std::string & name, const std::string & new_name); /// copy an existing node group void copyNodeGroup(const std::string & name, const std::string & new_name); /* ------------------------------------------------------------------------ */ /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); // /// destroy all element groups and the associated node groups // void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, const std::string & cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) template class dump_type> std::shared_ptr createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, ElementKind kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> std::shared_ptr createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, ElementKind kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> std::shared_ptr createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, ElementKind kind, ElementTypeMap nb_data_per_elem); template class ftype> std::shared_ptr createNodalField(const ftype * field, const std::string & group_name, UInt padding_size = 0); template class ftype> std::shared_ptr createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); + void onNodesAdded(const Array & new_nodes, const NewNodesEvent & event); + protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template inline std::shared_ptr createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, ElementKind kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template inline std::shared_ptr createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, ElementKind kind, ElementTypeMap nb_data_per_elem); /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: // AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; bool elementGroupExists(const std::string & name) { return element_groups.find(name) != element_groups.end(); } bool nodeGroupExists(const std::string & name) { return node_groups.find(name) != node_groups.end(); } private: template void renameGroup(GroupsType & groups, const std::string & name, const std::string & new_name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AKANTU_GROUP_MANAGER_HH_ */ diff --git a/src/mesh/mesh_inline_impl.hh b/src/mesh/mesh_inline_impl.hh index 22d27cfd4..6d9160394 100644 --- a/src/mesh/mesh_inline_impl.hh +++ b/src/mesh/mesh_inline_impl.hh @@ -1,765 +1,766 @@ /** * @file mesh_inline_impl.hh * * @author Guillaume Anciaux * @author Dana Christen * @author Mohit Pundir * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Thu Jul 15 2010 * @date last modification: Fri Dec 11 2020 * * @brief Implementation of the inline functions of the mesh class * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "element_class.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MESH_INLINE_IMPL_HH_ #define AKANTU_MESH_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline ElementKind Element::kind() const { return Mesh::getKind(type); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const { return connectivities.elementTypes(_pack...); } /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh, const std::string & origin) : MeshEvent(origin), new_numbering(mesh.getNbNodes(), 1, "new_numbering") {} /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh, const ID & new_numbering_id, const std::string & origin) : MeshEvent(origin), new_numbering(new_numbering_id, mesh.getID()) {} /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewElementsEvent & event) { this->fillNodesToElements(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewNodesEvent & event) { this->computeBoundingBox(); this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal); + GroupManager::onNodesAdded(event.getList(), event); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { this->connectivities.onElementsRemoved(event.getNewNumbering()); this->fillNodesToElements(); this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { const auto & new_numbering = event.getNewNumbering(); this->removeNodesFromArray(*nodes, new_numbering); if (nodes_global_ids and not is_mesh_facets) { this->removeNodesFromArray(*nodes_global_ids, new_numbering); } if (not is_mesh_facets) { this->removeNodesFromArray(*nodes_flags, new_numbering); } if (not nodes_to_elements.empty()) { std::vector>> tmp( nodes_to_elements.size()); auto it = nodes_to_elements.begin(); UInt new_nb_nodes = 0; for (auto new_i : new_numbering) { if (new_i != UInt(-1)) { tmp[new_i] = std::move(*it); ++new_nb_nodes; } ++it; } tmp.resize(new_nb_nodes); std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); } computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromArray(Array & vect, const Array & new_numbering) { Array tmp(vect.size(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.size(); ++i) { UInt new_i = new_numbering(i); if (new_i != UInt(-1)) { T * to_copy = vect.storage() + i * nb_component; std::uninitialized_copy(to_copy, to_copy + nb_component, tmp.storage() + new_i * nb_component); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if (not nodes_global_ids) { nodes_global_ids = std::make_shared>( nodes->size(), 1, getID() + ":nodes_global_ids"); for (auto && global_ids : enumerate(*nodes_global_ids)) { std::get<1>(global_ids) = std::get<0>(global_ids); } } AKANTU_DEBUG_OUT(); return *nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getConnectivityPointer(ElementType type, GhostType ghost_type) { if (connectivities.exists(type, ghost_type)) { return connectivities(type, ghost_type); } if (ghost_type != _not_ghost) { ghosts_counters.alloc(0, 1, type, ghost_type, 1); } AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); return connectivities.alloc(0, nb_nodes_per_element, type, ghost_type); } /* -------------------------------------------------------------------------- */ inline Array> & Mesh::getElementToSubelementPointer(ElementType type, GhostType ghost_type) { return getDataPointer>("element_to_subelement", type, ghost_type, 1, true); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getSubelementToElementPointer(ElementType type, GhostType ghost_type) { auto & array = getDataPointer( "subelement_to_element", type, ghost_type, getNbFacetsPerElement(type), false, is_mesh_facets, ElementNull); return array; } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement() const { return getData>("element_to_subelement"); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelementNC() { return getData>("element_to_subelement"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(ElementType type, GhostType ghost_type) const { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelementNC(ElementType type, GhostType ghost_type) { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const Element & element) const { return getData>("element_to_subelement")(element, 0); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelementNC(const Element & element) { return getData>("element_to_subelement")(element, 0); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement() const { return getData("subelement_to_element"); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getSubelementToElementNC() { return getData("subelement_to_element"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement(ElementType type, GhostType ghost_type) const { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getSubelementToElementNC(ElementType type, GhostType ghost_type) { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) const { return this->getSubelementToElement().get(element); } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElementNC(const Element & element) const { return this->getSubelementToElement().get(element); } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, ElementType el_type, GhostType ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent) { Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) { tmp.resize(mesh_parent->getNbElement(el_type, ghost_type)); } else { tmp.resize(this->getNbElement(el_type, ghost_type)); } } return tmp; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, ElementType el_type, GhostType ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_) { Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) { tmp.resize(mesh_parent->getNbElement(el_type, ghost_type), defaul_); } else { tmp.resize(this->getNbElement(el_type, ghost_type), defaul_); } } return tmp; } /* -------------------------------------------------------------------------- */ template inline const Array & Mesh::getData(const ID & data_name, ElementType el_type, GhostType ghost_type) const { return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getData(const ID & data_name, ElementType el_type, GhostType ghost_type) { return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline const ElementTypeMapArray & Mesh::getData(const ID & data_name) const { return this->getElementalData(data_name); } /* -------------------------------------------------------------------------- */ template inline ElementTypeMapArray & Mesh::getData(const ID & data_name) { return this->getElementalData(data_name); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(ElementType type, GhostType ghost_type) const { try { const Array & conn = connectivities(type, ghost_type); return conn.size(); } catch (...) { return 0; } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const UInt spatial_dimension, GhostType ghost_type, ElementKind kind) const { AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1), "spatial_dimension is " << spatial_dimension << " and is greater than 3 !"); UInt nb_element = 0; for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) { nb_element += getNbElement(type, ghost_type); } return nb_element; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector & barycenter) const { Vector conn = getConnectivity(element); Matrix local_coord(spatial_dimension, conn.size()); auto node_begin = make_view(*nodes, spatial_dimension).begin(); for (auto && node : enumerate(conn)) { local_coord(std::get<0>(node)) = Vector(node_begin[std::get<1>(node)]); } Math::barycenter(local_coord.storage(), conn.size(), spatial_dimension, barycenter.storage()); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(ElementType type) { UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT); #undef GET_NB_NODES_PER_ELEMENT return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(ElementType type) { ElementType p1_type = _not_defined; #define GET_P1_TYPE(type) p1_type = ElementClass::getP1ElementType() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE); #undef GET_P1_TYPE return p1_type; } /* -------------------------------------------------------------------------- */ inline ElementKind Mesh::getKind(ElementType type) { ElementKind kind = _ek_not_defined; #define GET_KIND(type) kind = ElementClass::getKind() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND); #undef GET_KIND return kind; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(ElementType type) { UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION); #undef GET_SPATIAL_DIMENSION return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNaturalSpaceDimension(const ElementType & type) { UInt natural_dimension = 0; #define GET_NATURAL_DIMENSION(type) \ natural_dimension = ElementClass::getNaturalSpaceDimension() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_DIMENSION); #undef GET_NATURAL_DIMENSION return natural_dimension; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetTypes(ElementType type, UInt /*t*/) { UInt nb = 0; #define GET_NB_FACET_TYPE(type) nb = ElementClass::getNbFacetTypes() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE); #undef GET_NB_FACET_TYPE return nb; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getFacetType(ElementType type, UInt t) { #define GET_FACET_TYPE(type) return ElementClass::getFacetType(t); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return _not_defined; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getAllFacetTypes(ElementType type) { #define GET_FACET_TYPE(type) return ElementClass::getFacetTypes(); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return ElementClass<_not_defined>::getFacetTypes(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(ElementType type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) n_facet = ElementClass::getNbFacetsPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(ElementType type, UInt t) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetLocalConnectivity(ElementType type, UInt t) { AKANTU_DEBUG_IN(); #define GET_FACET_CON(type) \ AKANTU_DEBUG_OUT(); \ return ElementClass::getFacetLocalConnectivityPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON); #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return ElementClass<_not_defined>::getFacetLocalConnectivityPerElement(0); // This avoid a compilation warning but will certainly // also cause a segfault if reached } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetConnectivity(const Element & element, UInt t) const { AKANTU_DEBUG_IN(); Matrix local_facets(getFacetLocalConnectivity(element.type, t)); Matrix facets(local_facets.rows(), local_facets.cols()); const Array & conn = connectivities(element.type, element.ghost_type); for (UInt f = 0; f < facets.rows(); ++f) { for (UInt n = 0; n < facets.cols(); ++n) { facets(f, n) = conn(element.element, local_facets(f, n)); } } AKANTU_DEBUG_OUT(); return facets; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) const { return connectivities.get(element); } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivityNC(const Element & element) { return connectivities.get(element); } /* -------------------------------------------------------------------------- */ template inline void Mesh::extractNodalValuesFromElement( const Array & nodal_values, T * local_coord, const UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const { for (UInt n = 0; n < n_nodes; ++n) { memcpy(local_coord + n * nb_degree_of_freedom, nodal_values.storage() + connectivity[n] * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(T)); } } /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(ElementType type, GhostType ghost_type) { getConnectivityPointer(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_pure_ghost; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_local_master_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_master; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_slave; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPeriodicSlave(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) == NodeFlag::_periodic_slave; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPeriodicMaster(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) == NodeFlag::_periodic_master; } /* -------------------------------------------------------------------------- */ inline NodeFlag Mesh::getNodeFlag(UInt local_id) const { return (*nodes_flags)(local_id); } /* -------------------------------------------------------------------------- */ inline Int Mesh::getNodePrank(UInt local_id) const { auto it = nodes_prank.find(local_id); return it == nodes_prank.end() ? -1 : it->second; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeLocalId(UInt global_id) const { if (nodes_global_ids == nullptr) { return global_id; } return nodes_global_ids->find(global_id); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->size(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElementList(const Array & elements) { UInt nb_nodes_per_element = 0; UInt nb_nodes = 0; ElementType current_element_type = _not_defined; for (const auto & el : elements) { if (el.type != current_element_type) { current_element_type = el.type; nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type); } nb_nodes += nb_nodes_per_element; } return nb_nodes; } /* -------------------------------------------------------------------------- */ inline Mesh & Mesh::getMeshFacets() { if (this->mesh_facets == nullptr) { AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); } return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshFacets() const { if (this->mesh_facets == nullptr) { AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); } return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshParent() const { if (this->mesh_parent == nullptr) { AKANTU_SILENT_EXCEPTION( "No parent mesh is defined! This is only valid in a mesh_facets"); } return *this->mesh_parent; } /* -------------------------------------------------------------------------- */ void Mesh::addPeriodicSlave(UInt slave, UInt master) { if (master == slave) { return; } // if pair already registered auto master_slaves = periodic_master_slave.equal_range(master); auto slave_it = std::find_if(master_slaves.first, master_slaves.second, [&](auto & pair) { return pair.second == slave; }); if (slave_it == master_slaves.second) { // no duplicates periodic_master_slave.insert(std::make_pair(master, slave)); AKANTU_DEBUG_INFO("adding periodic slave, slave gid:" << getNodeGlobalId(slave) << " [lid: " << slave << "]" << ", master gid:" << getNodeGlobalId(master) << " [lid: " << master << "]"); // std::cout << "adding periodic slave, slave gid:" << // getNodeGlobalId(slave) // << " [lid: " << slave << "]" // << ", master gid:" << getNodeGlobalId(master) // << " [lid: " << master << "]" << std::endl; } periodic_slave_master[slave] = master; auto set_flag = [&](auto node, auto flag) { (*nodes_flags)[node] &= ~NodeFlag::_periodic_mask; // clean periodic flags (*nodes_flags)[node] |= flag; }; set_flag(slave, NodeFlag::_periodic_slave); set_flag(master, NodeFlag::_periodic_master); } /* -------------------------------------------------------------------------- */ UInt Mesh::getPeriodicMaster(UInt slave) const { return periodic_slave_master.at(slave); } /* -------------------------------------------------------------------------- */ class Mesh::PeriodicSlaves { using internal_iterator = std::unordered_multimap::const_iterator; std::pair pair; public: PeriodicSlaves(const Mesh & mesh, UInt master) : pair(mesh.getPeriodicMasterSlaves().equal_range(master)) {} PeriodicSlaves(const PeriodicSlaves & other) = default; PeriodicSlaves(PeriodicSlaves && other) = default; PeriodicSlaves & operator=(const PeriodicSlaves & other) = default; class const_iterator { internal_iterator it; public: const_iterator(internal_iterator it) : it(it) {} const_iterator operator++() { ++it; return *this; } bool operator!=(const const_iterator & other) { return other.it != it; } bool operator==(const const_iterator & other) { return other.it == it; } auto operator*() { return it->second; } }; auto begin() const { return const_iterator(pair.first); } auto end() const { return const_iterator(pair.second); } }; /* -------------------------------------------------------------------------- */ inline decltype(auto) Mesh::getPeriodicSlaves(UInt master) const { return PeriodicSlaves(*this, master); } /* -------------------------------------------------------------------------- */ inline Vector Mesh::getConnectivityWithPeriodicity(const Element & element) const { Vector conn = getConnectivity(element); if (not isPeriodic()) { return conn; } for (auto && node : conn) { if (isPeriodicSlave(node)) { node = getPeriodicMaster(node); } } return conn; } } // namespace akantu #endif /* AKANTU_MESH_INLINE_IMPL_HH_ */ diff --git a/src/model/common/boundary_condition/boundary_condition_tmpl.hh b/src/model/common/boundary_condition/boundary_condition_tmpl.hh index 07ba189f1..ecbcfad58 100644 --- a/src/model/common/boundary_condition/boundary_condition_tmpl.hh +++ b/src/model/common/boundary_condition/boundary_condition_tmpl.hh @@ -1,233 +1,233 @@ /** * @file boundary_condition_tmpl.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Mon Oct 28 2019 * * @brief implementation of the applyBC * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "element_group.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_BOUNDARY_CONDITION_TMPL_HH_ #define AKANTU_BOUNDARY_CONDITION_TMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template void BoundaryCondition::initBC(ModelType & model, Array & primal, Array & dual) { this->model = &model; this->primal = &primal; this->dual = &dual; } /* -------------------------------------------------------------------------- */ template void BoundaryCondition::initBC(ModelType & model, Array & primal, Array & primal_increment, Array & dual) { this->initBC(model, primal, dual); this->primal_increment = &primal_increment; } /* -------------------------------------------------------------------------- */ /* Partial specialization for DIRICHLET functors */ template template struct BoundaryCondition::TemplateFunctionWrapper< FunctorType, BC::Functor::_dirichlet> { static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance) { auto & model = bc_instance.getModel(); auto & primal = bc_instance.getPrimal(); const auto & coords = model.getMesh().getNodes(); auto & boundary_flags = model.getBlockedDOFs(); UInt dim = model.getMesh().getSpatialDimension(); auto primal_iter = primal.begin(primal.getNbComponent()); auto coords_iter = coords.begin(dim); auto flags_iter = boundary_flags.begin(boundary_flags.getNbComponent()); for (auto n : group.getNodeGroup()) { Vector flag(flags_iter[n]); Vector primal(primal_iter[n]); Vector coords(coords_iter[n]); func(n, flag, primal, coords); } } }; /* -------------------------------------------------------------------------- */ /* Partial specialization for NEUMANN functors */ template template struct BoundaryCondition::TemplateFunctionWrapper< FunctorType, BC::Functor::_neumann> { static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance) { UInt dim = bc_instance.getModel().getSpatialDimension(); switch (dim) { case 1: { AKANTU_TO_IMPLEMENT(); break; } case 2: case 3: { applyBC(func, group, bc_instance, _not_ghost); applyBC(func, group, bc_instance, _ghost); break; } } } static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance, GhostType ghost_type) { auto & model = bc_instance.getModel(); auto & dual = bc_instance.getDual(); const auto & mesh = model.getMesh(); const auto & nodes_coords = mesh.getNodes(); const auto & fem_boundary = model.getFEEngineBoundary(); UInt dim = model.getSpatialDimension(); UInt nb_degree_of_freedom = dual.getNbComponent(); IntegrationPoint quad_point; quad_point.ghost_type = ghost_type; // Loop over the boundary element types for (auto && type : group.elementTypes(dim - 1, ghost_type)) { const auto & element_ids = group.getElements(type, ghost_type); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type, ghost_type); UInt nb_elements = element_ids.size(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array dual_before_integ(nb_elements * nb_quad_points, nb_degree_of_freedom, 0.); Array quad_coords(nb_elements * nb_quad_points, dim); const auto & normals_on_quad = fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type); fem_boundary.interpolateOnIntegrationPoints( nodes_coords, quad_coords, dim, type, ghost_type, element_ids); auto normals_begin = normals_on_quad.begin(dim); decltype(normals_begin) normals_iter; auto quad_coords_iter = quad_coords.begin(dim); auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom); quad_point.type = type; for (auto el : element_ids) { quad_point.element = el; normals_iter = normals_begin + el * nb_quad_points; for (auto q : arange(nb_quad_points)) { quad_point.num_point = q; func(quad_point, *dual_iter, *quad_coords_iter, *normals_iter); ++dual_iter; ++quad_coords_iter; ++normals_iter; } } Array dual_by_shapes(nb_elements * nb_quad_points, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, ghost_type, element_ids); Array dual_by_shapes_integ(nb_elements, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ, nb_degree_of_freedom * nb_nodes_per_element, type, ghost_type, element_ids); // assemble the result into force vector model.getDOFManager().assembleElementalArrayLocalArray( dual_by_shapes_integ, dual, type, ghost_type, 1., element_ids); } } }; /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func) { auto bit = model->getMesh().getGroupManager().element_group_begin(); auto bend = model->getMesh().getGroupManager().element_group_end(); for (; bit != bend; ++bit) { applyBC(func, *bit); } } /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func, const std::string & group_name) { try { const ElementGroup & element_group = model->getMesh().getElementGroup(group_name); applyBC(func, element_group); } catch (akantu::debug::Exception & e) { AKANTU_EXCEPTION("Error applying a boundary condition onto \"" << group_name << "\"! [" << e.what() << "]"); } } /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func, const ElementGroup & element_group) { #if !defined(AKANTU_NDEBUG) if (element_group.getDimension() != model->getSpatialDimension() - 1) { - AKANTU_DEBUG_WARNING("The group " - << element_group.getName() - << " does not contain only boundaries elements"); + AKANTU_DEBUG_INFO("The group " + << element_group.getName() + << " does not contain only boundaries elements"); } #endif TemplateFunctionWrapper::applyBC(func, element_group, *this); } #endif /* AKANTU_BOUNDARY_CONDITION_TMPL_HH_ */ } // namespace akantu diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 39ef230c1..817c09f72 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,704 +1,704 @@ /** * @file contact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Jul 28 2021 * * @brief Contact mechanics model * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "contact_mechanics_model.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, - Model::spatial_dimension, _not_ghost, + Model::spatial_dimension - 1, _not_ghost, _ek_regular); this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, id + ":contact_detector"); registerFEEngineObject("ContactFacetsFEEngine", mesh, Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (not this->parser.getLastParsedFile().empty()) { this->instantiateResolutions(); this->initResolutions(); } this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::instantiateResolutions() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); for (const auto & section : model_resolutions) { this->registerNewResolution(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); for (const auto & section : sub_sections) { this->registerNewResolution(section); } if (resolutions.empty()) { AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" << getID()); } are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ParserSection & section) { std::string res_name; std::string res_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); res_name = tmp; /** this can seem weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A contact resolution of type \'" << res_type << "\' in the input file has been defined without a name!"); } Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); res.parseSection(section); return res; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution( const ID & res_name, const ID & res_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == resolutions_names_to_id.end(), "A resolution with this name '" << res_name << "' has already been registered. " << "Please use unique names for resolutions"); UInt res_count = resolutions.size(); resolutions_names_to_id[res_name] = res_count; std::stringstream sstr_res; sstr_res << this->id << ":" << res_count << ":" << res_type; ID res_id = sstr_res.str(); std::unique_ptr resolution = ResolutionFactory::getInstance().allocate(res_type, spatial_dimension, opt_param, *this, res_id); resolutions.push_back(std::move(resolution)); return *(resolutions.back()); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initResolutions() { AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); if (!are_resolutions_instantiated) { instantiateResolutions(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { AKANTU_DEBUG_IN(); getFEEngine("ContactMechanicsModel").initShapeFunctions(_not_ghost); getFEEngine("ContactMechanicsModel").initShapeFunctions(_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver( TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType /*unused*/) { // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->normal_force, spatial_dimension, "normal_force"); this->allocNodalField(this->tangential_force, spatial_dimension, "tangential_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->contact_state, 1, "contact_state"); this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements"); this->allocNodalField(this->normals, spatial_dimension, "normals"); auto surface_dimension = spatial_dimension - 1; this->allocNodalField(this->tangents, surface_dimension * spatial_dimension, "tangents"); this->allocNodalField(this->projections, surface_dimension, "projections"); this->allocNodalField(this->previous_projections, surface_dimension, "previous_projections"); this->allocNodalField(this->previous_tangents, surface_dimension * spatial_dimension, "previous_tangents"); this->allocNodalField(this->tangential_tractions, surface_dimension, "tangential_tractions"); this->allocNodalField(this->previous_tangential_tractions, surface_dimension, "previous_tangential_tractions"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } 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 ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; 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"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); UInt nb_nodes = mesh.getNbNodes(); this->internal_force->clear(); this->normal_force->clear(); this->tangential_force->clear(); internal_force->resize(nb_nodes, 0.); normal_force->resize(nb_nodes, 0.); tangential_force->resize(nb_nodes, 0.); // assemble the forces due to contact auto assemble = [&](auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements // AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); // assemble(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->resize(nb_nodes); internal_array->zero(); }; resize_arrays(gaps); resize_arrays(normals); resize_arrays(tangents); resize_arrays(projections); resize_arrays(tangential_tractions); resize_arrays(contact_state); resize_arrays(nodal_area); resize_arrays(external_force); this->detector->search(contact_elements, *gaps, *normals, *tangents, *projections); - // intepenetration value must be positive for contact mechanics + // interpenetration value must be positive for contact mechanics // model to work by default the gap value from detector is negative std::for_each((*gaps).begin(), (*gaps).end(), [](Real & gap) { gap *= -1.; }); if (!contact_elements.empty()) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections (*previous_projections).copy(*projections); // saving previous tangents (*previous_tangents).copy(*tangents); // saving previous tangential traction (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); previous_master_elements->set(ElementNull); for (auto & element : contact_elements) { (*previous_master_elements)[element.slave] = element.master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas(GhostType ghost_type) { UInt nb_nodes = mesh.getNbNodes(); nodal_area->resize(nb_nodes); nodal_area->zero(); external_force->resize(nb_nodes); external_force->zero(); auto & fem_boundary = getFEEngineClassBoundary("ContactMechanicsModel"); fem_boundary.initShapeFunctions(getPositions(), _not_ghost); fem_boundary.initShapeFunctions(getPositions(), _ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); IntegrationPoint quad_point; quad_point.ghost_type = ghost_type; auto & group = mesh.getElementGroup("contact_surface"); UInt nb_degree_of_freedom = external_force->getNbComponent(); for (auto && type : group.elementTypes(spatial_dimension - 1, ghost_type)) { const auto & element_ids = group.getElements(type, ghost_type); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type, ghost_type); UInt nb_elements = element_ids.size(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array dual_before_integ(nb_elements * nb_quad_points, nb_degree_of_freedom, 0.); Array quad_coords(nb_elements * nb_quad_points, spatial_dimension); const auto & normals_on_quad = fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type); auto normals_begin = normals_on_quad.begin(spatial_dimension); decltype(normals_begin) normals_iter; auto quad_coords_iter = quad_coords.begin(spatial_dimension); auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom); quad_point.type = type; Element subelement; subelement.type = type; subelement.ghost_type = ghost_type; for (auto el : element_ids) { subelement.element = el; const auto & element_to_subelement = mesh.getElementToSubelement(type)(el); Vector outside(spatial_dimension); mesh.getBarycenter(subelement, outside); Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector inside_to_outside(spatial_dimension); inside_to_outside = outside - inside; normals_iter = normals_begin + el * nb_quad_points; quad_point.element = el; for (auto q : arange(nb_quad_points)) { quad_point.num_point = q; auto ddot = inside_to_outside.dot(*normals_iter); Vector normal(*normals_iter); if (ddot < 0) { normal *= -1.0; } (*dual_iter) .mul(Matrix::eye(spatial_dimension, 1), normal); ++dual_iter; ++quad_coords_iter; ++normals_iter; } } Array dual_by_shapes(nb_elements * nb_quad_points, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, ghost_type, element_ids); Array dual_by_shapes_integ(nb_elements, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ, nb_degree_of_freedom * nb_nodes_per_element, type, ghost_type, element_ids); this->getDOFManager().assembleElementalArrayLocalArray( dual_by_shapes_integ, *external_force, type, ghost_type, 1., element_ids); } for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); Vector force(std::get<1>(tuple)); area = force.norm(); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Contact Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + resolutions [" << std::endl; for (const auto & resolution : resolutions) { resolution->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { for (auto & resolution : resolutions) { resolution->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep(bool converged) { for (auto & resolution : resolutions) { resolution->afterSolveStep(converged); } } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string & /*unused*/, const std::string & /*unused*/, bool /*unused*/) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["contact_force"] = this->internal_force.get(); real_nodal_fields["normal_force"] = this->normal_force.get(); real_nodal_fields["tangential_force"] = this->tangential_force.get(); real_nodal_fields["blocked_dofs"] = this->blocked_dofs.get(); real_nodal_fields["normals"] = this->normals.get(); real_nodal_fields["tangents"] = this->tangents.get(); real_nodal_fields["gaps"] = this->gaps.get(); real_nodal_fields["areas"] = this->nodal_area.get(); real_nodal_fields["tangential_traction"] = this->tangential_tractions.get(); std::shared_ptr field; if (padding_flag) { field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); } else { field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); } return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldUInt(const std::string & field_name, const std::string & group_name, bool /*padding_flag*/) { std::shared_ptr field; if (field_name == "contact_state") { auto && func = std::make_unique>(); field = mesh.createNodalField(this->contact_state.get(), group_name); field = dumpers::FieldComputeProxy::createFieldCompute(field, std::move(func)); } return field; } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & elements, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & dofs, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } } // namespace akantu