diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 40cd5826c..1eb483f1c 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,202 +1,201 @@ /** * @file element_group.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Wed Nov 08 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element_type_map.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_GROUP_HH_ #define AKANTU_ELEMENT_GROUP_HH_ namespace akantu { class Mesh; class Element; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ElementGroup : private Memory, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ElementGroup(const std::string & name, const Mesh & mesh, NodeGroup & node_group, UInt dimension = _all_dimensions, const std::string & id = "element_group", const MemoryID & mem_id = 0); ElementGroup(const ElementGroup & /*unused*/); /* ------------------------------------------------------------------------ */ /* Type definitions */ /* ------------------------------------------------------------------------ */ public: using ElementList = ElementTypeMapArray; using NodeList = Array; /* ------------------------------------------------------------------------ */ /* Element iterator */ /* ------------------------------------------------------------------------ */ using type_iterator = ElementList::type_iterator; [[deprecated("Use elementTypes instead")]] inline type_iterator - firstType(UInt dim = _all_dimensions, - GhostType ghost_type = _not_ghost, + 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, + 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; + 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(); 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); + GhostType ghost_type = _not_ghost, bool add_nodes = true, + bool check_for_duplicate = true); inline void addNode(UInt node_id, bool check_for_duplicate = true); inline void removeNode(UInt node_id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the elements based on the underlying node group. virtual void fillFromNodeGroup(); // sort and remove duplicated values void optimize(); /// change the dimension if needed void addDimension(UInt dimension); private: inline void addElement(ElementType elem_type, UInt elem_id, GhostType ghost_type); friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - const Array & - getElements(ElementType type, - GhostType ghost_type = _not_ghost) const; + const Array & getElements(ElementType type, + GhostType ghost_type = _not_ghost) const; AKANTU_GET_MACRO(Elements, elements, const ElementTypeMapArray &); AKANTU_GET_MACRO_NOT_CONST(Elements, elements, ElementTypeMapArray &); + template auto size(Args &&... pack) const { + return elements.size(std::forward(pack)...); + } + // AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); AKANTU_GET_MACRO(NodeGroup, node_group, const NodeGroup &); AKANTU_GET_MACRO_NOT_CONST(NodeGroup, node_group, NodeGroup &); AKANTU_GET_MACRO(Dimension, dimension, UInt); AKANTU_GET_MACRO(Name, name, std::string); inline UInt getNbNodes() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Mesh to which this group belongs const Mesh & mesh; /// name of the group std::string name; /// list of elements composing the group ElementList elements; /// sub list of nodes which are composing the elements NodeGroup & node_group; /// group dimension UInt dimension{_all_dimensions}; /// empty arry for the iterator to work when an element type not present Array empty_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ElementGroup & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "element.hh" #include "element_group_inline_impl.hh" #endif /* AKANTU_ELEMENT_GROUP_HH_ */ diff --git a/src/mesh_utils/cohesive_element_inserter.cc b/src/mesh_utils/cohesive_element_inserter.cc index b6c8476b8..938ac5c35 100644 --- a/src/mesh_utils/cohesive_element_inserter.cc +++ b/src/mesh_utils/cohesive_element_inserter.cc @@ -1,324 +1,324 @@ /** * @file cohesive_element_inserter.cc * * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Mon Feb 19 2018 * * @brief Cohesive element inserter functions * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "cohesive_element_inserter.hh" #include "cohesive_element_inserter_helper.hh" #include "communicator.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "global_ids_updater.hh" #include "mesh_accessor.hh" #include "mesh_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, const ID & id) : Parsable(ParserType::_cohesive_inserter), id(id), mesh(mesh), mesh_facets(mesh.initMeshFacets()), insertion_facets("insertion_facets", id), insertion_limits(mesh.getSpatialDimension(), 2), check_facets("check_facets", id) { this->registerParam("cohesive_surfaces", physical_surfaces, _pat_parsable, "List of groups to consider for insertion"); this->registerParam("cohesive_zones", physical_zones, _pat_parsable, "List of groups to consider for insertion"); this->registerParam("bounding_box", insertion_limits, _pat_parsable, "Global limit for insertion"); UInt spatial_dimension = mesh.getSpatialDimension(); /// init insertion limits for (UInt dim = 0; dim < spatial_dimension; ++dim) { insertion_limits(dim, 0) = std::numeric_limits::max() * Real(-1.); insertion_limits(dim, 1) = std::numeric_limits::max(); } insertion_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = false); } /* -------------------------------------------------------------------------- */ CohesiveElementInserter::~CohesiveElementInserter() = default; /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::parseSection(const ParserSection & section) { Parsable::parseSection(section); if (is_extrinsic) { limitCheckFacets(this->check_facets); } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets() { limitCheckFacets(this->check_facets); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::setLimit(SpatialDirection axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_ASSERT( axis < mesh.getSpatialDimension(), "You are trying to limit insertion in a direction that doesn't exist"); insertion_limits(axis, 0) = std::min(first_limit, second_limit); insertion_limits(axis, 1) = std::max(first_limit, second_limit); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertIntrinsicElements() { limitCheckFacets(insertion_facets); return insertElements(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets( ElementTypeMapArray & check_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); check_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = true); check_facets.set(true); // remove the pure ghost elements for_each_element( mesh_facets, [&](auto && facet) { const auto & element_to_facet = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); auto & left = element_to_facet[0]; auto & right = element_to_facet[1]; if (right == ElementNull || (left.ghost_type == _ghost && right.ghost_type == _ghost)) { check_facets(facet) = false; return; } #ifndef AKANTU_NDEBUG if (left == ElementNull) { AKANTU_DEBUG_WARNING("By convention element should not have " "ElementNull on there first side: " << facet); } #endif if (left.kind() == _ek_cohesive or right.kind() == _ek_cohesive) { check_facets(facet) = false; } }, _spatial_dimension = spatial_dimension - 1); auto tolerance = Math::getTolerance(); Vector bary_facet(spatial_dimension); // set the limits to the bounding box for_each_element( mesh_facets, [&](auto && facet) { - auto & need_check = check_facets(facet, 0); + auto & need_check = check_facets(facet); if (not need_check) { return; } mesh_facets.getBarycenter(facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension and bary_facet(coord_in_limit) > (insertion_limits(coord_in_limit, 0) - tolerance) and bary_facet(coord_in_limit) < (insertion_limits(coord_in_limit, 1) + tolerance)) { ++coord_in_limit; } if (coord_in_limit != spatial_dimension) { need_check = false; } }, _spatial_dimension = spatial_dimension - 1); // remove the physical zones if (mesh.hasData("physical_names") and not physical_zones.empty()) { auto && physical_names = mesh.getData("physical_names"); for_each_element( mesh_facets, [&](auto && facet) { const auto & element_to_facet = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); auto count = 0; for (auto i : arange(2)) { const auto & element = element_to_facet[i]; if (element == ElementNull) { continue; } const auto & name = physical_names(element); count += find(physical_zones.begin(), physical_zones.end(), name) != physical_zones.end(); } if (count != 2) { check_facets(facet) = false; } }, _spatial_dimension = spatial_dimension - 1); } if (physical_surfaces.empty()) { AKANTU_DEBUG_OUT(); return; } if (not mesh_facets.hasData("physical_names")) { AKANTU_DEBUG_ASSERT( physical_surfaces.empty(), "No physical names in the mesh but insertion limited to a group"); AKANTU_DEBUG_OUT(); return; } const auto & physical_ids = mesh_facets.getData("physical_names"); // set the limits to the physical surfaces for_each_element( mesh_facets, [&](auto && facet) { auto & need_check = check_facets(facet, 0); if (not need_check) { return; } const auto & physical_id = physical_ids(facet); auto it = find(physical_surfaces.begin(), physical_surfaces.end(), physical_id); need_check = (it != physical_surfaces.end()); }, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertElements(bool only_double_facets) { CohesiveNewNodesEvent node_event(AKANTU_CURRENT_FUNCTION); NewElementsEvent element_event(AKANTU_CURRENT_FUNCTION); if (mesh_facets.isDistributed()) { mesh_facets.getElementSynchronizer().synchronizeOnce( *this, SynchronizationTag::_ce_groups); } CohesiveElementInserterHelper cohesive_element_inserter_helper( mesh, insertion_facets); UInt nb_new_elements{0}; if (only_double_facets) { nb_new_elements = cohesive_element_inserter_helper.insertFacetsOnly(); } else { nb_new_elements = cohesive_element_inserter_helper.insertCohesiveElement(); element_event.getList().copy( cohesive_element_inserter_helper.getNewElements()); } auto && doubled_nodes = cohesive_element_inserter_helper.getDoubledNodes(); auto nb_new_nodes = doubled_nodes.size(); node_event.getList().reserve(nb_new_nodes); node_event.getOldNodesList().reserve(nb_new_nodes); for (auto && doubled_node : make_view(doubled_nodes, 2)) { node_event.getList().push_back(doubled_node(1)); node_event.getOldNodesList().push_back(doubled_node(0)); } if (nb_new_elements > 0) { updateInsertionFacets(); } MeshAccessor mesh_accessor(mesh); std::tie(nb_new_nodes, nb_new_elements) = mesh_accessor.updateGlobalData(node_event, element_event); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::updateInsertionFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (auto && facet_gt : ghost_types) { for (auto && facet_type : mesh_facets.elementTypes(spatial_dimension - 1, facet_gt)) { auto & ins_facets = insertion_facets(facet_type, facet_gt); // this is the intrinsic case if (not is_extrinsic) { continue; } auto & f_check = check_facets(facet_type, facet_gt); for (auto && pair : zip(ins_facets, f_check)) { bool & ins = std::get<0>(pair); bool & check = std::get<1>(pair); if (ins) { ins = check = false; } } } } // resize for the newly added facets insertion_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = false); // resize for the newly added facets if (is_extrinsic) { check_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = false); } else { insertion_facets.set(false); } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/mesh_utils/cohesive_element_inserter_helper.cc b/src/mesh_utils/cohesive_element_inserter_helper.cc index d3ca59b51..dd304a2a2 100644 --- a/src/mesh_utils/cohesive_element_inserter_helper.cc +++ b/src/mesh_utils/cohesive_element_inserter_helper.cc @@ -1,914 +1,936 @@ /** * @file cohesive_element_inserter_helper.cc * * @author Nicolas Richart * * @date creation jeu sep 03 2020 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "cohesive_element_inserter_helper.hh" /* -------------------------------------------------------------------------- */ #include "data_accessor.hh" #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ CohesiveElementInserterHelper::CohesiveElementInserterHelper( Mesh & mesh, const ElementTypeMapArray & facet_insertion) : doubled_nodes(0, 2, "doubled_nodes"), mesh(mesh), mesh_facets(mesh.getMeshFacets()) { auto spatial_dimension = mesh_facets.getSpatialDimension(); for (auto gt : ghost_types) { for (auto type : mesh_facets.elementTypes(_ghost_type = gt)) { nb_new_facets(type, gt) = mesh_facets.getNbElement(type, gt); } } std::array nb_facet_to_insert{0, 0}; // creates a vector of the facets to insert std::vector potential_facets_to_double; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { const auto & f_insertion = facet_insertion(type_facet, gt_facet); auto & counter = nb_facet_to_insert[gt_facet == _not_ghost ? 0 : 1]; for (auto && data : enumerate(f_insertion)) { if (std::get<1>(data)) { UInt el = std::get<0>(data); potential_facets_to_double.push_back({type_facet, el, gt_facet}); ++counter; } } } } // Defines a global order of insertion oof cohesive element to ensure node // doubling appens in the smae order, this is necessary for the global node // numbering if (mesh_facets.isDistributed()) { const auto & comm = mesh_facets.getCommunicator(); ElementTypeMapArray global_orderings; global_orderings.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = -1); auto starting_index = nb_facet_to_insert[0]; comm.exclusiveScan(starting_index); // define the local numbers for facet to insert for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { for (auto data : zip(facet_insertion(type_facet, gt_facet), global_orderings(type_facet, gt_facet))) { if (std::get<0>(data)) { std::get<1>(data) = starting_index; ++starting_index; } } } } // retreives the oorder number from neighoring processors auto && synchronizer = mesh_facets.getElementSynchronizer(); SimpleElementDataAccessor data_accessor( global_orderings, SynchronizationTag::_ce_insertion_order); synchronizer.synchronizeOnce(data_accessor, SynchronizationTag::_ce_insertion_order); // sort the facets to double based on this global ordering std::sort(potential_facets_to_double.begin(), potential_facets_to_double.end(), [&global_orderings](auto && el1, auto && el2) { return global_orderings(el1) < global_orderings(el2); }); } for (auto d : arange(spatial_dimension)) { facets_to_double_by_dim[d] = std::make_unique>( 0, 2, "facets_to_double_" + std::to_string(d)); } auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; MeshAccessor mesh_accessor(mesh_facets); auto & elements_to_subelements = mesh_accessor.getElementToSubelement(); for (auto && facet_to_double : potential_facets_to_double) { auto gt_facet = facet_to_double.ghost_type; auto type_facet = facet_to_double.type; auto & elements_to_facets = elements_to_subelements(type_facet, gt_facet); auto & elements_to_facet = elements_to_facets(facet_to_double.element); if (elements_to_facet[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || elements_to_facet[1].kind() == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } auto new_facet = nb_new_facets(type_facet, gt_facet)++; facets_to_double.push_back(Vector{ facet_to_double, Element{type_facet, new_facet, gt_facet}}); /// update facet_to_element vector auto & element_to_update = elements_to_facet[1]; auto el = element_to_update.element; auto & facets_to_elements = mesh_facets.getSubelementToElement( element_to_update.type, element_to_update.ghost_type); auto facets_to_element = Vector( make_view(facets_to_elements, facets_to_elements.getNbComponent()) .begin()[el]); auto facet_to_update = std::find(facets_to_element.begin(), facets_to_element.end(), facet_to_double); AKANTU_DEBUG_ASSERT(facet_to_update != facets_to_element.end(), "Facet not found"); facet_to_update->element = new_facet; /// update elements connected to facet const auto & first_facet_list = elements_to_facet; elements_to_facets.push_back(first_facet_list); /// set new and original facets as boundary facets elements_to_facets(new_facet)[0] = elements_to_facets(new_facet)[1]; elements_to_facets(new_facet)[1] = ElementNull; elements_to_facets(facet_to_double.element)[1] = ElementNull; } } /* -------------------------------------------------------------------------- */ -inline bool +inline auto CohesiveElementInserterHelper::hasElement(const Vector & nodes_element, - const Vector & nodes) { + const Vector & nodes) -> bool { // if one of the nodes of "nodes" is not in "nodes_element" it stops auto it = std::mismatch(nodes.begin(), nodes.end(), nodes_element.begin(), [&](auto && node, auto && /*node2*/) -> bool { auto it = std::find(nodes_element.begin(), nodes_element.end(), node); return (it != nodes_element.end()); }); // true if all "nodes" where found in "nodes_element" return (it.first == nodes.end()); } /* -------------------------------------------------------------------------- */ - -/* -------------------------------------------------------------------------- */ -inline bool CohesiveElementInserterHelper::removeElementsInVector( +inline auto CohesiveElementInserterHelper::removeElementsInVector( const std::vector & elem_to_remove, - std::vector & elem_list) { - if (elem_list.size() <= elem_to_remove.size()) + std::vector & elem_list) -> bool { + if (elem_list.size() <= elem_to_remove.size()) { return true; + } auto el_it = elem_to_remove.begin(); auto el_last = elem_to_remove.end(); std::vector::iterator el_del; UInt deletions = 0; for (; el_it != el_last; ++el_it) { el_del = std::find(elem_list.begin(), elem_list.end(), *el_it); if (el_del != elem_list.end()) { elem_list.erase(el_del); ++deletions; } } AKANTU_DEBUG_ASSERT(deletions == 0 || deletions == elem_to_remove.size(), "Not all elements have been erased"); return deletions == 0; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, const std::vector * facet_list) { AKANTU_DEBUG_IN(); for (auto & element : element_list) { - if (element.type == _not_defined) + if (element.type == _not_defined) { continue; + } Vector connectivity = mesh.getConnectivity(element); if (element.kind() == _ek_cohesive) { AKANTU_DEBUG_ASSERT( facet_list != nullptr, "Provide a facet list in order to update cohesive elements"); const Vector facets = mesh_facets.getSubelementToElement(element); auto facet_nb_nodes = connectivity.size() / 2; /// loop over cohesive element's facets for (const auto & facet : enumerate(facets)) { /// skip facets if not present in the list if (std::find(facet_list->begin(), facet_list->end(), std::get<1>(facet)) == facet_list->end()) { continue; } auto n = std::get<0>(facet); - auto begin = connectivity.begin() + n * facet_nb_nodes; + auto begin = connectivity.begin() + static_cast(n * facet_nb_nodes); auto end = begin + facet_nb_nodes; auto it = std::find(begin, end, old_node); AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element"); *it = new_node; } } else { auto it = std::find(connectivity.begin(), connectivity.end(), old_node); AKANTU_DEBUG_ASSERT(it != connectivity.end(), "Node not found in current element"); /// update connectivity *it = new_node; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateSubelementToElement(UInt dim, bool facet_mode) { auto & facets_to_double = *facets_to_double_by_dim[dim]; auto & facets_to_subfacets = - elementsOfDimToElementsOfDim(dim + facet_mode, dim); + elementsOfDimToElementsOfDim(dim + static_cast(facet_mode), dim); for (auto && data : zip(make_view(facets_to_double, 2), facets_to_subfacets)) { const auto & old_subfacet = std::get<0>(data)[0]; const auto & new_subfacet = std::get<0>(data)[1]; auto & facet_to_subfacets = std::get<1>(data); MeshAccessor mesh_accessor(mesh_facets); for (auto & facet : facet_to_subfacets) { Vector subfacets = mesh_accessor.getSubelementToElement(facet); auto && subfacet_to_update_it = std::find(subfacets.begin(), subfacets.end(), old_subfacet); AKANTU_DEBUG_ASSERT(subfacet_to_update_it != subfacets.end(), "Subfacet not found"); *subfacet_to_update_it = new_subfacet; } } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateElementToSubelement(UInt dim, bool facet_mode) { auto & facets_to_double = *facets_to_double_by_dim[dim]; - auto & facets_to_subfacets = - elementsOfDimToElementsOfDim(dim + facet_mode, dim); + auto & facets_to_subfacets = elementsOfDimToElementsOfDim( + dim + static_cast(facet_mode), dim); MeshAccessor mesh_accessor(mesh_facets); // resize the arrays mesh_accessor.getElementToSubelement().initialize( mesh_facets, _spatial_dimension = dim, _with_nb_element = true); for (auto && data : zip(make_view(facets_to_double, 2), facets_to_subfacets)) { const auto & new_facet = std::get<0>(data)[1]; mesh_accessor.getElementToSubelement(new_facet) = std::get<1>(data); } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateQuadraticSegments(UInt dim) { AKANTU_DEBUG_IN(); auto spatial_dimension = mesh.getSpatialDimension(); auto & facets_to_double = *facets_to_double_by_dim[dim]; MeshAccessor mesh_accessor(mesh_facets); auto & connectivities = mesh_accessor.getConnectivities(); /// this ones matter only for segments in 3D Array> * element_to_subfacet_double = nullptr; Array> * facet_to_subfacet_double = nullptr; if (dim == spatial_dimension - 2) { element_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 2, dim); facet_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 1, dim); } auto & element_to_subelement = mesh_facets.getElementToSubelement(); std::vector middle_nodes; for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double[0]; - if (old_facet.type != _segment_3) + if (old_facet.type != _segment_3) { continue; + } auto node = connectivities(old_facet, 2); - if (not mesh.isPureGhostNode(node)) + if (not mesh.isPureGhostNode(node)) { middle_nodes.push_back(node); + } } auto n = doubled_nodes.size(); doubleNodes(middle_nodes); for (auto && data : enumerate(make_view(facets_to_double, 2))) { auto facet = std::get<0>(data); auto & old_facet = std::get<1>(data)[0]; - if (old_facet.type != _segment_3) + if (old_facet.type != _segment_3) { continue; + } auto old_node = connectivities(old_facet, 2); - if (mesh.isPureGhostNode(old_node)) + if (mesh.isPureGhostNode(old_node)) { continue; + } auto new_node = doubled_nodes(n, 1); auto & new_facet = std::get<1>(data)[1]; connectivities(new_facet, 2) = new_node; if (dim == spatial_dimension - 2) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_subelement(new_facet, 0)); updateElementalConnectivity(mesh, old_node, new_node, (*element_to_subfacet_double)(facet), &(*facet_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_subelement(new_facet, 0)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserterHelper::insertCohesiveElement() { auto nb_new_facets = insertFacetsOnly(); - if (nb_new_facets == 0) + if (nb_new_facets == 0) { return 0; + } updateCohesiveData(); return nb_new_facets; } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserterHelper::insertFacetsOnly() { UInt spatial_dimension = mesh.getSpatialDimension(); - if (facets_to_double_by_dim[spatial_dimension - 1]->size() == 0) + if (facets_to_double_by_dim[spatial_dimension - 1]->empty()) { return 0; + } if (spatial_dimension == 1) { doublePointFacet(); } else if (spatial_dimension == 2) { doubleFacets<1>(); findSubfacetToDouble<1>(); doubleSubfacet<2>(); } else if (spatial_dimension == 3) { doubleFacets<2>(); findSubfacetToDouble<2>(); doubleFacets<1>(); findSubfacetToDouble<1>(); doubleSubfacet<3>(); } return facets_to_double_by_dim[spatial_dimension - 1]->size(); } /* -------------------------------------------------------------------------- */ template void CohesiveElementInserterHelper::doubleFacets() { AKANTU_DEBUG_IN(); NewElementsEvent new_facets; auto spatial_dimension = mesh_facets.getSpatialDimension(); auto & facets_to_double = *facets_to_double_by_dim[dim]; MeshAccessor mesh_accessor(mesh_facets); for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto && old_facet = facet_to_double[0]; auto && new_facet = facet_to_double[1]; auto & facets_connectivities = mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type); facets_connectivities.resize( nb_new_facets(old_facet.type, old_facet.ghost_type)); auto facets_connectivities_begin = make_view(facets_connectivities, facets_connectivities.getNbComponent()) .begin(); // copy the connectivities Vector new_conn(facets_connectivities_begin[new_facet.element]); Vector old_conn(facets_connectivities_begin[old_facet.element]); new_conn = old_conn; // this will fail if multiple facet types exists for a given element type // \TODO handle multiple sub-facet types auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(old_facet.type); auto & subfacets_to_facets = mesh_accessor.getSubelementToElementNC( old_facet.type, old_facet.ghost_type); subfacets_to_facets.resize( nb_new_facets(old_facet.type, old_facet.ghost_type), ElementNull); auto subfacets_to_facets_begin = make_view(subfacets_to_facets, nb_subfacet_per_facet).begin(); // copy the subfacet to facets information Vector old_subfacets_to_facet( subfacets_to_facets_begin[old_facet.element]); Vector new_subfacet_to_facet( subfacets_to_facets_begin[new_facet.element]); new_subfacet_to_facet = old_subfacets_to_facet; for (auto & subfacet : old_subfacets_to_facet) { - if (subfacet == ElementNull) + if (subfacet == ElementNull) { continue; + } /// update facet_to_subfacet array mesh_accessor.getElementToSubelement(subfacet).push_back(new_facet); } new_facets.getList().push_back(new_facet); } /// update facet_to_subfacet and _segment_3 facets if any - if (not(dim == spatial_dimension - 1)) { - updateSubelementToElement(dim - 1, true); - updateElementToSubelement(dim - 1, true); + if (dim != spatial_dimension - 1) { + updateSubelementToElement(dim, true); + updateElementToSubelement(dim, true); updateQuadraticSegments(dim); } else { updateQuadraticSegments(dim); } mesh_facets.sendEvent(new_facets); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void CohesiveElementInserterHelper::findSubfacetToDouble() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); MeshAccessor mesh_accessor(mesh_facets); auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; auto & subfacets_to_facets = mesh_facets.getSubelementToElement(); auto & elements_to_facets = mesh_accessor.getElementToSubelement(); /// loop on every facet for (auto f : arange(2)) { for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto old_facet = facet_to_double[f]; auto new_facet = facet_to_double[1 - f]; const auto & starting_element = elements_to_facets(new_facet, 0)[0]; auto current_facet = old_facet; Vector subfacets_to_facet = subfacets_to_facets.get(old_facet); /// loop on every subfacet for (auto & subfacet : subfacets_to_facet) { - if (subfacet == ElementNull) + if (subfacet == ElementNull) { continue; + } if (dim == spatial_dimension - 2) { Vector subsubfacets_to_subfacet = subfacets_to_facets.get(subfacet); /// loop on every subsubfacet for (auto & subsubfacet : subsubfacets_to_subfacet) { - if (subsubfacet == ElementNull) + if (subsubfacet == ElementNull) { continue; + } Vector subsubfacet_connectivity = mesh_facets.getConnectivity(subsubfacet); std::vector element_list; std::vector facet_list; std::vector subfacet_list; /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet( starting_element, current_facet, subsubfacet_connectivity, - element_list, facet_list, &subfacet_list) == false && + element_list, facet_list, &subfacet_list) == false and removeElementsInVector( subfacet_list, elements_to_facets(subsubfacet)) == false) { Element new_subsubfacet{ subsubfacet.type, nb_new_facets(subsubfacet.type, subsubfacet.ghost_type)++, subsubfacet.ghost_type}; facets_to_double_by_dim[dim - 1]->push_back( Vector{subsubfacet, new_subsubfacet}); - elementsOfDimToElementsOfDim(dim, dim - 1) + elementsOfDimToElementsOfDim(dim - 1, dim - 1) .push_back(subfacet_list); - elementsOfDimToElementsOfDim(dim + 1, dim - 1) + elementsOfDimToElementsOfDim(dim, dim - 1) .push_back(facet_list); - elementsOfDimToElementsOfDim(dim + 2, dim - 1) + elementsOfDimToElementsOfDim(dim + 1, dim - 1) .push_back(element_list); } } } else { std::vector element_list; std::vector facet_list; Vector subfacet_connectivity = mesh_facets.getConnectivity(subfacet); /// check if subfacet is to be doubled if (findElementsAroundSubfacet(starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false and removeElementsInVector(facet_list, elements_to_facets(subfacet)) == false) { Element new_subfacet{ subfacet.type, nb_new_facets(subfacet.type, subfacet.ghost_type)++, subfacet.ghost_type}; facets_to_double_by_dim[dim - 1]->push_back( Vector{subfacet, new_subfacet}); elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list); elementsOfDimToElementsOfDim(dim + 1, dim - 1) .push_back(element_list); } } } } } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::doubleNodes( const std::vector & old_nodes) { auto & position = mesh.getNodes(); auto spatial_dimension = mesh.getSpatialDimension(); auto old_nb_nodes = position.size(); position.reserve(old_nb_nodes + old_nodes.size()); doubled_nodes.reserve(doubled_nodes.size() + old_nodes.size()); auto position_begin = position.begin(spatial_dimension); for (auto && data : enumerate(old_nodes)) { auto n = std::get<0>(data); auto old_node = std::get<1>(data); decltype(old_node) new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes.push_back(Vector{old_node, new_node}); /// update position Vector coords = Vector(position_begin[old_node]); position.push_back(coords); } } /* -------------------------------------------------------------------------- */ bool CohesiveElementInserterHelper::findElementsAroundSubfacet( const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & element_list, std::vector & facet_list, std::vector * subfacet_list) { bool facet_matched = false; element_list.push_back(starting_element); std::queue elements_to_check; elements_to_check.push(starting_element); /// keep going as long as there are elements to check while (not elements_to_check.empty()) { /// check current element Element & current_element = elements_to_check.front(); const Vector facets_to_element = mesh_facets.getSubelementToElement(current_element); // for every facet of the element for (auto & current_facet : facets_to_element) { - if (current_facet == ElementNull) + if (current_facet == ElementNull) { continue; + } - if (current_facet == end_facet) + if (current_facet == end_facet) { facet_matched = true; + } auto find_facet = std::find(facet_list.begin(), facet_list.end(), current_facet); - // facet already listed - if (find_facet != facet_list.end()) - continue; - // subfacet_connectivity is not in the connectivity of current_facet; + // facet already listed or subfacet_connectivity is not in the + // connectivity of current_facet; if ((find_facet != facet_list.end()) or not hasElement(mesh_facets.getConnectivity(current_facet), - subfacet_connectivity)) + subfacet_connectivity)) { continue; + } facet_list.push_back(current_facet); if (subfacet_list != nullptr) { const Vector subfacets_of_facet = mesh_facets.getSubelementToElement(current_facet); /// check subfacets for (const auto & current_subfacet : subfacets_of_facet) { - if (current_subfacet == ElementNull) + if (current_subfacet == ElementNull) { continue; + } + + auto find_subfacet = std::find( + subfacet_list->begin(), subfacet_list->end(), current_subfacet); + if ((find_subfacet != subfacet_list->end()) or + not hasElement(mesh_facets.getConnectivity(current_subfacet), + subfacet_connectivity)) { + continue; + } - if ((std::find(subfacet_list->begin(), subfacet_list->end(), - current_subfacet) == subfacet_list->end()) and - hasElement(mesh_facets.getConnectivity(current_subfacet), - subfacet_connectivity)) - subfacet_list->push_back(current_subfacet); + subfacet_list->push_back(current_subfacet); } } /// consider opposing element const auto & elements_to_facet = mesh_facets.getElementToSubelement(current_facet); UInt opposing = 0; - if (elements_to_facet[0] == current_element) + if (elements_to_facet[0] == current_element) { opposing = 1; + } auto & opposing_element = elements_to_facet[opposing]; /// skip null elements since they are on a boundary - if (opposing_element == ElementNull) + if (opposing_element == ElementNull) { continue; + } /// skip this element if already added if (std::find(element_list.begin(), element_list.end(), - opposing_element) != element_list.end()) + opposing_element) != element_list.end()) { continue; + } /// only regular elements have to be checked - if (opposing_element.kind() == _ek_regular) + if (opposing_element.kind() == _ek_regular) { elements_to_check.push(opposing_element); + } element_list.push_back(opposing_element); AKANTU_DEBUG_ASSERT(hasElement(mesh.getConnectivity(opposing_element), subfacet_connectivity), "Subfacet doesn't belong to this element"); } /// erased checked element from the list elements_to_check.pop(); } return facet_matched; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::updateCohesiveData() { UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; MeshAccessor mesh_accessor(mesh); MeshAccessor mesh_facets_accessor(mesh_facets); for (auto ghost_type : ghost_types) { for (auto cohesive_type : mesh.elementTypes(_element_kind = _ek_cohesive)) { auto nb_cohesive_elements = mesh.getNbElement(cohesive_type, ghost_type); nb_new_facets(cohesive_type, ghost_type) = nb_cohesive_elements; mesh_facets_accessor.getSubelementToElement(cohesive_type, ghost_type); } } auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; new_elements.reserve(new_elements.size() + facets_to_double.size()); std::array facets; auto & element_to_facet = mesh_facets_accessor.getElementToSubelement(); auto & facets_to_cohesive_elements = mesh_facets_accessor.getSubelementToElement(); for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double[0]; /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension ? 1 : 0] = old_facet; facets[third_dimension ? 0 : 1] = facet_to_double[1]; auto type_cohesive = FEEngine::getCohesiveElementType(old_facet.type); auto & facet_connectivities = mesh_facets.getConnectivity(old_facet.type, old_facet.ghost_type); auto facet_connectivity_it = facet_connectivities.begin(facet_connectivities.getNbComponent()); auto cohesive_element = Element{ type_cohesive, nb_new_facets(type_cohesive, old_facet.ghost_type)++, old_facet.ghost_type}; auto & cohesives_connectivities = mesh_accessor.getConnectivity(type_cohesive, old_facet.ghost_type); Matrix connectivity(facet_connectivities.getNbComponent(), 2); Vector facets_to_cohesive_element(2); for (auto s : arange(2)) { /// store doubled facets facets_to_cohesive_element[s] = facets[s]; // update connectivities connectivity(s) = Vector(facet_connectivity_it[facets[s].element]); /// update element_to_facet vectors element_to_facet(facets[s], 0)[1] = cohesive_element; } cohesives_connectivities.push_back(connectivity); facets_to_cohesive_elements(type_cohesive, old_facet.ghost_type) .push_back(facets_to_cohesive_element); /// add cohesive element to the element event list new_elements.push_back(cohesive_element); } } /* -------------------------------------------------------------------------- */ void CohesiveElementInserterHelper::doublePointFacet() { UInt spatial_dimension = mesh.getSpatialDimension(); - if (spatial_dimension != 1) + if (spatial_dimension != 1) { return; + } auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1]; auto & element_to_facet = mesh_facets.getElementToSubelement(); auto & position = mesh.getNodes(); MeshAccessor mesh_accessor(mesh_facets); for (auto ghost_type : ghost_types) { for (auto facet_type : nb_new_facets.elementTypes( spatial_dimension - 1, ghost_type, _ek_regular)) { auto nb_new_element = nb_new_facets(facet_type, ghost_type); auto & connectivities = mesh_accessor.getConnectivity(facet_type, ghost_type); connectivities.resize(connectivities.size() + nb_new_element); } } for (auto facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double[0]; auto & new_facet = facet_to_double[1]; auto element = element_to_facet(new_facet)[0]; auto & facet_connectivities = mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type); auto old_node = facet_connectivities(old_facet.element); auto new_node = position.size(); /// update position position.push_back(position(old_node)); facet_connectivities(new_facet.element) = new_node; Vector segment_conectivity = mesh.getConnectivity(element); /// update facet connectivity auto it = std::find(segment_conectivity.begin(), segment_conectivity.end(), old_node); *it = new_node; doubled_nodes.push_back(Vector{old_node, new_node}); } } /* -------------------------------------------------------------------------- */ template void CohesiveElementInserterHelper::doubleSubfacet() { - if (spatial_dimension == 1) + if (spatial_dimension == 1) { return; + } std::vector nodes_to_double; MeshAccessor mesh_accessor(mesh_facets); auto & connectivities = mesh_accessor.getConnectivities(); auto & facets_to_double = *facets_to_double_by_dim[0]; ElementTypeMap nb_new_elements; for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto && old_element = facet_to_double[0]; nb_new_elements(old_element.type, old_element.ghost_type) = 0; } for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto && old_element = facet_to_double[0]; ++nb_new_elements(old_element.type, old_element.ghost_type); } for (auto ghost_type : ghost_types) { for (auto facet_type : nb_new_elements.elementTypes(0, ghost_type, _ek_regular)) { auto & connectivities = mesh_accessor.getConnectivity(facet_type, ghost_type); connectivities.resize(connectivities.size() + nb_new_elements(facet_type, ghost_type)); } } for (auto && facet_to_double : make_view(facets_to_double, 2)) { auto & old_facet = facet_to_double(0); // auto & new_facet = facet_to_double(1); AKANTU_DEBUG_ASSERT( old_facet.type == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); /// list nodes double nodes_to_double.push_back(connectivities(old_facet)); } auto old_nb_doubled_nodes = doubled_nodes.size(); doubleNodes(nodes_to_double); auto double_nodes_view = make_view(doubled_nodes, 2); for (auto && data : zip(make_view(facets_to_double, 2), range(double_nodes_view.begin() + old_nb_doubled_nodes, double_nodes_view.end()), arange(facets_to_double.size()))) { // auto & old_facet = std::get<0>(data)[0]; auto & new_facet = std::get<0>(data)[1]; auto & nodes = std::get<1>(data); auto old_node = nodes(0); auto new_node = nodes(1); auto f = std::get<2>(data); /// add new nodes in connectivity connectivities(new_facet) = new_node; updateElementalConnectivity(mesh, old_node, new_node, elementsOfDimToElementsOfDim(2, 0)(f), &elementsOfDimToElementsOfDim(1, 0)(f)); updateElementalConnectivity(mesh_facets, old_node, new_node, elementsOfDimToElementsOfDim(1, 0)(f)); - if (spatial_dimension == 3) + if (spatial_dimension == 3) { updateElementalConnectivity(mesh_facets, old_node, new_node, elementsOfDimToElementsOfDim(0, 0)(f)); + } } updateSubelementToElement(0, spatial_dimension == 2); updateElementToSubelement(0, spatial_dimension == 2); } } // namespace akantu diff --git a/src/mesh_utils/cohesive_element_inserter_helper.hh b/src/mesh_utils/cohesive_element_inserter_helper.hh index 28deddf96..dbc9ddf8f 100644 --- a/src/mesh_utils/cohesive_element_inserter_helper.hh +++ b/src/mesh_utils/cohesive_element_inserter_helper.hh @@ -1,115 +1,116 @@ /** * @file cohesive_element_inserter_helper.hh * * @author Nicolas Richart * * @date creation jeu sep 03 2020 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COHESIVE_ELEMENT_INSERTER_HELPER_HH__ #define __AKANTU_COHESIVE_ELEMENT_INSERTER_HELPER_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ class CohesiveElementInserterHelper { public: CohesiveElementInserterHelper( Mesh & mesh, const ElementTypeMapArray & facet_insertion); UInt insertCohesiveElement(); UInt insertFacetsOnly(); private: template UInt insertFacetsOnlyImpl(); template void doubleFacets(); template void findSubfacetToDouble(); void doubleNodes(const std::vector & old_nodes); bool findElementsAroundSubfacet( const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & element_list, std::vector & facet_list, std::vector * subfacet_list = nullptr); - inline bool hasElement(const Vector & nodes_element, - const Vector & nodes); - inline bool + static inline bool hasElement(const Vector & nodes_element, + const Vector & nodes); + static inline bool removeElementsInVector(const std::vector & elem_to_remove, std::vector & elem_list); void updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, const std::vector * facet_list = nullptr); // update functions void updateElementToSubelement(UInt dim, bool facet_mode); void updateSubelementToElement(UInt dim, bool facet_mode); void updateQuadraticSegments(UInt dim); void updateCohesiveData(); void doublePointFacet(); template void doubleSubfacet(); decltype(auto) elementsOfDimToElementsOfDim(Int dim1, Int dim2) { AKANTU_DEBUG_ASSERT(dim1 >= 0 and dim1 <= 3, "dimension of target element out of range"); AKANTU_DEBUG_ASSERT(dim2 >= 0 and dim2 <= 3, "dimension of source element out of range"); auto & array = dimelements_to_dimelements[dim1][dim2]; - if (not array) + if (not array) { array = std::make_unique>>(); + } return (*array); } public: - decltype(auto) getNewElements() const { return (new_elements); } - decltype(auto) getDoubledNodes() const { return (doubled_nodes); } + decltype(auto) getNewElements() const { return (new_elements); } + decltype(auto) getDoubledNodes() const { return (doubled_nodes); } private: std::array>, 3> facets_to_double_by_dim; - std::array>>, 2>, 3> + std::array>>, 2>, 4> dimelements_to_dimelements; Array doubled_nodes; Array new_elements; Mesh & mesh; Mesh & mesh_facets; ElementTypeMap nb_new_facets; }; } // namespace akantu #endif /* __AKANTU_COHESIVE_ELEMENT_INSERTER_HELPER_HH__ */ diff --git a/src/mesh_utils/cohesive_element_inserter_inline_impl.hh b/src/mesh_utils/cohesive_element_inserter_inline_impl.hh index 650fe9581..59540edb7 100644 --- a/src/mesh_utils/cohesive_element_inserter_inline_impl.hh +++ b/src/mesh_utils/cohesive_element_inserter_inline_impl.hh @@ -1,88 +1,88 @@ /** * @file cohesive_element_inserter_inline_impl.hh * * @author Marco Vocialta * * @date creation: Wed Nov 05 2014 * @date last modification: Fri Dec 08 2017 * * @brief Cohesive element inserter inline functions * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "cohesive_element_inserter.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_COHESIVE_ELEMENT_INSERTER_INLINE_IMPL_HH_ #define AKANTU_COHESIVE_ELEMENT_INSERTER_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt CohesiveElementInserter::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if (tag == SynchronizationTag::_ce_groups) { size = elements.size() * sizeof(bool); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void CohesiveElementInserter::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); if (tag == SynchronizationTag::_ce_groups) { for (const auto & el : elements) { - const bool & data = insertion_facets(el, 0); + const bool & data = insertion_facets(el); buffer << data; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void CohesiveElementInserter::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (tag == SynchronizationTag::_ce_groups) { for (const auto & el : elements) { - bool & data = insertion_facets(el, 0); + bool & data = insertion_facets(el); buffer >> data; } } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* AKANTU_COHESIVE_ELEMENT_INSERTER_INLINE_IMPL_HH_ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index 4dca799d0..cf6d9c732 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,809 +1,810 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Marco Vocialta * * @date creation: Fri Aug 20 2010 * @date last modification: Wed Feb 21 2018 * * @brief All mesh utils necessary for various tasks * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) { spatial_dimension = mesh.getSpatialDimension(); } /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); for_each_element( mesh, [&](auto && element) { Vector conn = mesh.getConnectivity(element); for (auto && node : conn) { ++node_to_elem.rowOffset(node); } }, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list // Element e; node_to_elem.beginInsertions(); for_each_element( mesh, [&](auto && element) { Vector conn = mesh.getConnectivity(element); for (auto && node : conn) { node_to_elem.insertInRow(node, element); } }, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).size(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { ++node_to_elem.rowOffset(conn_val[el_offset + n]); } } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (auto ghost_type : ghost_types) { for (const auto & type : mesh.elementTypes(spatial_dimension - 1, ghost_type)) { mesh.getConnectivity(type, ghost_type).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension) { AKANTU_DEBUG_IN(); to_dimension = std::max(to_dimension, UInt(0)); AKANTU_DEBUG_ASSERT( mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension); /// sort facets and generate sub-facets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension) { AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the // sub-elements // example: if the function is called with mesh = mesh_facets const Mesh * mesh_facets_parent = nullptr; try { mesh_facets_parent = &mesh_facets.getMeshParent(); } catch (...) { } mesh_facets.defineMeshParent(mesh); MeshAccessor mesh_accessor(mesh_facets); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const auto mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; std::vector connected_elements; NewElementsEvent event(AKANTU_CURRENT_FUNCTION); // init the SubelementToElement data to improve performance for (auto && ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(dimension, ghost_type)) { auto & subelement_to_element = mesh_accessor.getSubelementToElement(type, ghost_type); subelement_to_element.resize(mesh.getNbElement(type, ghost_type), ElementNull); auto facet_types = mesh.getAllFacetTypes(type); for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); mesh_accessor.getElementToSubelement(facet_type, ghost_type); mesh_accessor.getConnectivity(facet_type, ghost_type); } } } const ElementSynchronizer * synchronizer = nullptr; if (mesh.isDistributed()) { synchronizer = &(mesh.getElementSynchronizer()); } Element current_element; for (auto && ghost_type : ghost_types) { GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; for (auto && type : mesh.elementTypes(dimension, ghost_type)) { auto facet_types = mesh.getAllFacetTypes(type); current_element.type = type; for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); auto nb_element = mesh.getNbElement(type, ghost_type); auto && element_to_subelement = &mesh_accessor.getElementToSubelementNC(facet_type, ghost_type); auto && connectivity_facets = &mesh_accessor.getConnectivity(facet_type, ghost_type); auto nb_nodes_per_facet = connectivity_facets->getNbComponent(); // Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; auto && facets = mesh.getFacetConnectivity(current_element, ft).transpose(); for (auto facet : facets) { // facet = facets(f); UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.zero(); // loop over the other nodes to search intersecting elements, // which are the elements that share another node with the // starting element after first_node for (auto && data : enumerate(node_to_elem.getRow(facet(0)))) { auto && local_el = std::get<0>(data); auto && first_node = std::get<1>(data); for (auto n : arange(1, nb_nodes_per_facet)) { auto && node_elements = node_to_elem.getRow(facet(n)); counter(local_el) += std::count( node_elements.begin(), node_elements.end(), first_node); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; connected_elements.clear(); for (auto && data : enumerate(node_to_elem.getRow(facet(0)))) { if (not(counter(std::get<0>(data)) == nb_nodes_per_facet - 1)) { continue; } auto && real_el = std::get<1>(data); ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } if (minimum_el != current_element) { continue; } bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) { ++n; } if (n == nb_nodes_per_facet) { full_ghost_facet = true; } if (full_ghost_facet) { continue; } if (boundary_only and nb_element_connected_to_facet != 1) { continue; } std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.reserve(elements.size() + connected_elements.size()); for (auto && connected_element : connected_elements) { elements.push_back(connected_element); } if (nb_element_connected_to_facet == 1) { /// boundary facet elements.push_back(ElementNull); } else if (nb_element_connected_to_facet == 2) { /// internal facet /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = {_not_ghost, _not_ghost}; for (UInt el = 0; el < connected_elements.size(); ++el) { gt[el] = connected_elements[el].ghost_type; } if ((gt[0] == _not_ghost) xor (gt[1] == _not_ghost)) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { prank[el] = synchronizer->getRank(connected_elements[el]); } // ugly trick from Marco detected :P bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) { facet_ghost_type = _not_ghost; } else { facet_ghost_type = _ghost; } connectivity_facets = &mesh_accessor.getConnectivity( facet_type, facet_ghost_type); element_to_subelement = &mesh_accessor.getElementToSubelementNC( facet_type, facet_ghost_type); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->size() - 1; Element facet_element{facet_type, current_facet, facet_ghost_type}; event.getList().push_back(facet_element); /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (auto & loc_el : elements) { if (loc_el == ElementNull) { continue; } auto && subelements = mesh_accessor.getSubelementToElement(loc_el); for (auto & el : subelements) { if (el != ElementNull) { continue; } el = facet_element; break; } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = &mesh_accessor.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_accessor.getElementToSubelement( facet_type, facet_ghost_type); } } } } } } mesh_facets.sendEvent(event); // restore the parent of mesh_facet if (mesh_facets_parent != nullptr) { mesh_facets.defineMeshParent(*mesh_facets_parent); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, Array & local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.size(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element) * nb_nodes_per_element, renumbering_map); old_nodes_numbers.resize(renumbering_map.size()); for (auto & renumber_pair : renumbering_map) { old_nodes_numbers(renumber_pair.second) = renumber_pair.first; } renumbering_map.clear(); MeshAccessor mesh_accessor(mesh); /// copy the renumbered connectivity to the right place auto & local_conn = mesh_accessor.getConnectivity(type); local_conn.resize(nb_local_element); if (nb_local_element > 0) { memcpy(local_conn.storage(), local_connectivities.storage(), nb_local_element * nb_nodes_per_element * sizeof(UInt)); } auto & ghost_conn = mesh_accessor.getConnectivity(type, _ghost); ghost_conn.resize(nb_ghost_element); if (nb_ghost_element > 0) { std::memcpy(ghost_conn.storage(), local_connectivities.storage() + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); } auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost); ghost_counter.resize(nb_ghost_element, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity( Array & list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes.storage(); UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; auto it = renumbering_map.find(node); if (it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh, AKANTU_CURRENT_FUNCTION); Array & nodes_removed = remove_nodes.getList(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(_all_dimensions, ghost_type, _ek_not_defined)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); Array & connectivity = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity.size()); renumberNodesInConnectivity( connectivity, nb_element * nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); for (auto && pair : renumbering_map) { new_numbering(std::get<0>(pair)) = std::get<1>(pair); } for (UInt i = 0; i < new_numbering.size(); ++i) { if (new_numbering(i) == UInt(-1)) { nodes_removed.push_back(i); } } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets( Mesh & mesh_facets, const ElementTypeMapArray & remote_global_connectivities, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray local_global_connectivities( "local_global_connectivity", mesh_facets.getID(), mesh_facets.getMemoryID()); local_global_connectivities.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _ghost_type = gt_facet, _with_nb_nodes_per_element = true, _with_nb_element = true); mesh_facets.getGlobalConnectivity(local_global_connectivities); MeshAccessor mesh_accessor(mesh_facets); /// loop on every facet for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & connectivity = mesh_accessor.getConnectivity(type_facet, gt_facet); auto & local_global_connectivity = local_global_connectivities(type_facet, gt_facet); const auto & remote_global_connectivity = remote_global_connectivities(type_facet, gt_facet); auto & element_per_facet = mesh_accessor.getElementToSubelementNC(type_facet, gt_facet); auto & subfacet_to_facet = mesh_accessor.getSubelementToElementNC(type_facet, gt_facet); auto nb_nodes_per_facet = connectivity.getNbComponent(); auto nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); for (auto && data : zip(make_view(connectivity, nb_nodes_per_facet), make_view(local_global_connectivity, nb_nodes_per_facet), make_view(remote_global_connectivity, nb_nodes_per_facet), make_view(subfacet_to_facet, subfacet_to_facet.getNbComponent()), make_view(element_per_facet))) { auto & conn = std::get<0>(data); auto & local_gconn = std::get<1>(data); const auto & remote_gconn = std::get<2>(data); /// skip facet if connectivities are the same if (local_gconn == remote_gconn) { continue; } /// re-arrange connectivity auto conn_tmp = conn; auto begin = local_gconn.begin(); auto end = local_gconn.end(); AKANTU_DEBUG_ASSERT(std::is_permutation(begin, end, remote_gconn.begin()), "This facets are not just permutation of each other, " << local_gconn << " and " << remote_gconn); for (auto && data : enumerate(remote_gconn)) { auto it = std::find(begin, end, std::get<1>(data)); AKANTU_DEBUG_ASSERT(it != end, "Node not found"); UInt new_position = it - begin; conn(new_position) = conn_tmp(std::get<0>(data)); ; } // std::transform(remote_gconn.begin(), remote_gconn.end(), conn.begin(), // [&](auto && gnode) { // auto it = std::find(begin, end, gnode); // AKANTU_DEBUG_ASSERT(it != end, "Node not found"); // return conn_tmp(it - begin); // }); /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { auto begin = remote_gconn.begin(); /// find first node auto it = std::find(begin, remote_gconn.end(), local_gconn(0)); UInt n; UInt start = it - begin; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (n = 1; n < nb_nodes_per_P1_facet && local_gconn(n) == remote_gconn((start + n) % nb_nodes_per_P1_facet); ++n) { ; } /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) { continue; } } /// update data to invert facet auto & element_per_facet = std::get<4>(data); if (element_per_facet[1] != ElementNull) { // by convention the first facet // cannot be a ElementNull std::swap(element_per_facet[0], element_per_facet[1]); } auto & subfacets_of_facet = std::get<3>(data); std::swap(subfacets_of_facet(0), subfacets_of_facet(1)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call " "the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID(), mesh.getMemoryID()); barycenters.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = _all_dimensions); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (const auto & type : mesh.elementTypes(_all_dimensions, ghost_type)) { element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & barycenters_arr = barycenters(type, ghost_type); barycenters_arr.resize(nb_element); auto bary = barycenters_arr.begin(spatial_dimension); auto bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { element.element = el; mesh.getBarycenter(element, *bary); } } } MeshAccessor mesh_accessor(mesh); for (Int sp(spatial_dimension); sp >= 1; --sp) { if (mesh.getNbElement(sp) == 0) { continue; } for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(sp, ghost_type)) { auto & subelement_to_element = mesh_accessor.getSubelementToElement(type, ghost_type); subelement_to_element.resize(mesh.getNbElement(type, ghost_type)); subelement_to_element.set(ElementNull); } for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { auto & element_to_subelement = mesh_accessor.getElementToSubelement(type, ghost_type); element_to_subelement.resize(mesh.getNbElement(type, ghost_type)); element_to_subelement.clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (auto ghost_type : ghost_types) { facet_element.ghost_type = ghost_type; for (const auto & type : mesh.elementTypes(sp - 1, ghost_type)) { facet_element.type = type; auto & element_to_subelement = - mesh_accessor.getElementToSubelementNC(type, ghost_type); + mesh_accessor.getElementToSubelement(type, ghost_type); const auto & connectivity = mesh.getConnectivity(type, ghost_type); + //element_to_subelement.resize(connectivity.size()); for (auto && data : enumerate( make_view(connectivity, mesh.getNbNodesPerElement(type)))) { const auto & facet = std::get<1>(data); facet_element.element = std::get<0>(data); std::map element_seen_counter; auto nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(type)); // count the number of node in common between the facet and the // other element connected to the nodes of the facet for (auto node : arange(nb_nodes_per_facet)) { for (auto & elem : nodes_to_elements.getRow(facet(node))) { auto cit = element_seen_counter.find(elem); if (cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } // check which are the connected elements std::vector connected_elements; for (auto && cit : element_seen_counter) { if (cit.second == nb_nodes_per_facet) { connected_elements.push_back(cit.first); } } // add the connected elements as sub-elements for (auto & connected_element : connected_elements) { element_to_subelement(facet_element.element) .push_back(connected_element); } // add the element as sub-element to the connected elements for (auto & connected_element : connected_elements) { Vector subelements_to_element = mesh.getSubelementToElement(connected_element); // find the position where to insert the element auto it = std::find(subelements_to_element.begin(), subelements_to_element.end(), ElementNull); AKANTU_DEBUG_ASSERT( it != subelements_to_element.end(), "The element " << connected_element << " seems to have too many facets!! (" << (it - subelements_to_element.begin()) << " < " << mesh.getNbFacetsPerElement(connected_element.type) << ")"); *it = facet_element; } } } } } AKANTU_DEBUG_OUT(); } } // namespace akantu