diff --git a/src/fem/boundary.cc b/src/fem/boundary.cc index 854408199..cae52ee6c 100644 --- a/src/fem/boundary.cc +++ b/src/fem/boundary.cc @@ -1,417 +1,457 @@ /** * @file boundary.cc * * @author Dana Christen + * @author David Kammer * * @date Wed Mar 06 09:30:00 2013 * * @brief Stores information relevent to the notion of domain boundary and surfaces. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include #include "boundary.hh" #include "mesh.hh" #include "aka_csr.hh" #include "mesh_utils.hh" #include "sub_boundary.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -Boundary::Boundary(const Mesh & mesh, const ID & id, const ID & parent_id, const MemoryID & mem_id) -:memory_id(mem_id), mesh(mesh) +Boundary::Boundary(const Mesh & mesh, + const ID & id, + const ID & parent_id, + const MemoryID & mem_id) +: memory_id(mem_id), mesh(mesh) { AKANTU_DEBUG_IN(); std::stringstream sstr; if(parent_id != "") sstr << parent_id << ":"; sstr << id; this->id = sstr.str(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Boundary::~Boundary() { - for(BoundaryList::iterator boundaries_iter = boundaries.begin(); boundaries_iter != boundaries.end(); boundaries_iter++) - { - delete boundaries_iter->second; + AKANTU_DEBUG_IN(); + + this->removeAllSubBoundaries(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Boundary::removeSubBoundary(const std::string & name) { + AKANTU_DEBUG_IN(); + + BoundaryList::iterator b_it = this->boundaries.find(name); + AKANTU_DEBUG_ASSERT(b_it != this->boundaries.end(), + "Cannot remove SubBoundary: " << name << + " because it does not exist in this Boundary."); + + // delete the subboundary object + delete b_it->second; + + // delete it from the boundaries (list) + this->boundaries.erase(b_it); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Boundary::removeAllSubBoundaries() { + AKANTU_DEBUG_IN(); + + while (this->boundaries.size() != 0) { + this->removeSubBoundary(this->boundaries.begin()->first); } + + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Boundary::BoundaryTypeSet Boundary::getBoundaryElementTypes() { + AKANTU_DEBUG_IN(); // find surface element type that exists in this mesh (from David's code) UInt dim = mesh.getSpatialDimension(); BoundaryTypeSet surface_types; + const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); - for (Mesh::ConnectivityTypeList::const_iterator it = type_list.begin(); it != type_list.end(); ++it) - { + + for (Mesh::ConnectivityTypeList::const_iterator it = type_list.begin(); + it != type_list.end(); + ++it) { ElementType surface_type = mesh.getFacetType(*it); - if (mesh.getSpatialDimension(*it) == dim && mesh.getNbElement(surface_type) != 0) - { + if (mesh.getSpatialDimension(*it) == dim + && mesh.getNbElement(surface_type) != 0) { surface_types.insert(surface_type); } } + + AKANTU_DEBUG_OUT(); return surface_types; } /* -------------------------------------------------------------------------- */ void Boundary::addElementAndNodesToBoundaryAlloc(const std::string & boundary_name, const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(elem_type); BoundaryList::iterator boundaries_iter = boundaries.find(boundary_name); const Array & connectivity = mesh.getConnectivity(elem_type, ghost_type); if(boundaries_iter == boundaries.end()) { SubBoundary * sub_b = new SubBoundary(boundary_name, std::string(id + ":sub_boundary:" + boundary_name), memory_id); boundaries_iter = boundaries.insert(boundaries_iter, std::pair(boundary_name, sub_b)); } boundaries_iter->second->addElement(elem_type, elem_id, ghost_type); for (UInt n(0); n < nb_nodes_per_element; n++) { boundaries_iter->second->addNode(connectivity(elem_id, n)); } } /* -------------------------------------------------------------------------- */ void Boundary::createBoundariesFromGeometry() { AKANTU_DEBUG_IN(); CSR node_to_elem; /// Get list of surface elements UInt spatial_dimension = mesh.getSpatialDimension(); // buildNode2Elements(mesh, node_offset, node_to_elem, spatial_dimension-1); MeshUtils::buildNode2Elements(mesh, node_to_elem, spatial_dimension-1); /// Find which types of elements have been linearized const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); ElementType lin_element_type[nb_types]; UInt nb_lin_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types+1]; ElementType type_p1; for(it = type_list.begin(); it != type_list.end(); ++it) { //std::cout <<"Main type: " << *it << std::endl; ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) { continue; } ElementType facet_type = mesh.getFacetType(type); //std::cout <<"Facet type: " << facet_type << std::endl; lin_element_type[nb_lin_types] = facet_type; nb_nodes_per_element[nb_lin_types] = Mesh::getNbNodesPerElement(facet_type); type_p1 = Mesh::getP1ElementType(facet_type); nb_nodes_per_element_p1[nb_lin_types] = Mesh::getNbNodesPerElement(type_p1); conn_val[nb_lin_types] = mesh.getConnectivity(facet_type, _not_ghost).values; nb_element[nb_lin_types] = mesh.getNbElement(facet_type, _not_ghost); nb_lin_types++; } for (UInt i = 1; i < nb_lin_types; ++i) { nb_element[i] += nb_element[i+1]; } for (UInt i = nb_lin_types; i > 0; --i) { nb_element[i] = nb_element[i-1]; } nb_element[0] = 0; /// Find close surfaces Array surface_value_id(1, nb_element[nb_lin_types], -1); Int * surf_val = surface_value_id.storage(); UInt nb_boundaries = 0; UInt nb_cecked_elements; UInt nb_elements_to_ceck; UInt * elements_to_ceck = new UInt [nb_element[nb_lin_types]]; memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt)); for (UInt lin_el = 0; lin_el < nb_element[nb_lin_types]; ++lin_el) { if(surf_val[lin_el] != -1) continue; /* Surface id already assigned */ /* First element of new surface */ surf_val[lin_el] = nb_boundaries; nb_cecked_elements = 0; nb_elements_to_ceck = 1; memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt)); elements_to_ceck[0] = lin_el; // Find others elements belonging to this surface while(nb_cecked_elements < nb_elements_to_ceck) { UInt ceck_lin_el = elements_to_ceck[nb_cecked_elements]; // Transform linearized index of element into ElementType one UInt lin_type_nb = 0; while (ceck_lin_el >= nb_element[lin_type_nb+1]) { lin_type_nb++; } UInt ceck_el = ceck_lin_el - nb_element[lin_type_nb]; // Get connected elements UInt el_offset = ceck_el*nb_nodes_per_element[lin_type_nb]; for (UInt n = 0; n < nb_nodes_per_element_p1[lin_type_nb]; ++n) { UInt node_id = conn_val[lin_type_nb][el_offset + n]; CSR::iterator it_n; for (it_n = node_to_elem.begin(node_id); it_n != node_to_elem.end(node_id); ++it_n) { if(surf_val[*it_n] == -1) { /* Found new surface element */ surf_val[*it_n] = nb_boundaries; elements_to_ceck[nb_elements_to_ceck] = *it_n; nb_elements_to_ceck++; } } } nb_cecked_elements++; } nb_boundaries++; } delete [] elements_to_ceck; /// Transform local linearized element index in the global one for (UInt i = 0; i < nb_lin_types; ++i) { nb_element[i] = nb_element[i+1] - nb_element[i]; } UInt el_offset = 0; for(UInt type_it = 0; type_it < nb_lin_types; ++type_it) { ElementType type = lin_element_type[type_it]; //std::cout << "Second type: " << type << std::endl; for (UInt el = 0; el < nb_element[type_it]; ++el) { std::stringstream sstr; sstr << surf_val[el+el_offset]; std::string surf_name(sstr.str()); //std::cout << "Inserting elements and nodes to surface " << surf_name << std::endl; addElementAndNodesToBoundaryAlloc(surf_name, type, el); } //std::cout << "Offset +" << nb_element[type_it] << std::endl; el_offset += nb_element[type_it]; } BoundaryList::iterator subB_iter = boundaries.begin(); BoundaryList::iterator subB_iter_end = boundaries.end(); for(;subB_iter != subB_iter_end; ++subB_iter) { SubBoundary & sub = *subB_iter->second; sub.cleanUpNodeList(); sub.registerDumper("paraview_"+sub.getName(), sub.getID(), true); sub.addDumpFilteredMesh(mesh, sub.elements, sub.nodes, mesh.getSpatialDimension() - 1, _not_ghost, _ek_regular); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Boundary::createBoundariesFromMeshData(const std::string & dataset_name) { UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator type_it = mesh.firstType(spatial_dimension - 1, *gt); Mesh::type_iterator type_end = mesh.lastType(spatial_dimension - 1, *gt); for (; type_it != type_end; ++type_it) { // Dirty check to assert that the processor has indeed this type of element // associated to the dataset name mesh.getData(*type_it, dataset_name, *gt); const Array & dataset = mesh.getData(*type_it, dataset_name, *gt); UInt nb_element = mesh.getNbElement(*type_it, *gt); AKANTU_DEBUG_ASSERT(dataset.getSize() == nb_element, "Not the same number of elements in the map from MeshData and in the mesh!"); // FIXME This could be improved (performance...) for(UInt i(0); i < nb_element; ++i) { std::stringstream sstr; sstr << dataset(i); std::string boundary_name = sstr.str(); addElementAndNodesToBoundaryAlloc(boundary_name, *type_it, i, *gt); } } } BoundaryList::iterator subB_iter = boundaries.begin(); BoundaryList::iterator subB_iter_end = boundaries.end(); for(;subB_iter != subB_iter_end; ++subB_iter) { SubBoundary & sub = *subB_iter->second; sub.cleanUpNodeList(); sub.registerDumper("paraview_"+sub.getName(), sub.getID(), true); sub.addDumpFilteredMesh(mesh, sub.elements, sub.nodes, mesh.getSpatialDimension() - 1, _not_ghost, _ek_regular); } } template void Boundary::createBoundariesFromMeshData(const std::string & dataset_name); template void Boundary::createBoundariesFromMeshData(const std::string & dataset_name); /* -------------------------------------------------------------------------- */ void Boundary::createSubBoundaryFromNodeGroup(const std::string & name, const Array & node_group) { CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem, mesh.getSpatialDimension() - 1); std::set seen; BoundaryList::iterator boundaries_iter = boundaries.find(name); if(boundaries_iter == boundaries.end()) { // Manual construction of the sub-boundary ID SubBoundary * sub_b = new SubBoundary(name, std::string(id + "_sub_boundary_" + name), memory_id); boundaries_iter = boundaries.insert(boundaries_iter, std::pair(name, sub_b)); } SubBoundary & sub_bound = *boundaries_iter->second; Array::const_iterator<> itn = node_group.begin(); Array::const_iterator<> endn = node_group.end(); for (;itn != endn; ++itn) { CSR::iterator ite = node_to_elem.begin(*itn); CSR::iterator ende = node_to_elem.end(*itn); for (;ite != ende; ++ite) { const Element & elem = *ite; if(seen.find(elem) != seen.end()) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(elem.type); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(elem.type, elem.ghost_type).begin(nb_nodes_per_element); const Vector & conn = conn_it[elem.element]; UInt count = 0; for (UInt n = 0; n < conn.size(); ++n) { count += (node_group.find(conn(n)) != -1 ? 1 : 0); } if(count == nb_nodes_per_element) { sub_bound.addElement(elem.type, elem.element, elem.ghost_type); for (UInt n(0); n < nb_nodes_per_element; n++) { sub_bound.addNode(conn(n)); } } seen.insert(elem); } } sub_bound.cleanUpNodeList(); sub_bound.registerDumper("paraview_"+sub_bound.getName(), sub_bound.getID(), true); sub_bound.addDumpFilteredMesh(mesh, sub_bound.elements, sub_bound.nodes, mesh.getSpatialDimension() - 1, _not_ghost, _ek_regular); } /* -------------------------------------------------------------------------- */ void Boundary::createBoundariesFromMeshData(const std::string & dataset_name) { createBoundariesFromMeshData(dataset_name); } /* -------------------------------------------------------------------------- */ void Boundary::printself(std::ostream & stream) const { for(const_iterator it(begin()); it != end(); ++it) { stream << "-- Boundary \"" << it->getName() << "\":" << std::endl; // Loop over the nodes stream << "---- " << it->getNbNodes() << " nodes"; if(AKANTU_DEBUG_TEST(dblDump)) { stream << ":" << std::endl; stream << "------ "; for(SubBoundary::nodes_const_iterator nodes_it(it->nodes_begin()); nodes_it!= it->nodes_end(); ++nodes_it) { stream << *nodes_it << " "; } } stream << std::endl; // Loop over the element types stream << "---- Elements:" << std::endl; Mesh::type_iterator type_it = mesh.firstType(mesh.getSpatialDimension()-1); Mesh::type_iterator type_end = mesh.lastType(mesh.getSpatialDimension()-1); for(; type_it != type_end; ++type_it) { stream << "------ Type " << *type_it; const Array & element_ids = it->getElements(*type_it); if(AKANTU_DEBUG_TEST(dblDump)) { stream << ":" << std::endl; stream << "-------- "; // Loop over the corresponding elements in the SubBoundary for(UInt i(0); igetElements(*type_it); } } stream << std::endl; } /* -------------------------------------------------------------------------- */ void Boundary::dump() { for(iterator it(begin()); it != end(); ++it) { it->dump(); } } __END_AKANTU__ diff --git a/src/fem/boundary.hh b/src/fem/boundary.hh index 845d55edc..bd3645d92 100644 --- a/src/fem/boundary.hh +++ b/src/fem/boundary.hh @@ -1,143 +1,155 @@ /** * @file boundary.hh * * @author Dana Christen + * @author David Kammer * * @date Wed Mar 06 09:30:00 2013 * * @brief Stores information relevent to the notion of domain boundary and surfaces. * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_BOUNDARY_HH__ #define __AKANTU_BOUNDARY_HH__ #include #include "aka_common.hh" #include "by_element_type.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ class Mesh; /* -------------------------------------------------------------------------- */ class SubBoundary; class Boundary { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: public: typedef std::map BoundaryList; typedef std::set BoundaryTypeSet; /* ------------------------------------------------------------------------ */ /* SubBoundary iterator */ /* ------------------------------------------------------------------------ */ public: template class internal_iterator { public: inline internal_iterator(const internal_iterator &); private: explicit inline internal_iterator(const container_iterator &); public: inline T & operator*() const; inline T * operator->() const; inline bool operator==(const internal_iterator &) const; inline bool operator!=(const internal_iterator &) const; inline internal_iterator & operator=(const internal_iterator &); inline internal_iterator & operator++(); inline internal_iterator operator++(int); private: friend class Boundary; container_iterator iter; }; typedef internal_iterator const_iterator; typedef internal_iterator iterator; inline const_iterator begin() const; inline const_iterator find(std::string name) const; inline const_iterator end() const; inline iterator begin(); inline iterator find(std::string name); inline iterator end(); /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - Boundary(const Mesh & mesh, const ID & id = "boundary", const ID & parent_id = "", const MemoryID & memory_id = 0); + Boundary(const Mesh & mesh, + const ID & id = "boundary", + const ID & parent_id = "", + const MemoryID & memory_id = 0); ~Boundary(); /* ------------------------------------------------------------------------ */ /* Methods and accessors */ /* ------------------------------------------------------------------------ */ public: - template void createBoundariesFromMeshData(const std::string & dataset_name); + template + void createBoundariesFromMeshData(const std::string & dataset_name); + void createBoundariesFromMeshData(const std::string & dataset_name); void createBoundariesFromGeometry(); /// Create a SubBoundary based on a node group void createSubBoundaryFromNodeGroup(const std::string & name, const Array & node_group); + inline const SubBoundary & operator()(const std::string & name) const; + BoundaryTypeSet getBoundaryElementTypes(); + + void removeSubBoundary(const std::string & name); + void removeAllSubBoundaries(); + AKANTU_GET_MACRO(NbBoundaries, boundaries.size(), UInt); void printself(std::ostream & stream) const; void dump(); private: void addElementAndNodesToBoundaryAlloc(const std::string & boundary_name, const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: ID id; BoundaryList boundaries; MemoryID memory_id; const Mesh & mesh; }; #include "boundary_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_BOUNDARY_HH__ */ diff --git a/src/mesh_utils/mesh_partition.cc b/src/mesh_utils/mesh_partition.cc index 199d06e61..54298749e 100644 --- a/src/mesh_utils/mesh_partition.cc +++ b/src/mesh_utils/mesh_partition.cc @@ -1,380 +1,417 @@ /** * @file mesh_partition.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date Tue Aug 17 16:19:43 2010 * * @brief implementation of common part of all partitioner * * @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 "mesh_partition.hh" #include "mesh_utils.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MeshPartition::MeshPartition(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id) : Memory(memory_id), id("MeshPartitioner"), mesh(mesh), spatial_dimension(spatial_dimension), partitions ("partition" , id, memory_id), ghost_partitions ("ghost_partition" , id, memory_id), - ghost_partitions_offset("ghost_partition_offset", id, memory_id) { + ghost_partitions_offset("ghost_partition_offset", id, memory_id), + saved_connectivity("saved_connectivity", id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MeshPartition::~MeshPartition() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * conversion in c++ of the GENDUALMETIS (mesh.c) function wrote by George in * Metis (University of Minnesota) */ void MeshPartition::buildDualGraph(Array & dxadj, Array & dadjncy, Array & edge_loads, - const EdgeLoadFunctor & edge_load_func, - const Array & pairs) { + const EdgeLoadFunctor & edge_load_func) { AKANTU_DEBUG_IN(); // tweak mesh; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); UInt nb_types = type_list.size(); UInt nb_good_types = 0; - UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt magic_number[nb_types]; // UInt * conn_val[nb_types]; UInt nb_element[nb_types]; Array * conn[nb_types]; - Array * conn_tmp[nb_types]; Array lin_to_element; Element elem; elem.ghost_type = _not_ghost; const_cast(mesh).updateTypesOffsets(_not_ghost); const_cast(mesh).updateTypesOffsets(_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != mesh.getSpatialDimension()) continue; elem.type = type; ElementType type_p1 = Mesh::getP1ElementType(type); - nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); nb_nodes_per_element_p1[nb_good_types] = Mesh::getNbNodesPerElement(type_p1); nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); magic_number[nb_good_types] = Mesh::getNbNodesPerElement(Mesh::getFacetType(type_p1)); conn[nb_good_types] = &const_cast &>(mesh.getConnectivity(type, _not_ghost)); for (UInt i = 0; i < nb_element[nb_good_types]; ++i) { elem.element = i; lin_to_element.push_back(elem); } - - if(pairs.getSize() != 0) { - conn_tmp[nb_good_types] = new Array(mesh.getConnectivity(type, _not_ghost)); - for (UInt i = 0; i < pairs.getSize(); ++i) { - for (UInt el = 0; el < nb_element[nb_good_types]; ++el) { - for (UInt n = 0; n < nb_nodes_per_element[nb_good_types]; ++n) { - if(pairs(i, 1) == (*conn[nb_good_types])(el, n)) - (*conn[nb_good_types])(el, n) = pairs(i, 0); - } - } - } - } - nb_good_types++; } CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); UInt nb_total_element = 0; UInt nb_total_node_element = 0; for (UInt t = 0; t < nb_good_types; ++t) { nb_total_element += nb_element[t]; nb_total_node_element += nb_element[t]*nb_nodes_per_element_p1[t]; } dxadj.resize(nb_total_element + 1); /// initialize the dxadj array for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) dxadj(linerized_el) = nb_nodes_per_element_p1[t]; /// convert the dxadj_val array in a csr one for (UInt i = 1; i < nb_total_element; ++i) dxadj(i) += dxadj(i-1); for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i-1); dxadj(0) = 0; dadjncy.resize(2*dxadj(nb_total_element)); edge_loads.resize(2*dxadj(nb_total_element)); /// weight map to determine adjacency unordered_map::type weight_map; for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { /// fill the weight map for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) { UInt node = (*conn[t])(el, n); CSR::iterator k; for (k = node_to_elem.rbegin(node); k != node_to_elem.rend(node); --k) { UInt current_el = *k; if(current_el <= linerized_el) break; unordered_map::type::iterator it_w; it_w = weight_map.find(current_el); if(it_w == weight_map.end()) { weight_map[current_el] = 1; } else { it_w->second++; } } } /// each element with a weight of the size of a facet are adjacent unordered_map::type::iterator it_w; for(it_w = weight_map.begin(); it_w != weight_map.end(); ++it_w) { if(it_w->second == magic_number[t]) { UInt adjacent_el = it_w->first; #if defined(AKANTU_COHESIVE_ELEMENT) /// Patch in order to prevent neighboring cohesive elements /// from detecting each other ElementKind linearized_el_kind = mesh.linearizedToElement(linerized_el).kind; ElementKind adjacent_el_kind = mesh.linearizedToElement(adjacent_el).kind; if (linearized_el_kind == adjacent_el_kind && linearized_el_kind == _ek_cohesive) continue; #endif UInt index_adj = dxadj(adjacent_el )++; UInt index_lin = dxadj(linerized_el)++; dadjncy(index_lin) = adjacent_el; dadjncy(index_adj) = linerized_el; } } weight_map.clear(); } } - if(pairs.getSize() != 0) { - for (UInt i = 0; i < nb_good_types; ++i) { - conn[i]->copy(*conn_tmp[i]); - delete conn_tmp[i]; - } - } - Int k_start = 0; for (UInt t = 0, linerized_el = 0, j = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { for (Int k = k_start; k < dxadj(linerized_el); ++k, ++j) dadjncy(j) = dadjncy(k); dxadj(linerized_el) = j; k_start += nb_nodes_per_element_p1[t]; } for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i - 1); dxadj(0) = 0; UInt adj = 0; for (UInt i = 0; i < nb_total_element; ++i) { UInt nb_adj = dxadj(i + 1) - dxadj(i); for (UInt j = 0; j < nb_adj; ++j, ++adj) { Int el_adj_id = dadjncy(dxadj(i) + j); Element el = lin_to_element(i); Element el_adj = lin_to_element(el_adj_id); Int load = edge_load_func(el, el_adj); edge_loads(adj) = load; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::fillPartitionInformation(const Mesh & mesh, const Int * linearized_partitions) { AKANTU_DEBUG_IN(); CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined); UInt linearized_el = 0; for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); partitions .alloc(nb_element, 1, type, _not_ghost); CSR & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions .alloc(0, 1, type, _ghost); const Array & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt part = linearized_partitions[linearized_el]; partitions(type, _not_ghost)(el) = part; std::list list_adj_part; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity.values[el * nb_nodes_per_element + n]; CSR::iterator ne; for (ne = node_to_elem.begin(node); ne != node_to_elem.end(node); ++ne) { UInt adj_el = *ne; UInt adj_part = linearized_partitions[adj_el]; if(part != adj_part) { list_adj_part.push_back(adj_part); } } } list_adj_part.sort(); list_adj_part.unique(); for(std::list::iterator adj_it = list_adj_part.begin(); adj_it != list_adj_part.end(); ++adj_it) { ghost_part_csr.getRows().push_back(*adj_it); ghost_part_csr.rowOffset(el)++; ghost_partitions(type, _ghost).push_back(*adj_it); ghost_partitions_offset(type, _ghost)(el)++; } } ghost_part_csr.countToCSR(); /// convert the ghost_partitions_offset array in an offset array Array & ghost_partitions_offset_ptr = ghost_partitions_offset(type, _ghost); for (UInt i = 1; i < nb_element; ++i) ghost_partitions_offset_ptr(i) += ghost_partitions_offset_ptr(i-1); for (UInt i = nb_element; i > 0; --i) ghost_partitions_offset_ptr(i) = ghost_partitions_offset_ptr(i-1); ghost_partitions_offset_ptr(0) = 0; } // Facets Mesh::type_iterator fit = mesh.firstType(spatial_dimension - 1, _not_ghost, _ek_not_defined); Mesh::type_iterator fend = mesh.lastType(spatial_dimension - 1, _not_ghost, _ek_not_defined); for(; fit != fend; ++fit) { ElementType type = *fit; UInt nb_element = mesh.getNbElement(type); partitions .alloc(nb_element, 1, type, _not_ghost); AKANTU_DEBUG_WARNING("Allocating partitions for " << type); CSR & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions .alloc(0, 1, type, _ghost); AKANTU_DEBUG_WARNING("Allocating ghost_partitions for " << type); const Array< std::vector > & elem_to_subelem = mesh.getElementToSubelement(type, _not_ghost); for(UInt i(0); i < mesh.getNbElement(type, _not_ghost); ++i) { // Facet loop const std::vector & adjacent_elems = elem_to_subelem(i); Element min_elem; UInt min_part(std::numeric_limits::max()); std::set adjacent_parts; for(UInt j(0); j < adjacent_elems.size(); ++j) { UInt adjacent_elem_id = adjacent_elems[j].element; UInt adjacent_elem_part = partitions(adjacent_elems[j].type, adjacent_elems[j].ghost_type)(adjacent_elem_id); if(adjacent_elem_part < min_part) { min_part = adjacent_elem_part; min_elem = adjacent_elems[j]; } adjacent_parts.insert(adjacent_elem_part); } partitions(type, _not_ghost)(i) = min_part; CSR::iterator git = ghost_partitions_csr(min_elem.type, _not_ghost).begin(min_elem.element); CSR::iterator gend = ghost_partitions_csr(min_elem.type, _not_ghost).end(min_elem.element); for(; git != gend; ++git) { adjacent_parts.insert(*git); } adjacent_parts.erase(min_part); std::set::const_iterator pit = adjacent_parts.begin(); std::set::const_iterator pend = adjacent_parts.end(); for(; pit != pend; ++pit) { ghost_part_csr.getRows().push_back(*pit); ghost_part_csr.rowOffset(i)++; ghost_partitions(type, _ghost).push_back(*pit); } ghost_partitions_offset(type, _ghost)(i+1) = ghost_partitions_offset(type, _ghost)(i+1) + adjacent_elems.size(); } ghost_part_csr.countToCSR(); } AKANTU_DEBUG_OUT(); } + +/* -------------------------------------------------------------------------- */ +void MeshPartition::tweakConnectivity(const Array & pairs) { + AKANTU_DEBUG_IN(); + + if(pairs.getSize() == 0) return; + + Mesh::type_iterator it = mesh.firstType(spatial_dimension, + _not_ghost, + _ek_not_defined); + Mesh::type_iterator end = mesh.lastType(spatial_dimension, + _not_ghost, + _ek_not_defined); + + for(; it != end; ++it) { + ElementType type = *it; + + Array & conn = const_cast &>(mesh.getConnectivity(type, _not_ghost)); + UInt nb_nodes_per_element = conn.getNbComponent(); + UInt nb_element = conn.getSize(); + + Array & saved_conn = saved_connectivity.alloc(nb_element, nb_nodes_per_element, type, _not_ghost); + saved_conn.copy(conn); + + for (UInt i = 0; i < pairs.getSize(); ++i) { + for (UInt el = 0; el < nb_element; ++el) { + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + if(pairs(i, 1) == conn(el, n)) + conn(el, n) = pairs(i, 0); + } + } + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MeshPartition::restoreConnectivity() { + AKANTU_DEBUG_IN(); + + ByElementTypeUInt::type_iterator it = saved_connectivity.firstType(spatial_dimension, + _not_ghost, + _ek_not_defined); + ByElementTypeUInt::type_iterator end = saved_connectivity.lastType(spatial_dimension, + _not_ghost, + _ek_not_defined); + for(; it != end; ++it) { + ElementType type = *it; + + Array & conn = const_cast &>(mesh.getConnectivity(type, _not_ghost)); + Array & saved_conn = saved_connectivity(type, _not_ghost); + conn.copy(saved_conn); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + __END_AKANTU__ diff --git a/src/mesh_utils/mesh_partition.hh b/src/mesh_utils/mesh_partition.hh index 33ae76eb5..c2b7f6886 100644 --- a/src/mesh_utils/mesh_partition.hh +++ b/src/mesh_utils/mesh_partition.hh @@ -1,157 +1,164 @@ /** * @file mesh_partition.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date Mon Aug 16 13:17:22 2010 * * @brief tools to partitionate a mesh * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_PARTITION_HH__ #define __AKANTU_MESH_PARTITION_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "aka_csr.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class MeshPartition : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshPartition(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id = 0); virtual ~MeshPartition(); class EdgeLoadFunctor { public: virtual Int operator()(__attribute__((unused)) const Element & el1, __attribute__((unused)) const Element & el2) const = 0; }; class ConstEdgeLoadFunctor : public EdgeLoadFunctor { public: virtual inline Int operator()(__attribute__((unused)) const Element & el1, __attribute__((unused)) const Element & el2) const { return 1; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void partitionate(UInt nb_part, const EdgeLoadFunctor & edge_load_func = ConstEdgeLoadFunctor(), const Array & pairs = Array()) = 0; virtual void reorder() = 0; /// function to print the contain of the class //virtual void printself(std::ostream & stream, int indent = 0) const = 0; /// fill the partitions array with a given linearized partition information void fillPartitionInformation(const Mesh & mesh, const Int * linearized_partitions); protected: /// build the dual graph of the mesh, for all element of spatial_dimension void buildDualGraph(Array & dxadj, Array & dadjncy, Array & edge_loads, - const EdgeLoadFunctor & edge_load_func, - const Array & pairs); + const EdgeLoadFunctor & edge_load_func); + + /// tweak the mesh to handle the PBC pairs + void tweakConnectivity(const Array & pairs); + /// restore the mesh that has been tweaked + void restoreConnectivity(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Partitions, partitions, const ByElementTypeUInt &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Partition, partitions, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GhostPartition, ghost_partitions, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GhostPartitionOffset, ghost_partitions_offset, UInt); AKANTU_GET_MACRO(GhostPartitionCSR, ghost_partitions_csr, const ByElementType< CSR > &); AKANTU_GET_MACRO(NbPartition, nb_partitions, UInt); AKANTU_SET_MACRO(NbPartition, nb_partitions, UInt); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id std::string id; /// the mesh to partition const Mesh & mesh; /// dimension of the elements to consider in the mesh UInt spatial_dimension; /// number of partitions UInt nb_partitions; /// partition numbers ByElementTypeUInt partitions; ByElementType< CSR > ghost_partitions_csr; ByElementTypeUInt ghost_partitions; ByElementTypeUInt ghost_partitions_offset; Array * permutation; + + ByElementTypeUInt saved_connectivity; + }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "mesh_partition_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const MeshPartition & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #ifdef AKANTU_USE_SCOTCH #include "mesh_partition_scotch.hh" #endif #endif /* __AKANTU_MESH_PARTITION_HH__ */ diff --git a/src/mesh_utils/mesh_partition/mesh_partition_mesh_data.cc b/src/mesh_utils/mesh_partition/mesh_partition_mesh_data.cc index 5aba38af2..25182a340 100644 --- a/src/mesh_utils/mesh_partition/mesh_partition_mesh_data.cc +++ b/src/mesh_utils/mesh_partition/mesh_partition_mesh_data.cc @@ -1,110 +1,119 @@ /** * @file mesh_partition_mesh_data.cc * * @author Dana Christen * * @date Mon May 01 10:22:00 2013 (On Labor Day, such a shame!) * * @brief implementation of the MeshPartitionMeshData class * * @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 "mesh_partition_mesh_data.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MeshPartitionMeshData::MeshPartitionMeshData(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id) : MeshPartition(mesh, spatial_dimension, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MeshPartitionMeshData::MeshPartitionMeshData(const Mesh & mesh, const ByElementTypeArray & mapping, UInt spatial_dimension, const MemoryID & memory_id) : MeshPartition(mesh, spatial_dimension, memory_id), partition_mapping(&mapping) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartitionMeshData::partitionate(UInt nb_part, const EdgeLoadFunctor & edge_load_func, const Array & pairs) { AKANTU_DEBUG_IN(); + + tweakConnectivity(pairs); + nb_partitions = nb_part; GhostType ghost_type = _not_ghost; UInt spatial_dimension = mesh.getSpatialDimension(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type, _ek_not_defined); UInt linearized_el = 0; UInt nb_elements = mesh.getNbElement(mesh.getSpatialDimension(), ghost_type); Int * partition_list = new Int[nb_elements]; for(; it != end; ++it) { ElementType type = *it; const Array & partition_array = (*partition_mapping)(type, ghost_type); Array::const_iterator< Vector > p_it = partition_array.begin(1); Array::const_iterator< Vector > p_end = partition_array.end(1); - AKANTU_DEBUG_ASSERT(p_end-p_it == mesh.getNbElement(type, ghost_type), "The partition mapping does not have the right number of entries for type " << type << " and ghost type " << ghost_type << "."); + AKANTU_DEBUG_ASSERT(p_end-p_it == mesh.getNbElement(type, ghost_type), + "The partition mapping does not have the right number " + << "of entries for type " << type << " and ghost type " + << ghost_type << "." << " Tags=" << p_end-p_it + << " Mesh=" << mesh.getNbElement(type, ghost_type)); for(; p_it != p_end; ++p_it, ++linearized_el) { partition_list[linearized_el] = (*p_it)(0); } } fillPartitionInformation(mesh, partition_list); delete[] partition_list; + restoreConnectivity(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartitionMeshData::reorder() { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void MeshPartitionMeshData::setPartitionMapping(const ByElementTypeArray & mapping) { partition_mapping = &mapping; } /* -------------------------------------------------------------------------- */ void MeshPartitionMeshData::setPartitionMappingFromMeshData(const std::string & data_name) { partition_mapping = &(mesh.getData(data_name)); } __END_AKANTU__ diff --git a/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc b/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc index eb77fc516..c223a5d05 100644 --- a/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc +++ b/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc @@ -1,461 +1,466 @@ /** * @file mesh_partition_scotch.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date Tue Aug 17 16:19:43 2010 * * @brief implementation of the MeshPartitionScotch class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "mesh_partition_scotch.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MeshPartitionScotch::MeshPartitionScotch(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id) : MeshPartition(mesh, spatial_dimension, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SCOTCH_Mesh * MeshPartitionScotch::createMesh() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes = mesh.getNbNodes(); UInt total_nb_element = 0; UInt nb_edge = 0; Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator end = mesh.lastType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); total_nb_element += nb_element; nb_edge += nb_element * nb_nodes_per_element; } SCOTCH_Num vnodbas = 0; SCOTCH_Num vnodnbr = nb_nodes; SCOTCH_Num velmbas = vnodnbr; SCOTCH_Num velmnbr = total_nb_element; SCOTCH_Num * verttab = new SCOTCH_Num[vnodnbr + velmnbr + 1]; SCOTCH_Num * vendtab = verttab + 1; SCOTCH_Num * velotab = NULL; SCOTCH_Num * vnlotab = NULL; SCOTCH_Num * vlbltab = NULL; memset(verttab, 0, (vnodnbr + velmnbr + 1) * sizeof(SCOTCH_Num)); it = mesh.firstType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, _not_ghost); /// count number of occurrence of each node for (UInt el = 0; el < nb_element; ++el) { UInt * conn_val = connectivity.values + el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { verttab[*(conn_val++)]++; } } } /// convert the occurrence array in a csr one for (UInt i = 1; i < nb_nodes; ++i) verttab[i] += verttab[i-1]; for (UInt i = nb_nodes; i > 0; --i) verttab[i] = verttab[i-1]; verttab[0] = 0; /// rearrange element to get the node-element list SCOTCH_Num edgenbr = verttab[vnodnbr] + nb_edge; SCOTCH_Num * edgetab = new SCOTCH_Num[edgenbr]; UInt linearized_el = 0; it = mesh.firstType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt * conn_val = connectivity.values + el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) edgetab[verttab[*(conn_val++)]++] = linearized_el + velmbas; } } for (UInt i = nb_nodes; i > 0; --i) verttab[i] = verttab[i-1]; verttab[0] = 0; SCOTCH_Num * verttab_tmp = verttab + vnodnbr + 1; SCOTCH_Num * edgetab_tmp = edgetab + verttab[vnodnbr]; it = mesh.firstType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { *verttab_tmp = *(verttab_tmp - 1) + nb_nodes_per_element; verttab_tmp++; UInt * conn = connectivity.values + el * nb_nodes_per_element; for (UInt i = 0; i < nb_nodes_per_element; ++i) { *(edgetab_tmp++) = *(conn++) + vnodbas; } } } SCOTCH_Mesh * meshptr = new SCOTCH_Mesh; SCOTCH_meshInit(meshptr); SCOTCH_meshBuild(meshptr, velmbas, vnodbas, velmnbr, vnodnbr, verttab, vendtab, velotab, vnlotab, vlbltab, edgenbr, edgetab); /// Check the mesh AKANTU_DEBUG_ASSERT(SCOTCH_meshCheck(meshptr) == 0, "Scotch mesh is not consistent"); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { /// save initial graph FILE *fmesh = fopen("ScotchMesh.msh", "w"); SCOTCH_meshSave(meshptr, fmesh); fclose(fmesh); /// write geometry file std::ofstream fgeominit; fgeominit.open("ScotchMesh.xyz"); fgeominit << spatial_dimension << std::endl << nb_nodes << std::endl; const Array & nodes = mesh.getNodes(); Real * nodes_val = nodes.values; for (UInt i = 0; i < nb_nodes; ++i) { fgeominit << i << " "; for (UInt s = 0; s < spatial_dimension; ++s) fgeominit << *(nodes_val++) << " "; fgeominit << std::endl;; } fgeominit.close(); } #endif AKANTU_DEBUG_OUT(); return meshptr; } /* -------------------------------------------------------------------------- */ void MeshPartitionScotch::destroyMesh(SCOTCH_Mesh * meshptr) { AKANTU_DEBUG_IN(); SCOTCH_Num velmbas, vnodbas, vnodnbr, velmnbr, *verttab, *vendtab, *velotab, *vnlotab, *vlbltab, edgenbr, *edgetab, degrptr; SCOTCH_meshData(meshptr, &velmbas, &vnodbas, &velmnbr, &vnodnbr, &verttab, &vendtab, &velotab, &vnlotab, &vlbltab, &edgenbr, &edgetab, °rptr); delete [] verttab; delete [] edgetab; SCOTCH_meshExit(meshptr); delete meshptr; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartitionScotch::partitionate(UInt nb_part, const EdgeLoadFunctor & edge_load_func, const Array & pairs) { AKANTU_DEBUG_IN(); nb_partitions = nb_part; + tweakConnectivity(pairs); + AKANTU_DEBUG_INFO("Partitioning the mesh " << mesh.getID() << " in " << nb_part << " parts."); Array dxadj; Array dadjncy; Array edge_loads; - buildDualGraph(dxadj, dadjncy, edge_loads, edge_load_func, pairs); + buildDualGraph(dxadj, dadjncy, edge_loads, edge_load_func); /// variables that will hold our structures in scotch format SCOTCH_Graph scotch_graph; SCOTCH_Strat scotch_strat; /// description number and arrays for struct mesh for scotch SCOTCH_Num baseval = 0; //base numbering for element and //nodes (0 -> C , 1 -> fortran) SCOTCH_Num vertnbr = dxadj.getSize() - 1; //number of vertexes SCOTCH_Num *parttab; //array of partitions SCOTCH_Num edgenbr = dxadj(vertnbr); //twice the number of "edges" //(an "edge" bounds two nodes) SCOTCH_Num * verttab = dxadj.storage(); //array of start indices in edgetab SCOTCH_Num * vendtab = NULL; //array of after-last indices in edgetab SCOTCH_Num * velotab = NULL; //integer load associated with //every vertex ( optional ) SCOTCH_Num * edlotab = edge_loads.storage(); //integer load associated with //every edge ( optional ) SCOTCH_Num * edgetab = dadjncy.storage(); // adjacency array of every vertex SCOTCH_Num * vlbltab = NULL; // vertex label array (optional) /// Allocate space for Scotch arrays parttab = new SCOTCH_Num[vertnbr]; /// Initialize the strategy structure SCOTCH_stratInit (&scotch_strat); /// Initialize the graph structure SCOTCH_graphInit(&scotch_graph); /// Build the graph from the adjacency arrays SCOTCH_graphBuild(&scotch_graph, baseval, vertnbr, verttab, vendtab, velotab, vlbltab, edgenbr, edgetab, edlotab); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { /// save initial graph FILE *fgraphinit = fopen("GraphIniFile.grf", "w"); SCOTCH_graphSave(&scotch_graph, fgraphinit); fclose(fgraphinit); /// write geometry file std::ofstream fgeominit; fgeominit.open("GeomIniFile.xyz"); fgeominit << spatial_dimension << std::endl << vertnbr << std::endl; const Array & nodes = mesh.getNodes(); Mesh::type_iterator f_it = mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined); Mesh::type_iterator f_end = mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined); UInt out_linerized_el = 0; for(; f_it != f_end; ++f_it) { ElementType type = *f_it; UInt nb_element = mesh.getNbElement(*f_it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type); Real mid[spatial_dimension] ; for (UInt el = 0; el < nb_element; ++el) { memset(mid, 0, spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity.values[nb_nodes_per_element * el + n]; for (UInt s = 0; s < spatial_dimension; ++s) mid[s] += ((Real) nodes.values[node * spatial_dimension + s]) / ((Real) nb_nodes_per_element); } fgeominit << out_linerized_el++ << " "; for (UInt s = 0; s < spatial_dimension; ++s) fgeominit << mid[s] << " "; fgeominit << std::endl;; } } fgeominit.close(); } #endif /// Check the graph AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0, "Graph to partition is not consistent"); /// Partition the mesh SCOTCH_graphPart(&scotch_graph, nb_part, &scotch_strat, parttab); /// Check the graph AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0, "Partitioned graph is not consistent"); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { /// save the partitioned graph FILE *fgraph=fopen("GraphFile.grf", "w"); SCOTCH_graphSave(&scotch_graph, fgraph); fclose(fgraph); /// save the partition map std::ofstream fmap; fmap.open("MapFile.map"); fmap << vertnbr << std::endl; for (Int i = 0; i < vertnbr; i++) fmap << i << " " << parttab[i] << std::endl; fmap.close(); } #endif /// free the scotch data structures SCOTCH_stratExit(&scotch_strat); SCOTCH_graphFree(&scotch_graph); SCOTCH_graphExit(&scotch_graph); fillPartitionInformation(mesh, parttab); delete [] parttab; + + restoreConnectivity(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartitionScotch::reorder() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Reordering the mesh " << mesh.getID()); SCOTCH_Mesh * scotch_mesh = createMesh(); UInt nb_nodes = mesh.getNbNodes(); SCOTCH_Strat scotch_strat; // SCOTCH_Ordering scotch_order; SCOTCH_Num * permtab = new SCOTCH_Num[nb_nodes]; SCOTCH_Num * peritab = NULL; SCOTCH_Num cblknbr = 0; SCOTCH_Num * rangtab = NULL; SCOTCH_Num * treetab = NULL; /// Initialize the strategy structure SCOTCH_stratInit (&scotch_strat); SCOTCH_Graph scotch_graph; SCOTCH_graphInit(&scotch_graph); SCOTCH_meshGraph(scotch_mesh, &scotch_graph); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { FILE *fgraphinit = fopen("ScotchMesh.grf", "w"); SCOTCH_graphSave(&scotch_graph,fgraphinit); fclose(fgraphinit); } #endif /// Check the graph // AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0, // "Mesh to Graph is not consistent"); SCOTCH_graphOrder(&scotch_graph, &scotch_strat, permtab, peritab, &cblknbr, rangtab, treetab); SCOTCH_graphExit(&scotch_graph); SCOTCH_stratExit(&scotch_strat); destroyMesh(scotch_mesh); /// Renumbering UInt spatial_dimension = mesh.getSpatialDimension(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(_all_dimensions, gt); Mesh::type_iterator end = mesh.lastType(_all_dimensions, gt); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, gt); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, gt); UInt * conn = connectivity.values; for (UInt el = 0; el < nb_element * nb_nodes_per_element; ++el, ++conn) { *conn = permtab[*conn]; } } } /// \todo think of a in-place way to do it Real * new_coordinates = new Real[spatial_dimension * nb_nodes]; Real * old_coordinates = mesh.getNodes().values; for (UInt i = 0; i < nb_nodes; ++i) { memcpy(new_coordinates + permtab[i] * spatial_dimension, old_coordinates + i * spatial_dimension, spatial_dimension * sizeof(Real)); } memcpy(old_coordinates, new_coordinates, nb_nodes * spatial_dimension * sizeof(Real)); delete [] new_coordinates; delete [] permtab; AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index 2cdc7a117..8be944131 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,1583 +1,1584 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Leonardo Snozzi * @author Nicolas Richart * @author Marco Vocialta * * @date Fri Aug 20 12:19:44 2010 * * @brief All mesh utils necessary for various tasks * * @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 "mesh_utils.hh" #include "aka_safe_enum.hh" #include "fem.hh" __BEGIN_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(); AKANTU_DEBUG_ASSERT(mesh.firstType(spatial_dimension) != mesh.lastType(spatial_dimension), "Some elements must be found in right dimension to compute facets!"); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt); for (; first != last; ++first) { ElementType type = *first; UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(type, *gt).begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) for (UInt n = 0; n < conn_it->size(); ++n) ++node_to_elem.rowOffset((*conn_it)(n)); } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list Element e; node_to_elem.beginInsertions(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt); e.ghost_type = *gt; for (; first != last; ++first) { ElementType type = *first; e.type = type; UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(type, *gt).begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; for (UInt n = 0; n < conn_it->size(); ++n) node_to_elem.insertInRow((*conn_it)(n), e); } } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * This function should disappear in the future */ 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(); UInt nb_nodes = mesh.getNbNodes(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).values; nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); nb_good_types++; } AKANTU_DEBUG_ASSERT(nb_good_types != 0, "Some elements must be found in right dimension to compute facets!"); /// 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 t = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) { ++node_to_elem.rowOffset(conn_val[t][el_offset + n]); } } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// rearrange element to get the node-element list for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) node_to_elem.insertInRow(conn_val[t][el_offset + n], linearized_el); } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsByElementType(const Mesh & mesh, - ElementType type, - CSR & node_to_elem) { + CSR & node_to_elem, + const ElementType & type, + const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - UInt nb_elements = mesh.getConnectivity(type, _not_ghost).getSize(); + UInt nb_elements = mesh.getConnectivity(type, ghost_type).getSize(); - UInt * conn_val = mesh.getConnectivity(type, _not_ghost).values; + 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(); ByElementTypeReal barycenter; ByElementTypeUInt prank_to_element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { mesh.getConnectivity(*it, *gt).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension, barycenter, prank_to_element); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(Mesh & mesh, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); ByElementTypeUInt prank_to_element; buildAllFacetsParallel(mesh, mesh_facets, prank_to_element); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacetsParallel(Mesh & mesh, Mesh & mesh_facets, ByElementTypeUInt & prank_to_element) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); ByElementTypeReal barycenter; /// generate facets buildFacetsDimension(mesh, mesh_facets, false, spatial_dimension, barycenter, prank_to_element); /// compute their barycenters mesh_facets.initByElementTypeArray(barycenter, spatial_dimension, spatial_dimension - 1); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, *gt); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, *gt); for(; it != end; ++it) { UInt nb_element = mesh_facets.getNbElement(*it, *gt); barycenter(*it, *gt).resize(nb_element); Array::iterator< Vector > bary = barycenter(*it, *gt).begin(spatial_dimension); Array::iterator< Vector > bary_end = barycenter(*it, *gt).end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { mesh_facets.getBarycenter(el, *it, bary->storage(), *gt); } } } /// copy nodes type pointer mesh_facets.nodes_type = mesh.nodes_type; /// sort facets and generate subfacets for (UInt i = spatial_dimension - 1; i > 0; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i, barycenter, prank_to_element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension, ByElementTypeReal & barycenter, ByElementTypeUInt & prank_to_element){ AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const Array::const_iterator< Vector > mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator first = mesh.firstType(dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(dimension, ghost_type); for(; first != last; ++first) { ElementType type = *first; if(mesh.getSpatialDimension(type) != dimension) continue; ElementType facet_type = mesh.getFacetType(type); UInt nb_element = mesh.getNbElement(type, ghost_type); // getting connectivity of boundary facets Array * connectivity_facets = mesh_facets.getConnectivityPointer(facet_type, ghost_type); connectivity_facets->resize(0); Array< std::vector > * element_to_subelement = mesh_facets.getElementToSubelementPointer(facet_type, ghost_type); element_to_subelement->resize(0); Array * subelement_to_element = mesh_facets.getSubelementToElementPointer(type, ghost_type); subelement_to_element->resize(nb_element); } } CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; const Array * nodes_type = NULL; nodes_type = &(mesh.getNodesType()); Element current_element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; Mesh::type_iterator first = mesh.firstType(dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(dimension, ghost_type); for(; first != last; ++first) { ElementType type = *first; ElementType facet_type = mesh.getFacetType(type); current_element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array< std::vector > * element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, ghost_type); Array * connectivity_facets = &mesh_facets.getConnectivity(facet_type, ghost_type); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; Matrix facets = mesh.getFacetConnectivity(el, type, ghost_type); UInt nb_nodes_per_facet = facets.cols(); for (UInt f = 0; f < facets.rows(); ++f) { Vector facet(nb_nodes_per_facet); for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = facets(f, n); UInt first_node_nb_elements = node_to_elem.getNbCols(facets(f, 0)); counter.resize(first_node_nb_elements); counter.clear(); //loop over the other nodes to search intersecting elements, //which are the elements that share another node with the //starting element after first_node CSR::iterator first_node_elements = node_to_elem.begin(facet(0)); CSR::iterator first_node_elements_end = node_to_elem.end(facet(0)); UInt local_el = 0; for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { CSR::iterator node_elements_begin = node_to_elem.begin(facet(n)); CSR::iterator node_elements_end = node_to_elem.end (facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; Array connected_elements; for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (counter(el_f) == nb_nodes_per_facet - 1) { ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } } if (minimum_el == current_element) { bool full_ghost_facet = false; if (nodes_type != NULL) { UInt n = 0; while (n < nb_nodes_per_facet && (*nodes_type)(facet(n)) == -3) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; } if (!full_ghost_facet) { if (!boundary_only || (boundary_only && nb_element_connected_to_facet == 1)) { std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); /// boundary facet if (nb_element_connected_to_facet == 1) elements.push_back(ElementNull); /// internal facet else if (nb_element_connected_to_facet == 2) { elements.push_back(connected_elements(1)); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2]; for (UInt el = 0; el < connected_elements.getSize(); ++el) gt[el] = connected_elements(el).ghost_type; if (gt[0] + gt[1] == 1) { try { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { UInt current_el = connected_elements(el).element; ElementType current_type = connected_elements(el).type; GhostType current_gt = connected_elements(el).ghost_type; Array & prank_to_el = prank_to_element(current_type, current_gt); prank[el] = prank_to_el(current_el); } bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, facet_ghost_type); } catch (...) { } } } /// facet of facet else { for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { elements.push_back(connected_elements(i)); } /// check if sorting is needed: /// - in 3D to sort triangles around segments /// - in 2D to sort segments around points if (dimension == spatial_dimension - 1) MeshUtils::sortElements(elements, facet, mesh, mesh_facets, barycenter); } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->getSize() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type != _not_defined) { Array & subelement_to_element = mesh_facets.getSubelementToElement(type, loc_el.ghost_type); for (UInt f_in = 0; f_in < facets.rows(); ++f_in) { if (subelement_to_element(loc_el.element, f_in).type == _not_defined) { subelement_to_element(loc_el.element, f_in).type = facet_type; subelement_to_element(loc_el.element, f_in).element = current_facet; subelement_to_element(loc_el.element, f_in).ghost_type = facet_ghost_type; break; } } } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, facet_ghost_type); } } } } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, UInt * 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.getSize(); ++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); std::map::iterator it = renumbering_map.begin(); std::map::iterator end = renumbering_map.end(); old_nodes_numbers.resize(renumbering_map.size()); for (;it != end; ++it) { old_nodes_numbers(it->second) = it->first; } renumbering_map.clear(); /// copy the renumbered connectivity to the right place Array * local_conn = mesh.getConnectivityPointer(type); local_conn->resize(nb_local_element); memcpy(local_conn->values, local_connectivities, nb_local_element * nb_nodes_per_element * sizeof(UInt)); Array * ghost_conn = mesh.getConnectivityPointer(type, _ghost); ghost_conn->resize(nb_ghost_element); memcpy(ghost_conn->values, local_connectivities + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity(UInt * list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes; UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; std::map::iterator it = renumbering_map.find(node); if(it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity_vect = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity_vect.getSize()); UInt * connectivity = connectivity_vect.storage(); 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)); std::map::iterator it = renumbering_map.begin(); std::map::iterator end = renumbering_map.end(); for (; it != end; ++it) { new_numbering(it->first) = it->second; } for (UInt i = 0; i < new_numbering.getSize(); ++i) { if(new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::insertIntrinsicCohesiveElementsInArea(Mesh & mesh, const Array & limits) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(limits.getNbComponent() == 2, "Number of components for limits array must be 2"); UInt spatial_dimension = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(limits.getSize() == spatial_dimension, "Limits array size must be equal to spatial dimension"); Real * bary_facet = new Real[spatial_dimension]; Mesh mesh_facets(spatial_dimension, mesh.getNodes(), "mesh_facets"); buildAllFacets(mesh, mesh_facets); Array facet_insertion; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1); for(; it != end; ++it) { const ElementType type_facet = *it; UInt nb_facet = mesh_facets.getNbElement(type_facet); facet_insertion.resize(nb_facet); facet_insertion.clear(); for (UInt f = 0; f < nb_facet; ++f) { mesh_facets.getBarycenter(f, type_facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet[coord_in_limit] > limits(coord_in_limit, 0) && bary_facet[coord_in_limit] < limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit == spatial_dimension) facet_insertion(f) = true; } insertIntrinsicCohesiveElements(mesh, mesh_facets, type_facet, facet_insertion); } delete [] bary_facet; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::insertIntrinsicCohesiveElements(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, const Array & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// create dummy doubled data containers ByElementTypeUInt doubled_facets("doubled_facets", ""); mesh_facets.initByElementTypeArray(doubled_facets, 2, spatial_dimension - 1); /// setup byelementtype insertion data ByElementTypeArray facet_ins("facet_ins", ""); mesh_facets.initByElementTypeArray(facet_ins, 1, spatial_dimension - 1); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1); for(; it != end; ++it) { if (*it != type_facet) continue; Array & f_ins = facet_ins(type_facet); f_ins.copy(facet_insertion); } /// add cohesive connectivity type ElementType type_cohesive = FEM::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive); /// insert cohesive elements insertCohesiveElements(mesh, mesh_facets, facet_ins, doubled_facets, false); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::insertCohesiveElements(Mesh & mesh, Mesh & mesh_facets, const ByElementTypeArray & facet_insertion, ByElementTypeUInt & doubled_facets, const bool extrinsic) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); NewNodesEvent node_event; Array & doubled_nodes = node_event.getList(); doubled_nodes.extendComponentsInterlaced(2, 1); UInt total_nb_new_elements = 0; NewElementsEvent element_event; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & doubled_f = doubled_facets(type_facet, gt_facet); doubled_f.resize(0); const Array & f_insertion = facet_insertion(type_facet, gt_facet); Array facets_to_be_doubled; for (UInt f = 0; f < f_insertion.getSize(); ++f) { if (f_insertion(f)) facets_to_be_doubled.push_back(f); } if(facets_to_be_doubled.getSize() == 0) continue; Element facet(type_facet, 0, gt_facet); /// update mesh for (UInt f = 0; f < facets_to_be_doubled.getSize(); ++f) { facet.element = facets_to_be_doubled(f); doubleFacet(mesh, mesh_facets, facet, doubled_nodes, doubled_f); } /// double middle nodes if it's the case if (type_facet == _segment_3) doubleMiddleNode(mesh, mesh_facets, gt_facet, doubled_nodes, doubled_f); /// loop over doubled facets to insert cohesive elements ElementType type_cohesive = FEM::getCohesiveElementType(type_facet); Array & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet); const Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive.getSize(); UInt new_nb_cohesive_elements = old_nb_cohesive_elements + doubled_f.getSize(); conn_cohesive.resize(new_nb_cohesive_elements); Array * facet_to_coh_element = mesh.getSubelementToElementPointer(type_cohesive, gt_facet); facet_to_coh_element->resize(new_nb_cohesive_elements); Element first_facet(type_facet, 0, gt_facet, _ek_regular); Element second_facet(type_facet, 0, gt_facet, _ek_regular); for (UInt f = 0; f < doubled_f.getSize(); ++f) { UInt nb_cohesive_elements = old_nb_cohesive_elements + f; first_facet.element = doubled_f(f, 0); second_facet.element = doubled_f(f, 1); /// copy first facet's connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { conn_cohesive(nb_cohesive_elements, n) = conn_facet(first_facet.element, n); conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet) = conn_facet(second_facet.element, n); } /// update element_to_facet vectors Element cohesive_element(type_cohesive, nb_cohesive_elements, gt_facet, _ek_cohesive); element_to_facet(first_facet.element)[1] = cohesive_element; element_to_facet(second_facet.element)[1] = cohesive_element; /// update facet to cohesive element vector (*facet_to_coh_element)(nb_cohesive_elements, 0) = first_facet; (*facet_to_coh_element)(nb_cohesive_elements, 1) = second_facet; } total_nb_new_elements += new_nb_cohesive_elements; } } UInt nb_new_nodes = doubled_nodes.getSize(); UInt total_nb_new_nodes = nb_new_nodes; /// update node global id StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); if (psize > 1 && extrinsic) { Array & nodes_type = const_cast &>(mesh.getNodesType()); UInt nb_old_nodes = nodes_type.getSize(); nodes_type.resize(nb_old_nodes + nb_new_nodes); /// update nodes' type for (UInt n = 0; n < nb_new_nodes; ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); nodes_type(new_node) = nodes_type(old_node); } comm.allReduce(&total_nb_new_nodes, 1, _so_sum); comm.allReduce(&total_nb_new_elements, 1, _so_sum); Array & nodes_global_ids = const_cast &>(mesh.getGlobalNodesIds()); nodes_global_ids.resize(nb_old_nodes + nb_new_nodes); } if (total_nb_new_nodes > 0 || !extrinsic) { mesh.sendEvent(node_event); } if (total_nb_new_elements > 0 || !extrinsic) { mesh.updateTypesOffsets(_not_ghost); mesh.sendEvent(element_event); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleMiddleNode(Mesh & mesh, Mesh & mesh_facets, GhostType gt_facet, Array & doubled_nodes, const Array & doubled_facets) { AKANTU_DEBUG_IN(); ElementType type_facet = _segment_3; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); const Array< std::vector > & elem_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt nb_new_facets = doubled_facets.getSize(); UInt old_nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(old_nb_doubled_nodes + nb_new_facets); UInt old_nb_nodes = position.getSize(); position.resize(old_nb_nodes + nb_new_facets); for (UInt f = 0; f < nb_new_facets; ++f) { UInt facet_first = doubled_facets(f, 0); UInt facet_second = doubled_facets(f, 1); UInt new_node = old_nb_nodes + f; /// store doubled nodes UInt nb_doubled_nodes = old_nb_doubled_nodes + f; UInt old_node = conn_facet(facet_first, 2); doubled_nodes(nb_doubled_nodes, 0) = old_node; doubled_nodes(nb_doubled_nodes, 1) = new_node; /// update position for (UInt dim = 0; dim < spatial_dimension; ++dim) position(new_node, dim) = position(old_node, dim); /// update facet connectivity conn_facet(facet_second, 2) = new_node; /// update element connectivity for (UInt el = 0; el < elem_to_facet(facet_second).size(); ++el) { Element elem = elem_to_facet(facet_second)[el]; ElementType type_elem = elem.type; if (type_elem != _not_defined) { UInt elem_global = elem.element; GhostType gt_elem = elem.ghost_type; Array & conn_elem = mesh.getConnectivity(type_elem, gt_elem); for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) { if (conn_elem(elem_global, n) == old_node) { conn_elem(elem_global, n) = new_node; break; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, Element & facet, Array & doubled_nodes, Array & doubled_facets) { AKANTU_DEBUG_IN(); const UInt f_index = facet.element; const ElementType type_facet = facet.type; const GhostType gt_facet = facet.ghost_type; const ElementType type_subfacet = mesh.getFacetType(type_facet); const UInt nb_subfacet = mesh.getNbFacetsPerElement(type_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); if (element_to_facet(f_index)[1].type == _not_defined || element_to_facet(f_index)[1].kind == _ek_cohesive) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); return; } /// adding a new facet by copying original one /// create new connectivity Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); UInt nb_facet = conn_facet.getSize(); conn_facet.resize(nb_facet + 1); for (UInt n = 0; n < conn_facet.getNbComponent(); ++n) conn_facet(nb_facet, n) = conn_facet(f_index, n); /// store doubled facets UInt nb_doubled_facets = doubled_facets.getSize(); doubled_facets.resize(nb_doubled_facets + 1); doubled_facets(nb_doubled_facets, 0) = f_index; doubled_facets(nb_doubled_facets, 1) = nb_facet; /// update elements connected to facet std::vector first_facet_list = element_to_facet(f_index); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(nb_facet)[0] = element_to_facet(nb_facet)[1]; element_to_facet(f_index)[1] = ElementNull; element_to_facet(nb_facet)[1] = ElementNull; /// update facet_to_element vector ElementType type = element_to_facet(nb_facet)[0].type; UInt el = element_to_facet(nb_facet)[0].element; GhostType gt = element_to_facet(nb_facet)[0].ghost_type; Array & facet_to_element = mesh_facets.getSubelementToElement(type, gt); UInt i; for (i = 0; facet_to_element(el, i) != facet && i <= facet_to_element.getNbComponent(); ++i); facet_to_element(el, i).element = nb_facet; /// create new links to subfacets and update list of facets /// connected to subfacets if (type_facet == _point_1) { Array & position = mesh.getNodes(); UInt old_node = conn_facet(f_index); UInt new_node = position.getSize(); /// update position position.resize(new_node + 1); position(new_node) = position(old_node); conn_facet(nb_facet) = new_node; Array & conn_segment = mesh.getConnectivity(type, gt); UInt nb_nodes_per_segment = conn_facet.getNbComponent(); /// update facet connectivity UInt i; for (i = 0; conn_segment(el, i) != old_node && i <= nb_nodes_per_segment; ++i); conn_segment(el, i) = new_node; Vector d_nodes(2); d_nodes(0) = old_node; d_nodes(1) = new_node; doubled_nodes.push_back(d_nodes); } else { Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); subfacet_to_facet.resize(nb_facet + 1); for (UInt sf = 0; sf < nb_subfacet; ++sf) { subfacet_to_facet(nb_facet, sf) = subfacet_to_facet(f_index, sf); Element subfacet = subfacet_to_facet(f_index, sf); if (subfacet == ElementNull) continue; UInt sf_index = subfacet.element; GhostType sf_gt = subfacet.ghost_type; Array< std::vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, sf_gt); /// find index to start looping around facets connected to current /// subfacet UInt start = 0; UInt nb_connected_facets = facet_to_subfacet(sf_index).size(); while (facet_to_subfacet(sf_index)[start] != facet && start <= facet_to_subfacet(sf_index).size()) ++start; /// add the new facet to the list next to the original one ++nb_connected_facets; facet_to_subfacet(sf_index).resize(nb_connected_facets); for (UInt f = nb_connected_facets - 1; f > start; --f) { facet_to_subfacet(sf_index)[f] = facet_to_subfacet(sf_index)[f - 1]; } /// check if the new facet should be inserted before or after the /// original one: the second element connected to both original /// and new facet will be _not_defined, so I check if the first /// one is equal to one of the elements connected to the following /// facet in the facet_to_subfacet vector UInt f_start = facet_to_subfacet(sf_index)[start].element; Element facet_next; UInt index_next = start + 2; if (index_next == nb_connected_facets) index_next = 0; facet_next = facet_to_subfacet(sf_index)[index_next]; while (facet_next.type == _not_defined) { ++index_next; if (index_next == nb_connected_facets) index_next = 0; facet_next = facet_to_subfacet(sf_index)[index_next]; } Array< std::vector > & element_to_facet_next = mesh_facets.getElementToSubelement(facet_next.type, facet_next.ghost_type); UInt f_next = facet_next.element; if ((element_to_facet(f_start)[0] == element_to_facet_next(f_next)[0]) || ( element_to_facet(f_start)[0] == element_to_facet_next(f_next)[1])) facet_to_subfacet(sf_index)[start].element = nb_facet; else facet_to_subfacet(sf_index)[start + 1].element = nb_facet; /// loop on every facet connected to the current subfacet for (UInt f = start + 2; ; ++f) { /// reset f in order to continue looping from the beginning if (f == nb_connected_facets) f = 0; /// exit loop if it reaches the end if (f == start) break; /// if current loop facet is on the boundary, double subfacet Element f_global = facet_to_subfacet(sf_index)[f]; if (f_global.type == _not_defined) continue; Array< std::vector > & el_to_f = mesh_facets.getElementToSubelement(f_global.type, f_global.ghost_type); if (el_to_f(f_global.element)[1].type == _not_defined || el_to_f(f_global.element)[1].kind == _ek_cohesive) { doubleSubfacet(mesh, mesh_facets, subfacet, start, f, doubled_nodes); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, const Element & subfacet, UInt start, UInt end, Array & doubled_nodes) { AKANTU_DEBUG_IN(); const UInt sf_index = subfacet.element; const ElementType type_subfacet = subfacet.type; const GhostType gt_subfacet = subfacet.ghost_type; Array< std::vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); UInt nb_subfacet = facet_to_subfacet.getSize(); Array & conn_point = mesh_facets.getConnectivity(_point_1, gt_subfacet); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); /// add the new subfacet if (spatial_dimension == 3) AKANTU_DEBUG_TO_IMPLEMENT(); UInt new_node = position.getSize(); /// add new node in connectivity UInt new_subfacet = conn_point.getSize(); conn_point.resize(new_subfacet + 1); conn_point(new_subfacet) = new_node; /// store doubled nodes UInt nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(nb_doubled_nodes + 1); UInt old_node = doubled_nodes(nb_doubled_nodes, 0) = conn_point(sf_index); doubled_nodes(nb_doubled_nodes, 1) = new_node; /// update position position.resize(new_node + 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) position(new_node, dim) = position(old_node, dim); /// create a vector for the new subfacet in facet_to_subfacet facet_to_subfacet.resize(nb_subfacet + 1); UInt nb_connected_facets = facet_to_subfacet(sf_index).size(); /// loop over facets from start to end for (UInt f = start + 1; ; ++f) { /// reset f in order to continue looping from the beginning if (f == nb_connected_facets) f = 0; Element facet_el = facet_to_subfacet(sf_index)[f]; UInt f_global = facet_el.element; ElementType type_facet = facet_el.type; GhostType gt_facet = facet_el.ghost_type; if (type_facet == _not_defined) continue; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); /// update facet connectivity UInt i; for (i = 0; conn_facet(f_global, i) != old_node && i <= nb_nodes_per_facet; ++i); conn_facet(f_global, i) = new_node; /// update element connectivity const Array< std::vector > & elem_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); for (UInt el = 0; el < elem_to_facet(f_global).size(); ++el) { Element elem = elem_to_facet(f_global)[el]; ElementType type_elem = elem.type; if (type_elem != _not_defined) { UInt elem_global = elem.element; GhostType gt_elem = elem.ghost_type; Array & conn_elem = mesh.getConnectivity(type_elem, gt_elem); UInt nb_nodes_per_element = conn_elem.getNbComponent(); /// integer to do the final for loop UInt n = 0; /// cohesive elements: it's neccesary to identify the correct /// facet to be updated by looking for another node if (elem.kind == _ek_cohesive) { if (spatial_dimension == 3) AKANTU_DEBUG_TO_IMPLEMENT(); else if (spatial_dimension == 2) { Vector position_f_node(position.storage() + new_node * spatial_dimension, spatial_dimension); Vector position_coh_node(position.storage() + conn_elem(elem_global, 0) * spatial_dimension, spatial_dimension); if (position_f_node == position_coh_node) n = nb_nodes_per_facet; } } /// update connectivity for (; n < nb_nodes_per_element; ++n) { if (conn_elem(elem_global, n) == old_node) { conn_elem(elem_global, n) = new_node; break; } } } } /// update subfacet_to_facet vector Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); for (i = 0; subfacet_to_facet(f_global, i) != subfacet && i <= subfacet_to_facet.getNbComponent(); ++i); subfacet_to_facet(f_global, i).element = nb_subfacet; /// add current facet to facet_to_subfacet last position Element current_facet(type_facet, f_global, gt_facet); facet_to_subfacet(nb_subfacet).push_back(current_facet); /// exit loop if it reaches the end if (f == end) break; } /// rearrange the facets connected to the original subfacet and /// compute the new number of facets connected to it if (end < start) { for (UInt f = 0; f < start - end; ++f) facet_to_subfacet(sf_index)[f] = facet_to_subfacet(sf_index)[f + end + 1]; nb_connected_facets = start - end; } else { for (UInt f = 1; f < nb_connected_facets - end; ++f) facet_to_subfacet(sf_index)[start + f] = facet_to_subfacet(sf_index)[end + f]; nb_connected_facets -= end - start; } /// resize list of facets of the original subfacet facet_to_subfacet(sf_index).resize(nb_connected_facets); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::computeFacetNormals(const Mesh & mesh_facets, ByElementTypeReal & normals, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 3) AKANTU_DEBUG_TO_IMPLEMENT(); const Array & position = mesh_facets.getNodes(); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); for(; it != end; ++it) { ElementType type_facet = *it; Array & normal = normals(type_facet, ghost_type); const Array & conn = mesh_facets.getConnectivity(type_facet, ghost_type); UInt nb_nodes_per_facet = conn.getNbComponent(); UInt nb_element = conn.getSize(); normal.resize(nb_element); Array::iterator > normal_it = normal.begin(spatial_dimension); Array::iterator > normal_end = normal.end(spatial_dimension); Array::const_iterator > conn_it = conn.begin(nb_nodes_per_facet); Vector v(spatial_dimension); for (; normal_it != normal_end; ++normal_it, ++conn_it) { Vector first(position.storage() + (*conn_it)(0) * spatial_dimension, spatial_dimension); Vector second(position.storage() + (*conn_it)(1) * spatial_dimension, spatial_dimension); v = second; v -= first; Math::normal2(v.storage(), normal_it->storage()); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if(mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_WARNING("There are not facets, add them in the mesh file or call the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ByElementTypeReal barycenters; mesh.initByElementTypeArray(barycenters, spatial_dimension, _all_dimensions); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator it = mesh.firstType(_all_dimensions, *gt); Mesh::type_iterator end = mesh.lastType(_all_dimensions, *gt); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, *gt); Array & barycenters_arr = barycenters(*it, *gt); barycenters_arr.resize(nb_element); Array::iterator< Vector > bary = barycenters_arr.begin(spatial_dimension); Array::iterator< Vector > bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { mesh.getBarycenter(el, *it, bary->storage(), *gt); } } } for(Int sp(spatial_dimension); sp >= 1; --sp) { if(mesh.getNbElement(sp) == 0) continue; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(sp, *git); Mesh::type_iterator tend = mesh.lastType(sp, *git); for (;tit != tend; ++tit) { mesh.getSubelementToElementPointer(*tit, *git)->resize(mesh.getNbElement(*tit, *git)); mesh.getSubelementToElement(*tit, *git).clear(); } tit = mesh.firstType(sp - 1, *git); tend = mesh.lastType(sp - 1, *git); for (;tit != tend; ++tit) { mesh.getElementToSubelementPointer(*tit, *git)->resize(mesh.getNbElement(*tit, *git)); mesh.getElementToSubelement(*tit, *git).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(sp - 1, *git); Mesh::type_iterator tend = mesh.lastType(sp - 1, *git); facet_element.ghost_type = *git; for (;tit != tend; ++tit) { facet_element.type = *tit; Array< std::vector > & element_to_subelement = mesh.getElementToSubelement(*tit, *git); const Array & connectivity = mesh.getConnectivity(*tit, *git); Array::const_iterator< Vector > fit = connectivity.begin(mesh.getNbNodesPerElement(*tit)); Array::const_iterator< Vector > fend = connectivity.end(mesh.getNbNodesPerElement(*tit)); UInt fid = 0; for (;fit != fend; ++fit, ++fid) { const Vector & facet = *fit; facet_element.element = fid; std::map element_seen_counter; UInt nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(*tit)); for (UInt n(0); n < nb_nodes_per_facet; ++n) { CSR::iterator eit = nodes_to_elements.begin(facet(n)); CSR::iterator eend = nodes_to_elements.end(facet(n)); for(;eit != eend; ++eit) { Element & elem = *eit; std::map::iterator cit = element_seen_counter.find(elem); if(cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } std::vector connected_elements; std::map::iterator cit = element_seen_counter.begin(); std::map::iterator cend = element_seen_counter.end(); for(;cit != cend; ++cit) { if(cit->second == nb_nodes_per_facet) connected_elements.push_back(cit->first); } if(connected_elements.size() == 1) MeshUtils::sortElements(connected_elements, facet, mesh, mesh, barycenters); std::vector::iterator ceit = connected_elements.begin(); std::vector::iterator ceend = connected_elements.end(); for(;ceit != ceend; ++ceit) element_to_subelement(fid).push_back(*ceit); for (UInt ce = 0; ce < connected_elements.size(); ++ce) { Element & elem = connected_elements[ce]; Array & subelement_to_element = *(mesh.getSubelementToElementPointer(elem.type, elem.ghost_type)); UInt f(0); for(; f < mesh.getNbFacetsPerElement(elem.type) && subelement_to_element(elem.element, f) != ElementNull; ++f); AKANTU_DEBUG_ASSERT(f < mesh.getNbFacetsPerElement(elem.type), "The element "<< elem << " seems to have to many facets!!"); subelement_to_element(elem.element, f) = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Internal functions */ /* -------------------------------------------------------------------------- */ void MeshUtils::sortElements(std::vector & elements, const Vector facet, const Mesh & mesh, const Mesh & mesh_facets, const ByElementTypeReal & barycenters) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_element_connected_to_facet = elements.size(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const Array::const_iterator< Vector > mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); /// node around which the sorting is carried out is /// the first node of the current facet const Vector & first_node_coord = mesh_facets_nodes_it[facet(0)]; /// associate to each element a real value based on /// atan2 function (check wikipedia) std::map atan2; if (spatial_dimension == 3) { const Vector & second_node_coord = mesh_facets_nodes_it[facet(1)]; /// vector connecting facet first node to second Vector tangent(spatial_dimension); tangent = second_node_coord; tangent -= first_node_coord; tangent.normalize(); const Array::const_iterator< Vector > bar = barycenters(elements[0].type, elements[0].ghost_type).begin(spatial_dimension); /// vector connecting facet first node and /// barycenter of elements(0) Vector bary_coord(spatial_dimension); bary_coord.copy(bar[elements[0].element]); bary_coord -= first_node_coord; /// two normals to the segment facet to define the /// reference system Vector normal1(spatial_dimension); Vector normal2(spatial_dimension); /// get normal1 and normal2 normal1.crossProduct(tangent, bary_coord); normal1.normalize(); normal2.crossProduct(tangent, normal1); /// project the barycenter coordinates on the two /// normals to have them on the same plane atan2[elements[0]] = std::atan2(bary_coord.dot(normal2), bary_coord.dot(normal1)); for (UInt n = 1; n < nb_element_connected_to_facet; ++n) { const Array::const_iterator< Vector > bar_it = barycenters(elements[n].type, elements[n].ghost_type).begin(spatial_dimension); bary_coord.copy(bar_it[elements[n].element]); bary_coord -= first_node_coord; /// project the barycenter coordinates on the two /// normals to have them on the same plane atan2[elements[n]] = std::atan2(bary_coord.dot(normal2), bary_coord.dot(normal1)); } } else if (spatial_dimension == 2) { for (UInt n = 0; n < nb_element_connected_to_facet; ++n) { const Array::const_iterator< Vector > bar_it = barycenters(elements[n].type, elements[n].ghost_type).begin(spatial_dimension); Vector bary_coord(spatial_dimension); bary_coord.copy(bar_it[elements[n].element]); bary_coord -= first_node_coord; atan2[elements[n]] = std::atan2(bary_coord(1), bary_coord(0)); } } /// sort elements according to their atan2 values ElementSorter sorter(atan2); std::sort(elements.begin(), elements.end(), sorter); } __END_AKANTU__ // LocalWords: ElementType diff --git a/src/mesh_utils/mesh_utils.hh b/src/mesh_utils/mesh_utils.hh index 79742a078..11a7647aa 100644 --- a/src/mesh_utils/mesh_utils.hh +++ b/src/mesh_utils/mesh_utils.hh @@ -1,251 +1,254 @@ /** * @file mesh_utils.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Leonardo Snozzi * @author Nicolas Richart * @author Marco Vocialta * * @date Fri Aug 20 12:19:44 2010 * * @brief All mesh utils necessary for various tasks * * @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_common.hh" #include "mesh.hh" #include "aka_csr.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_UTILS_HH__ #define __AKANTU_MESH_UTILS_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class MeshUtils { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshUtils(); virtual ~MeshUtils(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build map from nodes to elements - static void buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension = _all_dimensions); - static void buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension = _all_dimensions); - - // static void buildNode2Elements(const Mesh & mesh, Array & node_offset, Array & node_to_elem, UInt spatial_dimension = _all_dimensions); + static void buildNode2Elements(const Mesh & mesh, + CSR & node_to_elem, + UInt spatial_dimension = _all_dimensions); + static void buildNode2Elements(const Mesh & mesh, + CSR & node_to_elem, + UInt spatial_dimension = _all_dimensions); /// build map from nodes to elements for a specific element type - static void buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, CSR & node_to_elem); - - // static void buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, Array & node_offset, Array & node_to_elem); + static void buildNode2ElementsByElementType(const Mesh & mesh, + CSR & node_to_elem, + const ElementType & type, + const GhostType & ghost_type = _not_ghost); /// build facets elements on boundary static void buildFacets(Mesh & mesh); /// build facets elements : boundary and internals (for serial simulations) static void buildAllFacets(Mesh & mesh, Mesh & mesh_facets); /// build facets elements : boundary and internals (for parallel simulations) static void buildAllFacetsParallel(Mesh & mesh, Mesh & mesh_facets, ByElementTypeUInt & prank_to_element); /// build facets for a given spatial dimension static void buildFacetsDimension(Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension, ByElementTypeReal & barycenter, ByElementTypeUInt & prank_to_element); /// build normal to some elements // static void buildNormals(Mesh & mesh, UInt spatial_dimension=0); /// take the local_connectivity array as the array of local and ghost /// connectivity, renumber the nodes and set the connectivity of the mesh static void renumberMeshNodes(Mesh & mesh, UInt * local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes); // static void setUIntData(Mesh & mesh, UInt * data, UInt nb_tags, const ElementType & type); /// Detect closed surfaces of the mesh and save the surface id /// of the surface elements in the array surface_id static void buildSurfaceID(Mesh & mesh); /// compute pbc pair for on given direction static void computePBCMap(const Mesh & mymesh,const UInt dir, std::map & pbc_pair); /// compute pbc pair for a surface pair static void computePBCMap(const Mesh & mymesh, const std::pair & surface_pair, const ElementType type, std::map & pbc_pair); // /// tweak mesh connectivity to activate pbc // static void tweakConnectivityForPBC(Mesh & mesh, // bool flag_x, // bool flag_y = false, // bool flag_z = false); /// create a multimap of nodes per surfaces static void buildNodesPerSurface(const Mesh & mesh, CSR & nodes_per_surface); /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; /// remove not connected nodes /!\ this functions renumbers the nodes. static void purifyMesh(Mesh & mesh); /// function to insert intrinsic cohesive elements in a zone /// delimited by provided limits static void insertIntrinsicCohesiveElementsInArea(Mesh & mesh, const Array & limits); /// function to insert intrinsic cohesive elements on the selected /// facets static void insertIntrinsicCohesiveElements(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, const Array & facet_insertion); /// function to insert cohesive elements on the selected facets static void insertCohesiveElements(Mesh & mesh, Mesh & mesh_facets, const ByElementTypeArray & facet_insertion, ByElementTypeUInt & doubled_facets, const bool extrinsic); /// fill the subelement to element and the elements to subelements data static void fillElementToSubElementsData(Mesh & mesh); /// compute normals on facets static void computeFacetNormals(const Mesh & mesh_facets, ByElementTypeReal & normals, GhostType ghost_type = _not_ghost); private: /// match pairs that are on the associated pbc's static void matchPBCPairs(const Mesh & mymesh, const UInt dir, std::vector & selected_left, std::vector & selected_right, std::map & pbc_pair); /// function used by all the renumbering functions static void renumberNodesInConnectivity(UInt * list_nodes, UInt nb_nodes, std::map & renumbering_map); /// function to double a given facet and update the list of doubled /// nodes static void doubleFacet(Mesh & mesh, Mesh & mesh_facets, Element & facet, Array & doubled_nodes, Array & doubled_facets); /// function to double a subfacet given start and end index for /// local facet_to_subfacet vector, and update the list of doubled /// nodes static void doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, const Element & subfacet, UInt start, UInt end, Array & doubled_nodes); /// double middle nodes if facets are _segment_3 static void doubleMiddleNode(Mesh & mesh, Mesh & mesh_facets, GhostType gt_facet, Array & doubled_nodes, const Array & doubled_facets); static void sortElements(std::vector & elements, const Vector facet, const Mesh & mesh, const Mesh & mesh_facets, const ByElementTypeReal & barycenters); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; class ElementSorter { public: ElementSorter(const ElementSorter & e) : atan2(e.atan2) {} ElementSorter(std::map & atan2) : atan2(atan2) {} inline bool operator()(const Element & first, const Element & second); private: std::map & atan2; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "mesh_utils_inline_impl.cc" #endif /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const MeshUtils & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_MESH_UTILS_HH__ */ diff --git a/src/model/solid_mechanics/contact.cc b/src/model/solid_mechanics/contact.cc index ff8ce6b9e..6ed62e53e 100644 --- a/src/model/solid_mechanics/contact.cc +++ b/src/model/solid_mechanics/contact.cc @@ -1,264 +1,264 @@ /** * @file contact.cc * * @author David Simon Kammer * @author Leonardo Snozzi * @author Nicolas Richart * * @date Fri Oct 08 15:20:20 2010 * * @brief Common part of contact classes * * @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 "contact.hh" #include "contact_search.hh" #include "aka_common.hh" #include "contact_3d_explicit.hh" #include "contact_search_explicit.hh" #include "contact_2d_explicit.hh" #include "contact_search_2d_explicit.hh" #include "contact_rigid.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Contact::Contact(const SolidMechanicsModel & model, const ContactType & type, const ContactID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), model(model), type(type), node_to_elements_offset("node_to_elements_offset", id, memory_id), node_to_elements("node_to_elements", id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Contact::~Contact() { AKANTU_DEBUG_IN(); delete contact_search; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::initContact(bool add_surfaces_flag) { AKANTU_DEBUG_IN(); Mesh & mesh = model.getFEM().getMesh(); /// Build surfaces if not done yet if(mesh.getNbSurfaces() == 0) { /* initialise nb_surfaces to zero in mesh_io */ MeshUtils::buildFacets(mesh); MeshUtils::buildSurfaceID(mesh); } /// Detect facet types const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_facet_types = 0; ElementType facet_type[_max_element_type]; UInt spatial_dimension = mesh.getSpatialDimension(); for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) != spatial_dimension) continue; facet_type[nb_facet_types] = mesh.getFacetType(type); nb_facet_types++; } CSR suface_nodes; MeshUtils::buildNodesPerSurface(mesh, suface_nodes); suface_nodes.copy(surface_to_nodes_offset, surface_to_nodes); for (UInt el_type = 0; el_type < nb_facet_types; ++el_type) { ElementType type = facet_type[el_type]; this->node_to_elements_offset.alloc(0, 1, type, _not_ghost); this->node_to_elements.alloc(0, 1, type, _not_ghost); CSR node_to_elem; - MeshUtils::buildNode2ElementsByElementType(mesh, type, node_to_elem); + MeshUtils::buildNode2ElementsByElementType(mesh, node_to_elem, type); node_to_elem.copy(node_to_elements_offset(type, _not_ghost), node_to_elements(type, _not_ghost)); } if(add_surfaces_flag) { for (UInt i = 0; i < mesh.getNbSurfaces(); ++i) { addMasterSurface(i); std::cout << "Master surface " << i << " has been added automatically" << std::endl; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::initSearch() { AKANTU_DEBUG_IN(); contact_search->initSearch(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::initNeighborStructure() { AKANTU_DEBUG_IN(); contact_search->initNeighborStructure(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::initNeighborStructure(__attribute__ ((unused)) const Surface & master_surface) { AKANTU_DEBUG_IN(); contact_search->initNeighborStructure(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::checkAndUpdate() { AKANTU_DEBUG_IN(); std::vector::iterator it; for (it = master_surfaces.begin(); it != master_surfaces.end(); ++it) { if(contact_search->checkIfUpdateStructureNeeded(*it)) { contact_search->updateStructure(*it); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::updateContact() { AKANTU_DEBUG_IN(); std::vector::iterator it; for (it = master_surfaces.begin(); it != master_surfaces.end(); ++it) { contact_search->updateStructure(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::addMasterSurface(const Surface & master_surface) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(std::find(master_surfaces.begin(), master_surfaces.end(), master_surface) == master_surfaces.end(), "Master surface already registered in the master surface list"); master_surfaces.push_back(master_surface); contact_search->addMasterSurface(master_surface); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Contact::removeMasterSurface(const Surface & master_surface) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(std::find(master_surfaces.begin(), master_surfaces.end(), master_surface) != master_surfaces.end(), "Master surface not registered in the master surface list"); std::remove(master_surfaces.begin(), master_surfaces.end(), master_surface); contact_search->removeMasterSurface(master_surface); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Contact * Contact::newContact(const SolidMechanicsModel & model, const ContactType & contact_type, const ContactSearchType & contact_search_type, const ContactNeighborStructureType & contact_neighbor_structure_type, const ContactID & id, const MemoryID & memory_id) { AKANTU_DEBUG_IN(); Contact * tmp_contact = NULL; ContactSearch * tmp_search __attribute__ ((unused)) = NULL; switch(contact_type) { case _ct_2d_expli: { tmp_contact = new Contact2dExplicit(model, contact_type, id, memory_id); break; } case _ct_3d_expli: { tmp_contact = new Contact3dExplicit(model, contact_type, id, memory_id); break; } case _ct_rigid: { tmp_contact = new ContactRigid(model, contact_type, id, memory_id); break; } case _ct_not_defined: { AKANTU_DEBUG_ERROR("Not a valid contact type : " << contact_type); break; } } std::stringstream sstr; sstr << id << ":contact_search"; switch(contact_search_type) { case _cst_2d_expli: { tmp_search = new ContactSearch2dExplicit(*tmp_contact, contact_neighbor_structure_type, contact_search_type, sstr.str()); break; } case _cst_expli: { tmp_search = new ContactSearchExplicit(*tmp_contact, contact_neighbor_structure_type, contact_search_type, sstr.str()); break; } case _cst_not_defined: { AKANTU_DEBUG_ERROR("Not a valid contact search type : " << contact_search_type); break; } } AKANTU_DEBUG_OUT(); return tmp_contact; } __END_AKANTU__