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 <dana.christen@gmail.com>
+ * @author David Kammer <david.kammer@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <sstream>
#include <algorithm>
#include <iterator>
#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<UInt> & 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<std::string, SubBoundary*>(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<UInt> 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<Int> 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<UInt>::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<DumperParaview>("paraview_"+sub.getName(),
sub.getID(),
true);
sub.addDumpFilteredMesh(mesh,
sub.elements,
sub.nodes,
mesh.getSpatialDimension() - 1,
_not_ghost,
_ek_regular);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<typename T>
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<T>(*type_it, dataset_name, *gt);
const Array<T> & dataset = mesh.getData<T>(*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<DumperParaview>("paraview_"+sub.getName(),
sub.getID(),
true);
sub.addDumpFilteredMesh(mesh,
sub.elements,
sub.nodes,
mesh.getSpatialDimension() - 1,
_not_ghost,
_ek_regular);
}
}
template void Boundary::createBoundariesFromMeshData<std::string>(const std::string & dataset_name);
template void Boundary::createBoundariesFromMeshData<UInt>(const std::string & dataset_name);
/* -------------------------------------------------------------------------- */
void Boundary::createSubBoundaryFromNodeGroup(const std::string & name,
const Array<UInt> & node_group) {
CSR<Element> node_to_elem;
MeshUtils::buildNode2Elements(mesh, node_to_elem, mesh.getSpatialDimension() - 1);
std::set<Element> 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<std::string, SubBoundary*>(name, sub_b));
}
SubBoundary & sub_bound = *boundaries_iter->second;
Array<UInt>::const_iterator<> itn = node_group.begin();
Array<UInt>::const_iterator<> endn = node_group.end();
for (;itn != endn; ++itn) {
CSR<Element>::iterator ite = node_to_elem.begin(*itn);
CSR<Element>::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<UInt>::const_iterator< Vector<UInt> > conn_it =
mesh.getConnectivity(elem.type, elem.ghost_type).begin(nb_nodes_per_element);
const Vector<UInt> & 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<DumperParaview>("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<std::string>(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<UInt> & 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); i<element_ids.getSize(); ++i) {
stream << element_ids(i) << " ";
}
}
stream << std::endl;
//std::cout << it->getElements(*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 <dana.christen@gmail.com>
+ * @author David Kammer <david.kammer@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_BOUNDARY_HH__
#define __AKANTU_BOUNDARY_HH__
#include <set>
#include "aka_common.hh"
#include "by_element_type.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
class Mesh;
/* -------------------------------------------------------------------------- */
class SubBoundary;
class Boundary {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
private:
public:
typedef std::map<std::string, SubBoundary *> BoundaryList;
typedef std::set<ElementType> BoundaryTypeSet;
/* ------------------------------------------------------------------------ */
/* SubBoundary iterator */
/* ------------------------------------------------------------------------ */
public:
template<typename T, typename container_iterator>
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 SubBoundary,
BoundaryList::const_iterator> const_iterator;
typedef internal_iterator<SubBoundary,
BoundaryList::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 <typename T> void createBoundariesFromMeshData(const std::string & dataset_name);
+ template <typename T>
+ 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<UInt> & 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 <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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<Int> & dxadj, Array<Int> & dadjncy,
Array<Int> & edge_loads,
- const EdgeLoadFunctor & edge_load_func,
- const Array<UInt> & 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<UInt> * conn[nb_types];
- Array<UInt> * conn_tmp[nb_types];
Array<Element> lin_to_element;
Element elem;
elem.ghost_type = _not_ghost;
const_cast<Mesh &>(mesh).updateTypesOffsets(_not_ghost);
const_cast<Mesh &>(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<Array<UInt> &>(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<UInt>(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<UInt> 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<UInt, UInt>::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<UInt>::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<UInt, UInt>::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<UInt, UInt>::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<UInt> 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<UInt> & 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<UInt> & 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<UInt> list_adj_part;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connectivity.values[el * nb_nodes_per_element + n];
CSR<UInt>::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<UInt>::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<UInt> & 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<UInt> & 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<Element> > & elem_to_subelem = mesh.getElementToSubelement(type, _not_ghost);
for(UInt i(0); i < mesh.getNbElement(type, _not_ghost); ++i) { // Facet loop
const std::vector<Element> & adjacent_elems = elem_to_subelem(i);
Element min_elem;
UInt min_part(std::numeric_limits<UInt>::max());
std::set<UInt> 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<UInt>::iterator git = ghost_partitions_csr(min_elem.type, _not_ghost).begin(min_elem.element);
CSR<UInt>::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<UInt>::const_iterator pit = adjacent_parts.begin();
std::set<UInt>::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<UInt> & 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<UInt> & conn = const_cast<Array<UInt> &>(mesh.getConnectivity(type, _not_ghost));
+ UInt nb_nodes_per_element = conn.getNbComponent();
+ UInt nb_element = conn.getSize();
+
+ Array<UInt> & 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<UInt> & conn = const_cast<Array<UInt> &>(mesh.getConnectivity(type, _not_ghost));
+ Array<UInt> & 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 <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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<UInt> & pairs = Array<UInt>()) = 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<Int> & dxadj, Array<Int> & dadjncy,
Array<Int> & edge_loads,
- const EdgeLoadFunctor & edge_load_func,
- const Array<UInt> & pairs);
+ const EdgeLoadFunctor & edge_load_func);
+
+ /// tweak the mesh to handle the PBC pairs
+ void tweakConnectivity(const Array<UInt> & 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<UInt> > &);
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<UInt> > ghost_partitions_csr;
ByElementTypeUInt ghost_partitions;
ByElementTypeUInt ghost_partitions_offset;
Array<UInt> * 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 <dana.christen@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#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<UInt> & 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<UInt> & 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<UInt> & partition_array = (*partition_mapping)(type, ghost_type);
Array<UInt>::const_iterator< Vector<UInt> > p_it = partition_array.begin(1);
Array<UInt>::const_iterator< Vector<UInt> > 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<UInt> & mapping) {
partition_mapping = &mapping;
}
/* -------------------------------------------------------------------------- */
void MeshPartitionMeshData::setPartitionMappingFromMeshData(const std::string & data_name) {
partition_mapping = &(mesh.getData<UInt>(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 <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <cstdio>
#include <fstream>
/* -------------------------------------------------------------------------- */
#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<UInt> & 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<UInt> & 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<UInt> & 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<Real> & 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<UInt> & pairs) {
AKANTU_DEBUG_IN();
nb_partitions = nb_part;
+ tweakConnectivity(pairs);
+
AKANTU_DEBUG_INFO("Partitioning the mesh " << mesh.getID()
<< " in " << nb_part << " parts.");
Array<Int> dxadj;
Array<Int> dadjncy;
Array<Int> 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<Real> & 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<UInt> & 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<UInt> & 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 <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh_utils.hh"
#include "aka_safe_enum.hh"
#include "fem.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
void MeshUtils::buildNode2Elements(const Mesh & mesh,
CSR<Element> & 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<UInt>::const_iterator< Vector<UInt> > 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<UInt>::const_iterator< Vector<UInt> > 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<UInt> & 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<UInt> & node_to_elem) {
+ CSR<UInt> & 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<Real>::iterator< Vector<Real> > bary
= barycenter(*it, *gt).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > 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<Real> & mesh_facets_nodes = mesh_facets.getNodes();
const Array<Real>::const_iterator< Vector<Real> > 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<UInt> * connectivity_facets =
mesh_facets.getConnectivityPointer(facet_type, ghost_type);
connectivity_facets->resize(0);
Array< std::vector<Element> > * element_to_subelement =
mesh_facets.getElementToSubelementPointer(facet_type, ghost_type);
element_to_subelement->resize(0);
Array<Element> * subelement_to_element =
mesh_facets.getSubelementToElementPointer(type, ghost_type);
subelement_to_element->resize(nb_element);
}
}
CSR<Element> node_to_elem;
buildNode2Elements(mesh, node_to_elem, dimension);
Array<UInt> counter;
const Array<Int> * 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> > * element_to_subelement =
&mesh_facets.getElementToSubelement(facet_type, ghost_type);
Array<UInt> * connectivity_facets =
&mesh_facets.getConnectivity(facet_type, ghost_type);
for (UInt el = 0; el < nb_element; ++el) {
current_element.element = el;
Matrix<UInt> facets = mesh.getFacetConnectivity(el, type, ghost_type);
UInt nb_nodes_per_facet = facets.cols();
for (UInt f = 0; f < facets.rows(); ++f) {
Vector<UInt> 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<Element>::iterator first_node_elements = node_to_elem.begin(facet(0));
CSR<Element>::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<Element>::iterator node_elements_begin = node_to_elem.begin(facet(n));
CSR<Element>::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<Element> 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<Element> 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<UInt> & 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<Element> & 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<UInt> & old_nodes_numbers) {
AKANTU_DEBUG_IN();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
std::map<UInt, UInt> 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<UInt, UInt>::iterator it = renumbering_map.begin();
std::map<UInt, UInt>::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<UInt> * 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<UInt> * 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<UInt, UInt> & 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<UInt, UInt>::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<UInt, UInt> renumbering_map;
RemovedNodesEvent remove_nodes(mesh);
Array<UInt> & 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<UInt> & 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<UInt> & new_numbering = remove_nodes.getNewNumbering();
std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1));
std::map<UInt, UInt>::iterator it = renumbering_map.begin();
std::map<UInt, UInt>::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<Real> & 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<bool> 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<bool> & 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<bool> 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<bool> & 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<bool> & facet_insertion,
ByElementTypeUInt & doubled_facets,
const bool extrinsic) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
NewNodesEvent node_event;
Array<UInt> & 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<UInt> & doubled_f = doubled_facets(type_facet, gt_facet);
doubled_f.resize(0);
const Array<bool> & f_insertion = facet_insertion(type_facet, gt_facet);
Array<UInt> 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<UInt> & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet);
const Array<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet,
gt_facet);
UInt nb_nodes_per_facet = conn_facet.getNbComponent();
Array< std::vector<Element> > & 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<Element> * 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<Int> & nodes_type = const_cast<Array<Int> &>(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<UInt> & nodes_global_ids =
const_cast<Array<UInt> &>(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<UInt> & doubled_nodes,
const Array<UInt> & doubled_facets) {
AKANTU_DEBUG_IN();
ElementType type_facet = _segment_3;
Array<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet);
Array<Real> & position = mesh.getNodes();
UInt spatial_dimension = mesh.getSpatialDimension();
const Array< std::vector<Element> > & 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<UInt> & 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<UInt> & doubled_nodes,
Array<UInt> & 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> > & 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<UInt> & 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<Element> 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<Element> & 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<Real> & 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<UInt> & 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<UInt> d_nodes(2);
d_nodes(0) = old_node;
d_nodes(1) = new_node;
doubled_nodes.push_back(d_nodes);
}
else {
Array<Element> & 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<Element> > & 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> > & 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<Element> > & 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<UInt> & 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<Element> > & facet_to_subfacet
= mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet);
UInt nb_subfacet = facet_to_subfacet.getSize();
Array<UInt> & conn_point = mesh_facets.getConnectivity(_point_1, gt_subfacet);
Array<Real> & 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<UInt> & 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<Element> > & 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<UInt> & 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<Real> position_f_node(position.storage() +
new_node * spatial_dimension,
spatial_dimension);
Vector<Real> 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<Element> & 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<Real> & 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<Real> & normal = normals(type_facet, ghost_type);
const Array<UInt> & 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<Real>::iterator<Vector<Real> > normal_it = normal.begin(spatial_dimension);
Array<Real>::iterator<Vector<Real> > normal_end = normal.end(spatial_dimension);
Array<UInt>::const_iterator<Vector<UInt> > conn_it = conn.begin(nb_nodes_per_facet);
Vector<Real> v(spatial_dimension);
for (; normal_it != normal_end; ++normal_it, ++conn_it) {
Vector<Real> first(position.storage() + (*conn_it)(0) * spatial_dimension,
spatial_dimension);
Vector<Real> 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<Real> & barycenters_arr = barycenters(*it, *gt);
barycenters_arr.resize(nb_element);
Array<Real>::iterator< Vector<Real> > bary = barycenters_arr.begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > 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<Element> 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> > & element_to_subelement = mesh.getElementToSubelement(*tit, *git);
const Array<UInt> & connectivity = mesh.getConnectivity(*tit, *git);
Array<UInt>::const_iterator< Vector<UInt> > fit = connectivity.begin(mesh.getNbNodesPerElement(*tit));
Array<UInt>::const_iterator< Vector<UInt> > fend = connectivity.end(mesh.getNbNodesPerElement(*tit));
UInt fid = 0;
for (;fit != fend; ++fit, ++fid) {
const Vector<UInt> & facet = *fit;
facet_element.element = fid;
std::map<Element, UInt> element_seen_counter;
UInt nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(*tit));
for (UInt n(0); n < nb_nodes_per_facet; ++n) {
CSR<Element>::iterator eit = nodes_to_elements.begin(facet(n));
CSR<Element>::iterator eend = nodes_to_elements.end(facet(n));
for(;eit != eend; ++eit) {
Element & elem = *eit;
std::map<Element, UInt>::iterator cit = element_seen_counter.find(elem);
if(cit != element_seen_counter.end()) {
cit->second++;
} else {
element_seen_counter[elem] = 1;
}
}
}
std::vector<Element> connected_elements;
std::map<Element, UInt>::iterator cit = element_seen_counter.begin();
std::map<Element, UInt>::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<Element>::iterator ceit = connected_elements.begin();
std::vector<Element>::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<Element> & 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<Element> & elements, const Vector<UInt> 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<Real> & mesh_facets_nodes = mesh_facets.getNodes();
const Array<Real>::const_iterator< Vector<Real> > 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<Real> & first_node_coord = mesh_facets_nodes_it[facet(0)];
/// associate to each element a real value based on
/// atan2 function (check wikipedia)
std::map<Element, Real, CompElementLess> atan2;
if (spatial_dimension == 3) {
const Vector<Real> & second_node_coord = mesh_facets_nodes_it[facet(1)];
/// vector connecting facet first node to second
Vector<Real> tangent(spatial_dimension);
tangent = second_node_coord;
tangent -= first_node_coord;
tangent.normalize();
const Array<Real>::const_iterator< Vector<Real> > bar =
barycenters(elements[0].type, elements[0].ghost_type).begin(spatial_dimension);
/// vector connecting facet first node and
/// barycenter of elements(0)
Vector<Real> 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<Real> normal1(spatial_dimension);
Vector<Real> 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<Real>::const_iterator< Vector<Real> > 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<Real>::const_iterator< Vector<Real> > bar_it =
barycenters(elements[n].type, elements[n].ghost_type).begin(spatial_dimension);
Vector<Real> 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 <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "aka_csr.hh"
/* -------------------------------------------------------------------------- */
#include <vector>
/* -------------------------------------------------------------------------- */
#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<UInt> & node_to_elem, UInt spatial_dimension = _all_dimensions);
- static void buildNode2Elements(const Mesh & mesh, CSR<Element> & node_to_elem, UInt spatial_dimension = _all_dimensions);
-
- // static void buildNode2Elements(const Mesh & mesh, Array<UInt> & node_offset, Array<UInt> & node_to_elem, UInt spatial_dimension = _all_dimensions);
+ static void buildNode2Elements(const Mesh & mesh,
+ CSR<UInt> & node_to_elem,
+ UInt spatial_dimension = _all_dimensions);
+ static void buildNode2Elements(const Mesh & mesh,
+ CSR<Element> & 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<UInt> & node_to_elem);
-
- // static void buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, Array<UInt> & node_offset, Array<UInt> & node_to_elem);
+ static void buildNode2ElementsByElementType(const Mesh & mesh,
+ CSR<UInt> & 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<UInt> & 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<UInt,UInt> & pbc_pair);
/// compute pbc pair for a surface pair
static void computePBCMap(const Mesh & mymesh,
const std::pair<Surface, Surface> & surface_pair,
const ElementType type,
std::map<UInt,UInt> & 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<UInt> & 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<Real> & limits);
/// function to insert intrinsic cohesive elements on the selected
/// facets
static void insertIntrinsicCohesiveElements(Mesh & mesh,
Mesh & mesh_facets,
ElementType type_facet,
const Array<bool> & facet_insertion);
/// function to insert cohesive elements on the selected facets
static void insertCohesiveElements(Mesh & mesh,
Mesh & mesh_facets,
const ByElementTypeArray<bool> & 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<UInt> & selected_left,
std::vector<UInt> & selected_right,
std::map<UInt,UInt> & pbc_pair);
/// function used by all the renumbering functions
static void renumberNodesInConnectivity(UInt * list_nodes,
UInt nb_nodes,
std::map<UInt, UInt> & 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<UInt> & doubled_nodes,
Array<UInt> & 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<UInt> & doubled_nodes);
/// double middle nodes if facets are _segment_3
static void doubleMiddleNode(Mesh & mesh,
Mesh & mesh_facets,
GhostType gt_facet,
Array<UInt> & doubled_nodes,
const Array<UInt> & doubled_facets);
static void sortElements(std::vector<Element> & elements, const Vector<UInt> 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<Element, Real, CompElementLess> & atan2) : atan2(atan2) {}
inline bool operator()(const Element & first, const Element & second);
private:
std::map<Element, Real, CompElementLess> & 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 <david.kammer@epfl.ch>
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <iostream>
/* -------------------------------------------------------------------------- */
__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<UInt> 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<UInt> 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<Surface>::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<Surface>::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__