diff --git a/src/common/aka_static_memory.cc b/src/common/aka_static_memory.cc
index b1e9d8efa..f592ba56a 100644
--- a/src/common/aka_static_memory.cc
+++ b/src/common/aka_static_memory.cc
@@ -1,160 +1,160 @@
/**
* @file aka_static_memory.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Nov 08 2017
*
* @brief Memory management
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <sstream>
#include <stdexcept>
/* -------------------------------------------------------------------------- */
#include "aka_static_memory.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
bool StaticMemory::is_instantiated = false;
StaticMemory * StaticMemory::single_static_memory = nullptr;
UInt StaticMemory::nb_reference = 0;
/* -------------------------------------------------------------------------- */
StaticMemory & StaticMemory::getStaticMemory() {
if (!single_static_memory) {
single_static_memory = new StaticMemory();
is_instantiated = true;
}
nb_reference++;
return *single_static_memory;
}
/* -------------------------------------------------------------------------- */
void StaticMemory::destroy() {
nb_reference--;
if (nb_reference == 0) {
delete single_static_memory;
}
}
/* -------------------------------------------------------------------------- */
StaticMemory::~StaticMemory() {
AKANTU_DEBUG_IN();
MemoryMap::iterator memory_it;
for (memory_it = memories.begin(); memory_it != memories.end(); ++memory_it) {
ArrayMap::iterator vector_it;
for (vector_it = (memory_it->second).begin();
vector_it != (memory_it->second).end(); ++vector_it) {
delete vector_it->second;
}
(memory_it->second).clear();
}
memories.clear();
is_instantiated = false;
StaticMemory::single_static_memory = nullptr;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StaticMemory::sfree(const MemoryID & memory_id, const ID & name) {
AKANTU_DEBUG_IN();
try {
auto & vectors = const_cast<ArrayMap &>(getMemory(memory_id));
ArrayMap::iterator vector_it;
vector_it = vectors.find(name);
if (vector_it != vectors.end()) {
AKANTU_DEBUG_INFO("Array " << name
<< " removed from the static memory number "
<< memory_id);
delete vector_it->second;
vectors.erase(vector_it);
AKANTU_DEBUG_OUT();
return;
}
- } catch (debug::Exception e) {
+ } catch (debug::Exception & e) {
AKANTU_EXCEPTION("The memory "
<< memory_id << " does not exist (perhaps already freed) ["
<< e.what() << "]");
AKANTU_DEBUG_OUT();
return;
}
AKANTU_DEBUG_INFO("The vector " << name
<< " does not exist (perhaps already freed)");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StaticMemory::printself(std::ostream & stream, int indent) const {
std::string space = "";
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
std::streamsize prec = stream.precision();
stream.precision(2);
stream << space << "StaticMemory [" << std::endl;
UInt nb_memories = memories.size();
stream << space << " + nb memories : " << nb_memories << std::endl;
Real tot_size = 0;
MemoryMap::const_iterator memory_it;
for (memory_it = memories.begin(); memory_it != memories.end(); ++memory_it) {
Real mem_size = 0;
stream << space << AKANTU_INDENT << "Memory [" << std::endl;
UInt mem_id = memory_it->first;
stream << space << AKANTU_INDENT << " + memory id : " << mem_id
<< std::endl;
UInt nb_vectors = (memory_it->second).size();
stream << space << AKANTU_INDENT << " + nb vectors : " << nb_vectors
<< std::endl;
stream.precision(prec);
ArrayMap::const_iterator vector_it;
for (vector_it = (memory_it->second).begin();
vector_it != (memory_it->second).end(); ++vector_it) {
(vector_it->second)->printself(stream, indent + 2);
mem_size += (vector_it->second)->getMemorySize();
}
stream << space << AKANTU_INDENT
<< " + total size : " << printMemorySize<char>(mem_size)
<< std::endl;
stream << space << AKANTU_INDENT << "]" << std::endl;
tot_size += mem_size;
}
stream << space << " + total size : " << printMemorySize<char>(tot_size)
<< std::endl;
stream << space << "]" << std::endl;
stream.precision(prec);
}
/* -------------------------------------------------------------------------- */
-} // akantu
+} // namespace akantu
diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc
index cbb7b64e0..f756132cd 100644
--- a/src/mesh/mesh.cc
+++ b/src/mesh/mesh.cc
@@ -1,585 +1,595 @@
/**
* @file mesh.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Tue Feb 20 2018
*
* @brief class handling meshes
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
/* -------------------------------------------------------------------------- */
#include "element_class.hh"
#include "group_manager_inline_impl.cc"
#include "mesh.hh"
#include "mesh_global_data_updater.hh"
#include "mesh_io.hh"
#include "mesh_iterators.hh"
#include "mesh_utils.hh"
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
#include "element_synchronizer.hh"
#include "facet_synchronizer.hh"
#include "mesh_utils_distribution.hh"
#include "node_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
#include "dumper_field.hh"
#include "dumper_internal_material_field.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <limits>
#include <sstream>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id,
Communicator & communicator)
: Memory(id, memory_id),
GroupManager(*this, id + ":group_manager", memory_id),
connectivities("connectivities", id, memory_id),
ghosts_counters("ghosts_counters", id, memory_id),
normals("normals", id, memory_id), spatial_dimension(spatial_dimension),
lower_bounds(spatial_dimension, 0.), upper_bounds(spatial_dimension, 0.),
size(spatial_dimension, 0.), local_lower_bounds(spatial_dimension, 0.),
local_upper_bounds(spatial_dimension, 0.),
mesh_data("mesh_data", id, memory_id), communicator(&communicator) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id,
const MemoryID & memory_id)
: Mesh(spatial_dimension, id, memory_id, communicator) {
AKANTU_DEBUG_IN();
this->nodes =
std::make_shared<Array<Real>>(0, spatial_dimension, id + ":coordinates");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id)
: Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id,
memory_id) {}
/* -------------------------------------------------------------------------- */
Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr<Array<Real>> & nodes,
const ID & id, const MemoryID & memory_id)
: Mesh(spatial_dimension, id, memory_id,
Communicator::getStaticCommunicator()) {
this->nodes = nodes;
this->nb_global_nodes = this->nodes->size();
this->nodes_to_elements.resize(nodes->size());
for (auto & node_set : nodes_to_elements) {
node_set = std::make_unique<std::set<Element>>();
}
this->computeBoundingBox();
}
/* -------------------------------------------------------------------------- */
void Mesh::getBarycenters(Array<Real> & barycenter, const ElementType & type,
const GhostType & ghost_type) const {
barycenter.resize(getNbElement(type, ghost_type));
for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) {
getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type},
std::get<1>(data));
}
}
/* -------------------------------------------------------------------------- */
Mesh & Mesh::initMeshFacets(const ID & id) {
AKANTU_DEBUG_IN();
if (!mesh_facets) {
mesh_facets = std::make_unique<Mesh>(spatial_dimension, this->nodes,
getID() + ":" + id, getMemoryID());
mesh_facets->mesh_parent = this;
mesh_facets->is_mesh_facets = true;
mesh_facets->nodes_type = this->nodes_type;
mesh_facets->nodes_global_ids = this->nodes_global_ids;
MeshUtils::buildAllFacets(*this, *mesh_facets, 0);
if (mesh.isDistributed()) {
mesh_facets->is_distributed = true;
mesh_facets->element_synchronizer = std::make_unique<FacetSynchronizer>(
*mesh_facets, mesh.getElementSynchronizer());
}
/// transfers the the mesh physical names to the mesh facets
if (not this->hasData("physical_names")) {
AKANTU_DEBUG_OUT();
return *mesh_facets;
}
if (not mesh_facets->hasData("physical_names")) {
mesh_facets->registerData<std::string>("physical_names");
}
auto & mesh_phys_data = this->getData<std::string>("physical_names");
auto & phys_data = mesh_facets->getData<std::string>("physical_names");
phys_data.initialize(*mesh_facets,
_spatial_dimension = spatial_dimension - 1,
_with_nb_element = true);
ElementTypeMapArray<Real> barycenters(getID(), "temporary_barycenters");
barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension,
_spatial_dimension = spatial_dimension - 1,
_with_nb_element = true);
for (auto && ghost_type : ghost_types) {
for (auto && type :
barycenters.elementTypes(spatial_dimension - 1, ghost_type)) {
mesh_facets->getBarycenters(barycenters(type, ghost_type), type,
ghost_type);
}
}
for_each_element(
mesh,
[&](auto && element) {
Vector<Real> barycenter(spatial_dimension);
mesh.getBarycenter(element, barycenter);
auto norm_barycenter = barycenter.norm();
auto tolerance = Math::getTolerance();
if (norm_barycenter > tolerance)
tolerance *= norm_barycenter;
const auto & element_to_facet = mesh_facets->getElementToSubelement(
element.type, element.ghost_type);
Vector<Real> barycenter_facet(spatial_dimension);
auto range =
enumerate(make_view(barycenters(element.type, element.ghost_type),
spatial_dimension));
#ifndef AKANTU_NDEBUG
auto min_dist = std::numeric_limits<Real>::max();
#endif
// this is a spacial search coded the most inefficient way.
auto facet =
std::find_if(range.begin(), range.end(), [&](auto && data) {
auto facet = std::get<0>(data);
if (element_to_facet(facet)[1] == ElementNull)
return false;
auto norm_distance = barycenter.distance(std::get<1>(data));
#ifndef AKANTU_NDEBUG
min_dist = std::min(min_dist, norm_distance);
#endif
return (norm_distance < tolerance);
});
if (facet == range.end()) {
AKANTU_DEBUG_INFO("The element "
<< element
<< " did not find its associated facet in the "
"mesh_facets! Try to decrease math tolerance. "
"The closest element was at a distance of "
<< min_dist);
return;
}
// set physical name
phys_data(Element{element.type, UInt(std::get<0>(*facet)),
element.ghost_type}) = mesh_phys_data(element);
},
_spatial_dimension = spatial_dimension - 1);
mesh_facets->createGroupsFromMeshData<std::string>("physical_names");
}
AKANTU_DEBUG_OUT();
return *mesh_facets;
}
/* -------------------------------------------------------------------------- */
void Mesh::defineMeshParent(const Mesh & mesh) {
AKANTU_DEBUG_IN();
this->mesh_parent = &mesh;
this->is_mesh_facets = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Mesh::~Mesh() = default;
/* -------------------------------------------------------------------------- */
void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) {
MeshIO mesh_io;
mesh_io.read(filename, *this, mesh_io_type);
type_iterator it =
this->firstType(spatial_dimension, _not_ghost, _ek_not_defined);
type_iterator last =
this->lastType(spatial_dimension, _not_ghost, _ek_not_defined);
if (it == last)
AKANTU_EXCEPTION(
"The mesh contained in the file "
<< filename << " does not seem to be of the good dimension."
<< " No element of dimension " << spatial_dimension << " where read.");
this->nb_global_nodes = this->nodes->size();
this->computeBoundingBox();
this->nodes_to_elements.resize(nodes->size());
for (auto & node_set : nodes_to_elements) {
node_set = std::make_unique<std::set<Element>>();
}
}
/* -------------------------------------------------------------------------- */
void Mesh::write(const std::string & filename,
const MeshIOType & mesh_io_type) {
MeshIO mesh_io;
mesh_io.write(filename, *this, mesh_io_type);
}
+/* -------------------------------------------------------------------------- */
+
+void Mesh::makeReady() {
+ this->nb_global_nodes = this->nodes->size();
+ this->computeBoundingBox();
+ // this->nodes_flags->resize(nodes->size(), NodeFlag::_normal);
+ this->nodes_to_elements.resize(nodes->size());
+ for (auto & node_set : nodes_to_elements) {
+ node_set = std::make_unique<std::set<Element>>();
+ }
+}
+
/* -------------------------------------------------------------------------- */
void Mesh::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "Mesh [" << std::endl;
stream << space << " + id : " << getID() << std::endl;
stream << space << " + spatial dimension : " << this->spatial_dimension
<< std::endl;
stream << space << " + nodes [" << std::endl;
nodes->printself(stream, indent + 2);
stream << space << " + connectivities [" << std::endl;
connectivities.printself(stream, indent + 2);
stream << space << " ]" << std::endl;
GroupManager::printself(stream, indent + 1);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void Mesh::computeBoundingBox() {
AKANTU_DEBUG_IN();
for (UInt k = 0; k < spatial_dimension; ++k) {
local_lower_bounds(k) = std::numeric_limits<double>::max();
local_upper_bounds(k) = -std::numeric_limits<double>::max();
}
for (UInt i = 0; i < nodes->size(); ++i) {
// if(!isPureGhostNode(i))
for (UInt k = 0; k < spatial_dimension; ++k) {
local_lower_bounds(k) = std::min(local_lower_bounds[k], (*nodes)(i, k));
local_upper_bounds(k) = std::max(local_upper_bounds[k], (*nodes)(i, k));
}
}
if (this->is_distributed) {
Matrix<Real> reduce_bounds(spatial_dimension, 2);
for (UInt k = 0; k < spatial_dimension; ++k) {
reduce_bounds(k, 0) = local_lower_bounds(k);
reduce_bounds(k, 1) = -local_upper_bounds(k);
}
communicator->allReduce(reduce_bounds, SynchronizerOperation::_min);
for (UInt k = 0; k < spatial_dimension; ++k) {
lower_bounds(k) = reduce_bounds(k, 0);
upper_bounds(k) = -reduce_bounds(k, 1);
}
} else {
this->lower_bounds = this->local_lower_bounds;
this->upper_bounds = this->local_upper_bounds;
}
size = upper_bounds - lower_bounds;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Mesh::initNormals() {
normals.initialize(*this, _nb_component = spatial_dimension,
_spatial_dimension = spatial_dimension,
_element_kind = _ek_not_defined);
}
/* -------------------------------------------------------------------------- */
void Mesh::getGlobalConnectivity(
ElementTypeMapArray<UInt> & global_connectivity) {
AKANTU_DEBUG_IN();
for (auto && ghost_type : ghost_types) {
for (auto type :
global_connectivity.elementTypes(_spatial_dimension = _all_dimensions,
_element_kind = _ek_not_defined, _ghost_type = ghost_type)) {
if (not connectivities.exists(type, ghost_type))
continue;
auto & local_conn = connectivities(type, ghost_type);
auto & g_connectivity = global_connectivity(type, ghost_type);
UInt nb_nodes = local_conn.size() * local_conn.getNbComponent();
std::transform(local_conn.begin_reinterpret(nb_nodes),
local_conn.end_reinterpret(nb_nodes),
g_connectivity.begin_reinterpret(nb_nodes),
- [&](UInt l) -> UInt {
- return this->getNodeGlobalId(l);
- });
+ [&](UInt l) -> UInt { return this->getNodeGlobalId(l); });
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name,
const std::string & group_name) {
if (group_name == "all")
return this->getDumper(dumper_name);
else
return element_groups[group_name]->getDumper(dumper_name);
}
/* -------------------------------------------------------------------------- */
template <typename T>
ElementTypeMap<UInt> Mesh::getNbDataPerElem(ElementTypeMapArray<T> & arrays,
const ElementKind & element_kind) {
ElementTypeMap<UInt> nb_data_per_elem;
for (auto type : elementTypes(spatial_dimension, _not_ghost, element_kind)) {
UInt nb_elements = this->getNbElement(type);
auto & array = arrays(type);
nb_data_per_elem(type) = array.getNbComponent() * array.size();
nb_data_per_elem(type) /= nb_elements;
}
return nb_data_per_elem;
}
/* -------------------------------------------------------------------------- */
template ElementTypeMap<UInt>
Mesh::getNbDataPerElem(ElementTypeMapArray<Real> & array,
const ElementKind & element_kind);
template ElementTypeMap<UInt>
Mesh::getNbDataPerElem(ElementTypeMapArray<UInt> & array,
const ElementKind & element_kind);
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
template <typename T>
dumper::Field *
Mesh::createFieldFromAttachedData(const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind) {
dumper::Field * field = nullptr;
ElementTypeMapArray<T> * internal = nullptr;
try {
internal = &(this->getData<T>(field_id));
} catch (...) {
return nullptr;
}
ElementTypeMap<UInt> nb_data_per_elem =
this->getNbDataPerElem(*internal, element_kind);
field = this->createElementalField<T, dumper::InternalMaterialField>(
*internal, group_name, this->spatial_dimension, element_kind,
nb_data_per_elem);
return field;
}
template dumper::Field *
Mesh::createFieldFromAttachedData<Real>(const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind);
template dumper::Field *
Mesh::createFieldFromAttachedData<UInt>(const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind);
#endif
/* -------------------------------------------------------------------------- */
void Mesh::distribute() {
this->distribute(Communicator::getStaticCommunicator());
}
/* -------------------------------------------------------------------------- */
void Mesh::distribute(Communicator & communicator) {
AKANTU_DEBUG_ASSERT(is_distributed == false,
"This mesh is already distribute");
this->communicator = &communicator;
this->element_synchronizer = std::make_unique<ElementSynchronizer>(
*this, this->getID() + ":element_synchronizer", this->getMemoryID(),
true);
this->node_synchronizer = std::make_unique<NodeSynchronizer>(
*this, this->getID() + ":node_synchronizer", this->getMemoryID(), true);
Int psize = this->communicator->getNbProc();
if (psize > 1) {
-// return;
-// }
+ // return;
+ // }
#ifdef AKANTU_USE_SCOTCH
Int prank = this->communicator->whoAmI();
if (prank == 0) {
MeshPartitionScotch partition(*this, spatial_dimension);
partition.partitionate(psize);
MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition);
} else {
MeshUtilsDistribution::distributeMeshCentralized(*this, 0);
}
#else
if (!(psize == 1)) {
AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool");
}
#endif
this->computeBoundingBox();
}
this->is_distributed = true;
}
/* -------------------------------------------------------------------------- */
void Mesh::getAssociatedElements(const Array<UInt> & node_list,
Array<Element> & elements) {
for (const auto & node : node_list)
for (const auto & element : *nodes_to_elements[node])
elements.push_back(element);
}
/* -------------------------------------------------------------------------- */
void Mesh::fillNodesToElements() {
Element e;
UInt nb_nodes = nodes->size();
for (UInt n = 0; n < nb_nodes; ++n) {
if (this->nodes_to_elements[n])
this->nodes_to_elements[n]->clear();
else
this->nodes_to_elements[n] = std::make_unique<std::set<Element>>();
}
for (auto ghost_type : ghost_types) {
e.ghost_type = ghost_type;
for (const auto & type :
elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) {
e.type = type;
UInt nb_element = this->getNbElement(type, ghost_type);
Array<UInt>::const_iterator<Vector<UInt>> conn_it =
connectivities(type, ghost_type)
.begin(Mesh::getNbNodesPerElement(type));
for (UInt el = 0; el < nb_element; ++el, ++conn_it) {
e.element = el;
const Vector<UInt> & conn = *conn_it;
for (UInt n = 0; n < conn.size(); ++n)
nodes_to_elements[conn(n)]->insert(e);
}
}
}
}
/* -------------------------------------------------------------------------- */
std::tuple<UInt, UInt>
Mesh::updateGlobalData(NewNodesEvent & nodes_event,
NewElementsEvent & elements_event) {
if (global_data_updater)
return this->global_data_updater->updateData(nodes_event, elements_event);
else {
return std::make_tuple(nodes_event.getList().size(),
elements_event.getList().size());
}
}
/* -------------------------------------------------------------------------- */
void Mesh::registerGlobalDataUpdater(
std::unique_ptr<MeshGlobalDataUpdater> && global_data_updater) {
this->global_data_updater = std::move(global_data_updater);
}
/* -------------------------------------------------------------------------- */
void Mesh::eraseElements(const Array<Element> & elements) {
ElementTypeMap<UInt> last_element;
RemovedElementsEvent event(*this);
auto & remove_list = event.getList();
auto & new_numbering = event.getNewNumbering();
for (auto && el : elements) {
if (el.ghost_type != _not_ghost) {
auto & count = ghosts_counters(el);
--count;
if (count > 0)
continue;
}
remove_list.push_back(el);
if (not last_element.exists(el.type, el.ghost_type)) {
UInt nb_element = mesh.getNbElement(el.type, el.ghost_type);
last_element(nb_element - 1, el.type, el.ghost_type);
auto & numbering =
new_numbering.alloc(nb_element, 1, el.type, el.ghost_type);
for (auto && pair : enumerate(numbering)) {
std::get<1>(pair) = std::get<0>(pair);
}
}
UInt & pos = last_element(el.type, el.ghost_type);
auto & numbering = new_numbering(el.type, el.ghost_type);
numbering(el.element) = UInt(-1);
numbering(pos) = el.element;
--pos;
}
this->sendEvent(event);
}
} // namespace akantu
diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh
index a9e0b4a2a..d7b167358 100644
--- a/src/mesh/mesh.hh
+++ b/src/mesh/mesh.hh
@@ -1,622 +1,625 @@
/**
* @file mesh.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Dana Christen <dana.christen@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Mon Feb 19 2018
*
* @brief the class representing the meshes
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_HH__
#define __AKANTU_MESH_HH__
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "aka_event_handler_manager.hh"
#include "aka_memory.hh"
#include "dumpable.hh"
#include "element.hh"
#include "element_class.hh"
#include "element_type_map.hh"
#include "group_manager.hh"
#include "mesh_data.hh"
#include "mesh_events.hh"
/* -------------------------------------------------------------------------- */
#include <set>
/* -------------------------------------------------------------------------- */
namespace akantu {
class Communicator;
class ElementSynchronizer;
class NodeSynchronizer;
class MeshGlobalDataUpdater;
} // namespace akantu
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Mesh */
/* -------------------------------------------------------------------------- */
/**
* @class Mesh this contain the coordinates of the nodes in the Mesh.nodes
* Array, and the connectivity. The connectivity are stored in by element
* types.
*
* In order to loop on all element you have to loop on all types like this :
* @code{.cpp}
for(auto & type : mesh.elementTypes()) {
UInt nb_element = mesh.getNbElement(type);
const Array<UInt> & conn = mesh.getConnectivity(type);
for(UInt e = 0; e < nb_element; ++e) {
...
}
}
or
for_each_element(mesh, [](Element & element) {
std::cout << element << std::endl
});
@endcode
*/
class Mesh : protected Memory,
public EventHandlerManager<MeshEventHandler>,
public GroupManager,
public Dumpable {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
private:
/// default constructor used for chaining, the last parameter is just to
/// differentiate constructors
Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id,
Communicator & communicator);
public:
/// constructor that create nodes coordinates array
Mesh(UInt spatial_dimension, const ID & id = "mesh",
const MemoryID & memory_id = 0);
/// mesh not distributed and not using the default communicator
Mesh(UInt spatial_dimension, Communicator & communicator,
const ID & id = "mesh", const MemoryID & memory_id = 0);
/// constructor that use an existing nodes coordinates array, by knowing its
/// ID
// Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id,
// const MemoryID & memory_id = 0);
/**
* constructor that use an existing nodes coordinates
* array, by getting the vector of coordinates
*/
Mesh(UInt spatial_dimension, const std::shared_ptr<Array<Real>> & nodes,
const ID & id = "mesh", const MemoryID & memory_id = 0);
~Mesh() override;
/// read the mesh from a file
void read(const std::string & filename,
const MeshIOType & mesh_io_type = _miot_auto);
/// write the mesh to a file
void write(const std::string & filename,
const MeshIOType & mesh_io_type = _miot_auto);
private:
+
+ void makeReady();
+
/// initialize the connectivity to NULL and other stuff
void init();
/// function that computes the bounding box (fills xmin, xmax)
void computeBoundingBox();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// patitionate the mesh among the processors involved in their computation
virtual void distribute(Communicator & communicator);
virtual void distribute();
/// function to print the containt of the class
void printself(std::ostream & stream, int indent = 0) const override;
/// extract coordinates of nodes from an element
template <typename T>
inline void extractNodalValuesFromElement(const Array<T> & nodal_values,
T * elemental_values,
UInt * connectivity, UInt n_nodes,
UInt nb_degree_of_freedom) const;
// /// extract coordinates of nodes from a reversed element
// inline void extractNodalCoordinatesFromPBCElement(Real * local_coords,
// UInt * connectivity,
// UInt n_nodes);
/// add a Array of connectivity for the type <type>.
inline void addConnectivityType(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
template <class Event> inline void sendEvent(Event & event) {
// if(event.getList().size() != 0)
EventHandlerManager<MeshEventHandler>::sendEvent<Event>(event);
}
/// prepare the event to remove the elements listed
void eraseElements(const Array<Element> & elements);
/* ------------------------------------------------------------------------ */
template <typename T>
inline void removeNodesFromArray(Array<T> & vect,
const Array<UInt> & new_numbering);
/// initialize normals
void initNormals();
/// init facets' mesh
Mesh & initMeshFacets(const ID & id = "mesh_facets");
/// define parent mesh
void defineMeshParent(const Mesh & mesh);
/// get global connectivity array
void getGlobalConnectivity(ElementTypeMapArray<UInt> & global_connectivity);
public:
void getAssociatedElements(const Array<UInt> & node_list,
Array<Element> & elements);
private:
/// fills the nodes_to_elements structure
void fillNodesToElements();
/// update the global ids, nodes type, ...
std::tuple<UInt, UInt> updateGlobalData(NewNodesEvent & nodes_event,
NewElementsEvent & elements_event);
void registerGlobalDataUpdater(
std::unique_ptr<MeshGlobalDataUpdater> && global_data_updater);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the id of the mesh
AKANTU_GET_MACRO(ID, Memory::id, const ID &);
/// get the id of the mesh
AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &);
/// get the spatial dimension of the mesh = number of component of the
/// coordinates
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// get the nodes Array aka coordinates
AKANTU_GET_MACRO(Nodes, *nodes, const Array<Real> &);
AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array<Real> &);
/// get the normals for the elements
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real);
/// get the number of nodes
AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt);
/// get the Array of global ids of the nodes (only used in parallel)
AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array<UInt> &);
AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array<UInt> &);
/// get the global id of a node
inline UInt getNodeGlobalId(UInt local_id) const;
/// get the global id of a node
inline UInt getNodeLocalId(UInt global_id) const;
/// get the global number of nodes
inline UInt getNbGlobalNodes() const;
/// get the nodes type Array
AKANTU_GET_MACRO(NodesType, *nodes_type, const Array<NodeType> &);
protected:
AKANTU_GET_MACRO_NOT_CONST(NodesType, *nodes_type, Array<NodeType> &);
public:
inline NodeType getNodeType(UInt local_id) const;
/// say if a node is a pure ghost node
inline bool isPureGhostNode(UInt n) const;
/// say if a node is pur local or master node
inline bool isLocalOrMasterNode(UInt n) const;
inline bool isLocalNode(UInt n) const;
inline bool isMasterNode(UInt n) const;
inline bool isSlaveNode(UInt n) const;
AKANTU_GET_MACRO(LowerBounds, lower_bounds, const Vector<Real> &);
AKANTU_GET_MACRO(UpperBounds, upper_bounds, const Vector<Real> &);
AKANTU_GET_MACRO(LocalLowerBounds, local_lower_bounds, const Vector<Real> &);
AKANTU_GET_MACRO(LocalUpperBounds, local_upper_bounds, const Vector<Real> &);
/// get the connectivity Array for a given type
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt);
AKANTU_GET_MACRO(Connectivities, connectivities,
const ElementTypeMapArray<UInt> &);
/// get the number of element of a type in the mesh
inline UInt getNbElement(const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
/// get the number of element for a given ghost_type and a given dimension
inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions,
const GhostType & ghost_type = _not_ghost,
const ElementKind & kind = _ek_not_defined) const;
/// compute the barycenter of a given element
inline void getBarycenter(const Element & element,
Vector<Real> & barycenter) const;
void getBarycenters(Array<Real> & barycenter, const ElementType & type,
const GhostType & ghost_type) const;
#ifndef SWIG
/// get the element connected to a subelement (element of lower dimension)
const auto & getElementToSubelement() const;
/// get the element connected to a subelement
const auto &
getElementToSubelement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get the element connected to a subelement
auto & getElementToSubelement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost);
/// get the elements connected to a subelement
const auto & getElementToSubelement(const Element & element) const;
/// get the subelement (element of lower dimension) connected to a element
const auto & getSubelementToElement() const;
/// get the subelement connected to an element
const auto &
getSubelementToElement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get the subelement connected to an element
auto & getSubelementToElement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost);
/// get the subelement (element of lower dimension) connected to a element
VectorProxy<Element> getSubelementToElement(const Element & element) const;
/// get connectivity of a given element
inline VectorProxy<UInt> getConnectivity(const Element & element) const;
protected:
inline auto & getElementToSubelement(const Element & element);
inline VectorProxy<Element> getSubelementToElement(const Element & element);
inline VectorProxy<UInt> getConnectivity(const Element & element);
#endif
public:
/// register a new ElementalTypeMap in the MeshData
template <typename T>
inline ElementTypeMapArray<T> & registerData(const std::string & data_name);
template <typename T>
inline bool hasData(const ID & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get a name field associated to the mesh
template <typename T>
inline const Array<T> &
getData(const ID & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get a name field associated to the mesh
template <typename T>
inline Array<T> & getData(const ID & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost);
/// tells if the data data_name is contained in the mesh or not
inline bool hasData(const ID & data_name) const;
/// get a name field associated to the mesh
template <typename T>
inline const ElementTypeMapArray<T> & getData(const ID & data_name) const;
/// get a name field associated to the mesh
template <typename T>
inline ElementTypeMapArray<T> & getData(const ID & data_name);
template <typename T>
ElementTypeMap<UInt> getNbDataPerElem(ElementTypeMapArray<T> & array,
const ElementKind & element_kind);
template <typename T>
dumper::Field * createFieldFromAttachedData(const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind);
/// templated getter returning the pointer to data in MeshData (modifiable)
template <typename T>
inline Array<T> &
getDataPointer(const std::string & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost,
UInt nb_component = 1, bool size_to_nb_element = true,
bool resize_with_parent = false);
template <typename T>
inline Array<T> & getDataPointer(const ID & data_name,
const ElementType & el_type,
const GhostType & ghost_type,
UInt nb_component, bool size_to_nb_element,
bool resize_with_parent, const T & defaul_);
/// Facets mesh accessor
inline const Mesh & getMeshFacets() const;
inline Mesh & getMeshFacets();
/// Parent mesh accessor
inline const Mesh & getMeshParent() const;
inline bool isMeshFacets() const { return this->is_mesh_facets; }
/// defines is the mesh is distributed or not
inline bool isDistributed() const { return this->is_distributed; }
#ifndef SWIG
/// return the dumper from a group and and a dumper name
DumperIOHelper & getGroupDumper(const std::string & dumper_name,
const std::string & group_name);
#endif
/* ------------------------------------------------------------------------ */
/* Wrappers on ElementClass functions */
/* ------------------------------------------------------------------------ */
public:
/// get the number of nodes per element for a given element type
static inline UInt getNbNodesPerElement(const ElementType & type);
/// get the number of nodes per element for a given element type considered as
/// a first order element
static inline ElementType getP1ElementType(const ElementType & type);
/// get the kind of the element type
static inline ElementKind getKind(const ElementType & type);
/// get spatial dimension of a type of element
static inline UInt getSpatialDimension(const ElementType & type);
/// get number of facets of a given element type
static inline UInt getNbFacetsPerElement(const ElementType & type);
/// get number of facets of a given element type
static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t);
#ifndef SWIG
/// get local connectivity of a facet for a given facet type
static inline auto getFacetLocalConnectivity(const ElementType & type,
UInt t = 0);
/// get connectivity of facets for a given element
inline auto getFacetConnectivity(const Element & element, UInt t = 0) const;
#endif
/// get the number of type of the surface element associated to a given
/// element type
static inline UInt getNbFacetTypes(const ElementType & type, UInt t = 0);
#ifndef SWIG
/// get the type of the surface element associated to a given element
static inline constexpr auto getFacetType(const ElementType & type,
UInt t = 0);
/// get all the type of the surface element associated to a given element
static inline constexpr auto getAllFacetTypes(const ElementType & type);
#endif
/// get the number of nodes in the given element list
static inline UInt getNbNodesPerElementList(const Array<Element> & elements);
/* ------------------------------------------------------------------------ */
/* Element type Iterator */
/* ------------------------------------------------------------------------ */
using type_iterator = ElementTypeMapArray<UInt, ElementType>::type_iterator;
using ElementTypesIteratorHelper =
ElementTypeMapArray<UInt, ElementType>::ElementTypesIteratorHelper;
template <typename... pack>
ElementTypesIteratorHelper elementTypes(pack &&... _pack) const;
inline type_iterator firstType(UInt dim = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_regular) const {
return connectivities.firstType(dim, ghost_type, kind);
}
inline type_iterator lastType(UInt dim = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_regular) const {
return connectivities.lastType(dim, ghost_type, kind);
}
AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer,
const ElementSynchronizer &);
AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer,
ElementSynchronizer &);
AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer,
const NodeSynchronizer &);
AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer,
NodeSynchronizer &);
// AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator
// &);
#ifndef SWIG
AKANTU_GET_MACRO(Communicator, *communicator, const auto &);
AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &);
#endif
/* ------------------------------------------------------------------------ */
/* Private methods for friends */
/* ------------------------------------------------------------------------ */
private:
friend class MeshAccessor;
friend class MeshUtils;
AKANTU_GET_MACRO(NodesPointer, *nodes, Array<Real> &);
/// get a pointer to the nodes_global_ids Array<UInt> and create it if
/// necessary
inline Array<UInt> & getNodesGlobalIdsPointer();
/// get a pointer to the nodes_type Array<Int> and create it if necessary
inline Array<NodeType> & getNodesTypePointer();
/// get a pointer to the connectivity Array for the given type and create it
/// if necessary
inline Array<UInt> &
getConnectivityPointer(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
/// get the ghost element counter
inline Array<UInt> &
getGhostsCounters(const ElementType & type,
const GhostType & ghost_type = _ghost) {
AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost,
"No ghost counter for _not_ghost elements");
return ghosts_counters(type, ghost_type);
}
/// get a pointer to the element_to_subelement Array for the given type and
/// create it if necessary
inline Array<std::vector<Element>> &
getElementToSubelementPointer(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
/// get a pointer to the subelement_to_element Array for the given type and
/// create it if necessary
inline Array<Element> &
getSubelementToElementPointer(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
AKANTU_GET_MACRO_NOT_CONST(MeshData, mesh_data, MeshData &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// array of the nodes coordinates
std::shared_ptr<Array<Real>> nodes;
/// global node ids
std::shared_ptr<Array<UInt>> nodes_global_ids;
/// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in
/// [0-N] slave node and master is proc i
std::shared_ptr<Array<NodeType>> nodes_type;
/// global number of nodes;
UInt nb_global_nodes{0};
/// all class of elements present in this mesh (for heterogenous meshes)
ElementTypeMapArray<UInt> connectivities;
/// count the references on ghost elements
ElementTypeMapArray<UInt> ghosts_counters;
/// map to normals for all class of elements present in this mesh
ElementTypeMapArray<Real> normals;
/// the spatial dimension of this mesh
UInt spatial_dimension{0};
/// min of coordinates
Vector<Real> lower_bounds;
/// max of coordinates
Vector<Real> upper_bounds;
/// size covered by the mesh on each direction
Vector<Real> size;
/// local min of coordinates
Vector<Real> local_lower_bounds;
/// local max of coordinates
Vector<Real> local_upper_bounds;
/// Extra data loaded from the mesh file
MeshData mesh_data;
/// facets' mesh
std::unique_ptr<Mesh> mesh_facets;
/// parent mesh (this is set for mesh_facets meshes)
const Mesh * mesh_parent{nullptr};
/// defines if current mesh is mesh_facets or not
bool is_mesh_facets{false};
/// defines if the mesh is centralized or distributed
bool is_distributed{false};
/// Communicator on which mesh is distributed
Communicator * communicator;
/// Element synchronizer
std::unique_ptr<ElementSynchronizer> element_synchronizer;
/// Node synchronizer
std::unique_ptr<NodeSynchronizer> node_synchronizer;
using NodesToElements = std::vector<std::unique_ptr<std::set<Element>>>;
/// class to update global data using external knowledge
std::unique_ptr<MeshGlobalDataUpdater> global_data_updater;
/// This info is stored to simplify the dynamic changes
NodesToElements nodes_to_elements;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* Inline functions */
/* -------------------------------------------------------------------------- */
#include "element_type_map_tmpl.hh"
#include "mesh_inline_impl.cc"
#endif /* __AKANTU_MESH_HH__ */
diff --git a/src/mesh/mesh_accessor.hh b/src/mesh/mesh_accessor.hh
index c841b5f1b..39abd42a1 100644
--- a/src/mesh/mesh_accessor.hh
+++ b/src/mesh/mesh_accessor.hh
@@ -1,150 +1,152 @@
/**
* @file mesh_accessor.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jun 30 2015
* @date last modification: Tue Sep 19 2017
*
* @brief this class allow to access some private member of mesh it is used for
* IO for examples
*
* @section LICENSE
*
* Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_ACCESSOR_HH__
#define __AKANTU_MESH_ACCESSOR_HH__
namespace akantu {
class NodeSynchronizer;
class ElementSynchronizer;
class MeshGlobalDataUpdater;
} // namespace akantu
namespace akantu {
class MeshAccessor {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
explicit MeshAccessor(Mesh & mesh) : _mesh(mesh) {}
virtual ~MeshAccessor() = default;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the global number of nodes
inline UInt getNbGlobalNodes() const { return this->_mesh.nb_global_nodes; }
/// set the global number of nodes
inline void setNbGlobalNodes(UInt nb_global_nodes) {
this->_mesh.nb_global_nodes = nb_global_nodes;
}
/// set the mesh as being distributed
inline void setDistributed() { this->_mesh.is_distributed = true; }
/// get a pointer to the nodes_global_ids Array<UInt> and create it if
/// necessary
inline auto & getNodesGlobalIds() {
return this->_mesh.getNodesGlobalIdsPointer();
}
/// get a pointer to the nodes_type Array<Int> and create it if necessary
inline auto & getNodesType() { return this->_mesh.getNodesTypePointer(); }
/// get a pointer to the coordinates Array
inline auto & getNodes() { return this->_mesh.getNodesPointer(); }
/// get a pointer to the coordinates Array
inline auto getNodesSharedPtr() { return this->_mesh.nodes; }
/// get a pointer to the connectivity Array for the given type and create it
/// if necessary
inline auto & getConnectivity(const ElementType & type,
const GhostType & ghost_type = _not_ghost) {
return this->_mesh.getConnectivityPointer(type, ghost_type);
}
/// get the ghost element counter
inline auto & getGhostsCounters(const ElementType & type,
const GhostType & ghost_type = _ghost) {
return this->_mesh.getGhostsCounters(type, ghost_type);
}
/// get a pointer to the element_to_subelement Array for the given type and
/// create it if necessary
inline auto &
getElementToSubelement(const ElementType & type,
const GhostType & ghost_type = _not_ghost) {
return this->_mesh.getElementToSubelementPointer(type, ghost_type);
}
/// get a pointer to the subelement_to_element Array for the given type and
/// create it if necessary
inline auto &
getSubelementToElement(const ElementType & type,
const GhostType & ghost_type = _not_ghost) {
return this->_mesh.getSubelementToElementPointer(type, ghost_type);
}
template <typename T>
inline auto &
getData(const std::string & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost, UInt nb_component = 1,
bool size_to_nb_element = true, bool resize_with_parent = false) {
return this->_mesh.getDataPointer<T>(data_name, el_type, ghost_type,
nb_component, size_to_nb_element,
resize_with_parent);
}
auto & getMeshData() { return this->_mesh.getMeshData(); }
/// get the node synchonizer
auto & getNodeSynchronizer() { return *this->_mesh.node_synchronizer; }
/// get the element synchonizer
auto & getElementSynchronizer() { return *this->_mesh.element_synchronizer; }
decltype(auto) updateGlobalData(NewNodesEvent & nodes_event,
NewElementsEvent & elements_event) {
return this->_mesh.updateGlobalData(nodes_event, elements_event);
}
void registerGlobalDataUpdater(
std::unique_ptr<MeshGlobalDataUpdater> && global_data_updater) {
this->_mesh.registerGlobalDataUpdater(
std::forward<std::unique_ptr<MeshGlobalDataUpdater>>(
global_data_updater));
}
+ void makeReady() { return this->_mesh.makeReady(); }
+
private:
Mesh & _mesh;
};
} // namespace akantu
#endif /* __AKANTU_MESH_ACCESSOR_HH__ */
diff --git a/src/model/boundary_condition_tmpl.hh b/src/model/boundary_condition_tmpl.hh
index f45d2ac88..4b5e6399f 100644
--- a/src/model/boundary_condition_tmpl.hh
+++ b/src/model/boundary_condition_tmpl.hh
@@ -1,256 +1,256 @@
/**
* @file boundary_condition_tmpl.hh
*
* @author Dana Christen <dana.christen@gmail.com>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri May 03 2013
* @date last modification: Tue Feb 20 2018
*
* @brief implementation of the applyBC
*
* @section LICENSE
*
* Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "element_group.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_BOUNDARY_CONDITION_TMPL_HH__
#define __AKANTU_BOUNDARY_CONDITION_TMPL_HH__
namespace akantu {
/* -------------------------------------------------------------------------- */
template <typename ModelType>
void BoundaryCondition<ModelType>::initBC(ModelType & model,
Array<Real> & primal,
Array<Real> & dual) {
this->model = &model;
this->primal = &primal;
this->dual = &dual;
}
/* -------------------------------------------------------------------------- */
template <typename ModelType>
void BoundaryCondition<ModelType>::initBC(ModelType & model,
Array<Real> & primal,
Array<Real> & primal_increment,
Array<Real> & dual) {
this->initBC(model, primal, dual);
this->primal_increment = &primal_increment;
}
/* -------------------------------------------------------------------------- */
/* Partial specialization for DIRICHLET functors */
template <typename ModelType>
template <typename FunctorType>
struct BoundaryCondition<ModelType>::TemplateFunctionWrapper<
FunctorType, BC::Functor::_dirichlet> {
static inline void applyBC(const FunctorType & func,
const ElementGroup & group,
BoundaryCondition<ModelType> & bc_instance) {
auto & model = bc_instance.getModel();
auto & primal = bc_instance.getPrimal();
const auto & coords = model.getMesh().getNodes();
auto & boundary_flags = model.getBlockedDOFs();
UInt dim = model.getMesh().getSpatialDimension();
auto primal_iter = primal.begin(primal.getNbComponent());
auto coords_iter = coords.begin(dim);
auto flags_iter = boundary_flags.begin(boundary_flags.getNbComponent());
for (auto n : group.getNodeGroup()) {
Vector<bool> flag(flags_iter[n]);
Vector<Real> primal(primal_iter[n]);
Vector<Real> coords(coords_iter[n]);
func(n, flag, primal, coords);
}
}
};
/* -------------------------------------------------------------------------- */
/* Partial specialization for NEUMANN functors */
template <typename ModelType>
template <typename FunctorType>
struct BoundaryCondition<ModelType>::TemplateFunctionWrapper<
FunctorType, BC::Functor::_neumann> {
static inline void applyBC(const FunctorType & func,
const ElementGroup & group,
BoundaryCondition<ModelType> & bc_instance) {
UInt dim = bc_instance.getModel().getSpatialDimension();
switch (dim) {
case 1: {
AKANTU_TO_IMPLEMENT();
break;
}
case 2:
case 3: {
applyBC(func, group, bc_instance, _not_ghost);
applyBC(func, group, bc_instance, _ghost);
break;
}
}
}
static inline void applyBC(const FunctorType & func,
const ElementGroup & group,
BoundaryCondition<ModelType> & bc_instance,
GhostType ghost_type) {
auto & model = bc_instance.getModel();
auto & dual = bc_instance.getDual();
const auto & mesh = model.getMesh();
const auto & nodes_coords = mesh.getNodes();
const auto & fem_boundary = model.getFEEngineBoundary();
UInt dim = model.getSpatialDimension();
UInt nb_degree_of_freedom = dual.getNbComponent();
IntegrationPoint quad_point;
quad_point.ghost_type = ghost_type;
// Loop over the boundary element types
for (auto && type : group.elementTypes(dim - 1, ghost_type)) {
const auto & element_ids = group.getElements(type, ghost_type);
UInt nb_quad_points =
fem_boundary.getNbIntegrationPoints(type, ghost_type);
UInt nb_elements = element_ids.size();
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
Array<Real> * dual_before_integ = new Array<Real>(
nb_elements * nb_quad_points, nb_degree_of_freedom, 0.);
Array<Real> * quad_coords =
new Array<Real>(nb_elements * nb_quad_points, dim);
const auto & normals_on_quad =
fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type);
fem_boundary.interpolateOnIntegrationPoints(
nodes_coords, *quad_coords, dim, type, ghost_type, element_ids);
auto normals_begin = normals_on_quad.begin(dim);
decltype(normals_begin) normals_iter;
auto quad_coords_iter = quad_coords->begin(dim);
auto dual_iter = dual_before_integ->begin(nb_degree_of_freedom);
quad_point.type = type;
for (auto el : element_ids) {
quad_point.element = el;
normals_iter = normals_begin + el * nb_quad_points;
for (auto && q : arange(nb_quad_points)) {
quad_point.num_point = q;
func(quad_point, *dual_iter, *quad_coords_iter, *normals_iter);
++dual_iter;
++quad_coords_iter;
++normals_iter;
}
}
delete quad_coords;
/* -------------------------------------------------------------------- */
// Initialization of iterators
auto dual_iter_mat = dual_before_integ->begin(nb_degree_of_freedom, 1);
auto shapes_iter_begin = fem_boundary.getShapes(type, ghost_type)
.begin(1, nb_nodes_per_element);
Array<Real> * dual_by_shapes =
new Array<Real>(nb_elements * nb_quad_points,
nb_degree_of_freedom * nb_nodes_per_element);
auto dual_by_shapes_iter =
dual_by_shapes->begin(nb_degree_of_freedom, nb_nodes_per_element);
/* -------------------------------------------------------------------- */
// Loop computing dual x shapes
for (auto el : element_ids) {
auto shapes_iter = shapes_iter_begin + el * nb_quad_points;
for (UInt q(0); q < nb_quad_points;
++q, ++dual_iter_mat, ++dual_by_shapes_iter, ++shapes_iter) {
dual_by_shapes_iter->template mul<false, false>(*dual_iter_mat,
*shapes_iter);
}
}
delete dual_before_integ;
Array<Real> * dual_by_shapes_integ = new Array<Real>(
nb_elements, nb_degree_of_freedom * nb_nodes_per_element);
fem_boundary.integrate(*dual_by_shapes, *dual_by_shapes_integ,
nb_degree_of_freedom * nb_nodes_per_element, type,
ghost_type, element_ids);
delete dual_by_shapes;
// assemble the result into force vector
model.getDOFManager().assembleElementalArrayLocalArray(
*dual_by_shapes_integ, dual, type, ghost_type, 1., element_ids);
delete dual_by_shapes_integ;
}
}
};
/* -------------------------------------------------------------------------- */
template <typename ModelType>
template <typename FunctorType>
inline void BoundaryCondition<ModelType>::applyBC(const FunctorType & func) {
auto bit = model->getMesh().getGroupManager().element_group_begin();
auto bend = model->getMesh().getGroupManager().element_group_end();
for (; bit != bend; ++bit)
applyBC(func, *bit);
}
/* -------------------------------------------------------------------------- */
template <typename ModelType>
template <typename FunctorType>
inline void
BoundaryCondition<ModelType>::applyBC(const FunctorType & func,
const std::string & group_name) {
try {
const ElementGroup & element_group =
model->getMesh().getElementGroup(group_name);
applyBC(func, element_group);
- } catch (akantu::debug::Exception e) {
+ } catch (akantu::debug::Exception & e) {
AKANTU_EXCEPTION("Error applying a boundary condition onto \""
<< group_name << "\"! [" << e.what() << "]");
}
}
/* -------------------------------------------------------------------------- */
template <typename ModelType>
template <typename FunctorType>
inline void
BoundaryCondition<ModelType>::applyBC(const FunctorType & func,
const ElementGroup & element_group) {
#if !defined(AKANTU_NDEBUG)
if (element_group.getDimension() != model->getSpatialDimension() - 1)
AKANTU_DEBUG_WARNING("The group "
<< element_group.getName()
<< " does not contain only boundaries elements");
#endif
TemplateFunctionWrapper<FunctorType>::applyBC(func, element_group, *this);
}
#endif /* __AKANTU_BOUNDARY_CONDITION_TMPL_HH__ */
} // namespace akantu
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
index 5c4404dbd..246c372c8 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
@@ -1,274 +1,274 @@
/**
* @file material_elastic_linear_anisotropic.cc
*
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Enrico Milanese <enrico.milanese@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Sep 25 2013
* @date last modification: Tue Feb 20 2018
*
* @brief Anisotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "material_elastic_linear_anisotropic.hh"
#include "solid_mechanics_model.hh"
#include <algorithm>
#include <sstream>
/* -------------------------------------------------------------------------- */
#define AKANTU_WARNING_UNDEFINED_VAR_TEMPLATE
#include "aka_warning.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template <UInt dim>
MaterialElasticLinearAnisotropic<dim>::MaterialElasticLinearAnisotropic(
SolidMechanicsModel & model, const ID & id, bool symmetric)
- : Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim),
- C(voigt_h::size, voigt_h::size), eigC(voigt_h::size),
+ : Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim, 0.),
+ C(voigt_h::size, voigt_h::size, 0.), eigC(voigt_h::size, 0.),
symmetric(symmetric), alpha(0), was_stiffness_assembled(false) {
AKANTU_DEBUG_IN();
- this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim));
+ this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim, 0.));
(*this->dir_vecs.back())[0] = 1.;
this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of main material axis");
if (dim > 1) {
- this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim));
+ this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim, 0.));
(*this->dir_vecs.back())[1] = 1.;
this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of secondary material axis");
}
if (dim > 2) {
- this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim));
+ this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim, 0.));
(*this->dir_vecs.back())[2] = 1.;
this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of tertiary material axis");
}
for (UInt i = 0; i < voigt_h::size; ++i) {
UInt start = 0;
if (this->symmetric) {
start = i;
}
for (UInt j = start; j < voigt_h::size; ++j) {
std::stringstream param("C");
param << "C" << i + 1 << j + 1;
this->registerParam(param.str(), this->Cprime(i, j), Real(0.),
_pat_parsmod, "Coefficient " + param.str());
}
}
this->registerParam("alpha", this->alpha, _pat_parsmod,
"Proportion of viscous stress");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim> void MaterialElasticLinearAnisotropic<dim>::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MaterialElasticLinearAnisotropic<dim>::updateInternalParameters() {
Material::updateInternalParameters();
if (this->symmetric) {
for (UInt i = 0; i < voigt_h::size; ++i) {
for (UInt j = i + 1; j < voigt_h::size; ++j) {
this->Cprime(j, i) = this->Cprime(i, j);
}
}
}
this->rotateCprime();
this->C.eig(this->eigC);
this->was_stiffness_assembled = false;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void MaterialElasticLinearAnisotropic<Dim>::rotateCprime() {
// start by filling the empty parts fo Cprime
UInt diff = Dim * Dim - voigt_h::size;
for (UInt i = voigt_h::size; i < Dim * Dim; ++i) {
for (UInt j = 0; j < Dim * Dim; ++j) {
this->Cprime(i, j) = this->Cprime(i - diff, j);
}
}
for (UInt i = 0; i < Dim * Dim; ++i) {
for (UInt j = voigt_h::size; j < Dim * Dim; ++j) {
this->Cprime(i, j) = this->Cprime(i, j - diff);
}
}
// construction of rotator tensor
// normalise rotation matrix
for (UInt j = 0; j < Dim; ++j) {
Vector<Real> rot_vec = this->rot_mat(j);
rot_vec = *this->dir_vecs[j];
rot_vec.normalize();
}
// make sure the vectors form a right-handed base
Vector<Real> test_axis(3);
Vector<Real> v1(3), v2(3), v3(3);
if (Dim == 2) {
for (UInt i = 0; i < Dim; ++i) {
v1[i] = this->rot_mat(0, i);
v2[i] = this->rot_mat(1, i);
v3[i] = 0.;
}
v3[2] = 1.;
v1[2] = 0.;
v2[2] = 0.;
} else if (Dim == 3) {
v1 = this->rot_mat(0);
v2 = this->rot_mat(1);
v3 = this->rot_mat(2);
}
test_axis.crossProduct(v1, v2);
test_axis -= v3;
if (test_axis.norm() > 8 * std::numeric_limits<Real>::epsilon()) {
AKANTU_ERROR("The axis vectors do not form a right-handed coordinate "
<< "system. I. e., ||n1 x n2 - n3|| should be zero, but "
<< "it is " << test_axis.norm() << ".");
}
// create the rotator and the reverse rotator
Matrix<Real> rotator(Dim * Dim, Dim * Dim);
Matrix<Real> revrotor(Dim * Dim, Dim * Dim);
for (UInt i = 0; i < Dim; ++i) {
for (UInt j = 0; j < Dim; ++j) {
for (UInt k = 0; k < Dim; ++k) {
for (UInt l = 0; l < Dim; ++l) {
UInt I = voigt_h::mat[i][j];
UInt J = voigt_h::mat[k][l];
rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j);
revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l);
}
}
}
}
// create the full rotated matrix
Matrix<Real> Cfull(Dim * Dim, Dim * Dim);
Cfull = rotator * Cprime * revrotor;
for (UInt i = 0; i < voigt_h::size; ++i) {
for (UInt j = 0; j < voigt_h::size; ++j) {
this->C(i, j) = Cfull(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MaterialElasticLinearAnisotropic<dim>::computeStress(
ElementType el_type, GhostType ghost_type) {
// Wikipedia convention:
// 2*eps_ij (i!=j) = voigt_eps_I
// http://en.wikipedia.org/wiki/Voigt_notation
AKANTU_DEBUG_IN();
Matrix<Real> strain(dim, dim);
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
this->computeStressOnQuad(strain, sigma);
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MaterialElasticLinearAnisotropic<dim>::computeTangentModuli(
const ElementType & el_type, Array<Real> & tangent_matrix,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
this->computeTangentModuliOnQuad(tangent);
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
this->was_stiffness_assembled = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MaterialElasticLinearAnisotropic<dim>::computePotentialEnergy(
ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Material::computePotentialEnergy(el_type, ghost_type);
AKANTU_DEBUG_ASSERT(!this->finite_deformation,
"finite deformation not possible in material anisotropic "
"(TO BE IMPLEMENTED)");
if (ghost_type != _not_ghost)
return;
Array<Real>::scalar_iterator epot =
this->potential_energy(el_type, ghost_type).begin();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
computePotentialEnergyOnQuad(grad_u, sigma, *epot);
++epot;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real MaterialElasticLinearAnisotropic<dim>::getCelerity(
__attribute__((unused)) const Element & element) const {
return std::sqrt(this->eigC(0) / rho);
}
/* -------------------------------------------------------------------------- */
INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic);
-} // akantu
+} // namespace akantu
/* -------------------------------------------------------------------------- */
#include "aka_warning_restore.hh"
/* -------------------------------------------------------------------------- */
diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh b/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh
index 0cea87042..819f3f568 100644
--- a/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh
+++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh
@@ -1,136 +1,136 @@
/**
* @file material_elastic_orthotropic.hh
*
* @author Till Junge <till.junge@epfl.ch>
* @author Enrico Milanese <enrico.milanese@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Feb 16 2018
*
* @brief Orthotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material_elastic_linear_anisotropic.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__
#define __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__
namespace akantu {
/**
* Orthotropic elastic material
*
* parameters in the material files :
* - n1 : direction of x-axis in material base, normalisation not necessary
* (default: {1, 0, 0})
* - n2 : direction of y-axis in material base, normalisation not necessary
* (default: {0, 1, 0})
* - n3 : direction of z-axis in material base, normalisation not necessary
* (default: {0, 0, 1})
* - rho : density (default: 0)
* - E1 : Young's modulus along n1 (default: 0)
* - E2 : Young's modulus along n2 (default: 0)
* - E3 : Young's modulus along n3 (default: 0)
* - nu12 : Poisson's ratio along 12 (default: 0)
* - nu13 : Poisson's ratio along 13 (default: 0)
* - nu23 : Poisson's ratio along 23 (default: 0)
* - G12 : Shear modulus along 12 (default: 0)
* - G13 : Shear modulus along 13 (default: 0)
* - G23 : Shear modulus along 23 (default: 0)
*/
template <UInt Dim>
class MaterialElasticOrthotropic
: public MaterialElasticLinearAnisotropic<Dim> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElasticOrthotropic(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial() override;
void updateInternalParameters() override;
void
computePotentialEnergyByElement(ElementType type, UInt index,
Vector<Real> & epot_on_quad_points) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(E1, E1, Real);
AKANTU_GET_MACRO(E2, E2, Real);
AKANTU_GET_MACRO(E3, E3, Real);
AKANTU_GET_MACRO(Nu12, nu12, Real);
AKANTU_GET_MACRO(Nu13, nu13, Real);
AKANTU_GET_MACRO(Nu23, nu23, Real);
AKANTU_GET_MACRO(G12, G12, Real);
AKANTU_GET_MACRO(G13, G13, Real);
AKANTU_GET_MACRO(G23, G23, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// the n1 young modulus
- Real E1;
+ Real E1{0};
/// the n2 young modulus
- Real E2;
+ Real E2{0};
/// the n3 young modulus
- Real E3;
+ Real E3{0};
/// 12 Poisson coefficient
- Real nu12;
+ Real nu12{0};
/// 13 Poisson coefficient
- Real nu13;
+ Real nu13{0};
/// 23 Poisson coefficient
- Real nu23;
+ Real nu23{0};
/// 12 shear modulus
- Real G12;
+ Real G12{0};
/// 13 shear modulus
- Real G13;
+ Real G13{0};
/// 23 shear modulus
- Real G23;
+ Real G23{0};
};
-} // akantu
+} // namespace akantu
#endif /* __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__ */
diff --git a/src/synchronizer/communicator.cc b/src/synchronizer/communicator.cc
index fc925697c..7d8685551 100644
--- a/src/synchronizer/communicator.cc
+++ b/src/synchronizer/communicator.cc
@@ -1,181 +1,179 @@
/**
* @file communicator.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Mon Feb 05 2018
*
* @brief implementation of the common part of the static communicator
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
#if defined(AKANTU_USE_MPI)
#include "mpi_communicator_data.hh"
#endif
/* -------------------------------------------------------------------------- */
namespace akantu {
#if defined(AKANTU_USE_MPI)
int MPICommunicatorData::is_externaly_initialized = 0;
#endif
UInt InternalCommunicationRequest::counter = 0;
/* -------------------------------------------------------------------------- */
InternalCommunicationRequest::InternalCommunicationRequest(UInt source,
UInt dest)
: source(source), destination(dest) {
this->id = counter++;
}
/* -------------------------------------------------------------------------- */
InternalCommunicationRequest::~InternalCommunicationRequest() = default;
/* -------------------------------------------------------------------------- */
void InternalCommunicationRequest::printself(std::ostream & stream,
int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "CommunicationRequest [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + source : " << source << std::endl;
stream << space << " + destination : " << destination << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
Communicator::~Communicator() {
auto * event = new FinalizeCommunicatorEvent(*this);
this->sendEvent(*event);
delete event;
}
/* -------------------------------------------------------------------------- */
Communicator & Communicator::getStaticCommunicator() {
AKANTU_DEBUG_IN();
if (!static_communicator) {
int nb_args = 0;
char ** null = nullptr;
- static_communicator =
- std::make_unique<Communicator>(nb_args, null, private_member{});
+ static_communicator = std::make_unique<Communicator>(nb_args, null);
}
AKANTU_DEBUG_OUT();
return *static_communicator;
}
/* -------------------------------------------------------------------------- */
Communicator & Communicator::getStaticCommunicator(int & argc, char **& argv) {
if (!static_communicator)
- static_communicator =
- std::make_unique<Communicator>(argc, argv, private_member{});
+ static_communicator = std::make_unique<Communicator>(argc, argv);
return getStaticCommunicator();
}
} // namespace akantu
#ifdef AKANTU_USE_MPI
#include "communicator_mpi_inline_impl.cc"
#else
#include "communicator_dummy_inline_impl.cc"
#endif
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Template instantiation */
/* -------------------------------------------------------------------------- */
#define AKANTU_COMM_INSTANTIATE(T) \
template void Communicator::probe<T>(Int sender, Int tag, \
CommunicationStatus & status) const; \
template bool Communicator::asyncProbe<T>( \
Int sender, Int tag, CommunicationStatus & status) const; \
template void Communicator::sendImpl<T>( \
const T * buffer, Int size, Int receiver, Int tag, \
const CommunicationMode & mode) const; \
template void Communicator::receiveImpl<T>(T * buffer, Int size, Int sender, \
Int tag) const; \
template CommunicationRequest Communicator::asyncSendImpl<T>( \
const T * buffer, Int size, Int receiver, Int tag, \
const CommunicationMode & mode) const; \
template CommunicationRequest Communicator::asyncReceiveImpl<T>( \
T * buffer, Int size, Int sender, Int tag) const; \
template void Communicator::allGatherImpl<T>(T * values, int nb_values) \
const; \
template void Communicator::allGatherVImpl<T>(T * values, int * nb_values) \
const; \
template void Communicator::gatherImpl<T>(T * values, int nb_values, \
int root) const; \
template void Communicator::gatherImpl<T>( \
T * values, int nb_values, T * gathered, int nb_gathered) const; \
template void Communicator::gatherVImpl<T>(T * values, int * nb_values, \
int root) const; \
template void Communicator::broadcastImpl<T>(T * values, int nb_values, \
int root) const; \
template void Communicator::allReduceImpl<T>( \
T * values, int nb_values, const SynchronizerOperation & op) const
#define MIN_MAX_REAL SCMinMaxLoc<Real, int>
AKANTU_COMM_INSTANTIATE(bool);
AKANTU_COMM_INSTANTIATE(Real);
AKANTU_COMM_INSTANTIATE(UInt);
AKANTU_COMM_INSTANTIATE(Int);
AKANTU_COMM_INSTANTIATE(char);
AKANTU_COMM_INSTANTIATE(NodeType);
AKANTU_COMM_INSTANTIATE(MIN_MAX_REAL);
#if AKANTU_INTEGER_SIZE > 4
AKANTU_COMM_INSTANTIATE(int);
#endif
// template void Communicator::send<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * buffer, Int size, Int receiver, Int tag);
// template void Communicator::receive<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * buffer, Int size, Int sender, Int tag);
// template CommunicationRequest
// Communicator::asyncSend<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * buffer, Int size, Int receiver, Int tag);
// template CommunicationRequest
// Communicator::asyncReceive<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * buffer, Int size, Int sender, Int tag);
// template void Communicator::probe<SCMinMaxLoc<Real, int>>(
// Int sender, Int tag, CommunicationStatus & status);
// template void Communicator::allGather<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * values, int nb_values);
// template void Communicator::allGatherV<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * values, int * nb_values);
// template void Communicator::gather<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * values, int nb_values, int root);
// template void Communicator::gatherV<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * values, int * nb_values, int root);
// template void Communicator::broadcast<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * values, int nb_values, int root);
// template void Communicator::allReduce<SCMinMaxLoc<Real, int>>(
// SCMinMaxLoc<Real, int> * values, int nb_values,
// const SynchronizerOperation & op);
-} // akantu
+} // namespace akantu
diff --git a/src/synchronizer/communicator.hh b/src/synchronizer/communicator.hh
index e3ac8f462..1d53440a2 100644
--- a/src/synchronizer/communicator.hh
+++ b/src/synchronizer/communicator.hh
@@ -1,436 +1,438 @@
/**
* @file communicator.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Nov 15 2017
*
* @brief Class handling the parallel communications
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "aka_common.hh"
#include "aka_event_handler_manager.hh"
#include "communication_buffer.hh"
#include "communication_request.hh"
#include "communicator_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_STATIC_COMMUNICATOR_HH__
#define __AKANTU_STATIC_COMMUNICATOR_HH__
namespace akantu {
namespace debug {
class CommunicationException : public Exception {
public:
CommunicationException()
: Exception("An exception happen during a communication process.") {}
};
} // namespace debug
/// @enum SynchronizerOperation reduce operation that the synchronizer can
/// perform
enum class SynchronizerOperation {
_sum,
_min,
_max,
_prod,
_land,
_band,
_lor,
_bor,
_lxor,
_bxor,
_min_loc,
_max_loc,
_null
};
enum class CommunicationMode { _auto, _synchronous, _ready };
namespace {
int _any_source = -1;
}
} // namespace akantu
namespace akantu {
struct CommunicatorInternalData {
virtual ~CommunicatorInternalData() = default;
};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
class Communicator : public EventHandlerManager<CommunicatorEventHandler> {
- struct private_member {};
+ struct private_member;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- Communicator(int & argc, char **& argv, const private_member &);
+ Communicator(int & argc, char **& argv);
+ Communicator(const private_member &);
+
~Communicator() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Point to Point */
/* ------------------------------------------------------------------------ */
template <typename T>
inline void receive(Array<T> & values, Int sender, Int tag) const {
return this->receiveImpl(
values.storage(), values.size() * values.getNbComponent(), sender, tag);
}
template <typename Tensor>
inline void
receive(Tensor & values, Int sender, Int tag,
std::enable_if_t<is_tensor<Tensor>::value> * = nullptr) const {
return this->receiveImpl(values.storage(), values.size(), sender, tag);
}
template <bool is_static>
inline void receive(CommunicationBufferTemplated<is_static> & values,
Int sender, Int tag) const {
return this->receiveImpl(values.storage(), values.size(), sender, tag);
}
template <typename T>
inline void
receive(T & values, Int sender, Int tag,
std::enable_if_t<std::is_arithmetic<T>::value> * = nullptr) const {
return this->receiveImpl(&values, 1, sender, tag);
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void
send(const Array<T> & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const {
return this->sendImpl(values.storage(),
values.size() * values.getNbComponent(), receiver,
tag, mode);
}
template <typename Tensor>
inline void
send(const Tensor & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto,
std::enable_if_t<is_tensor<Tensor>::value> * = nullptr) const {
return this->sendImpl(values.storage(), values.size(), receiver, tag, mode);
}
template <bool is_static>
inline void
send(const CommunicationBufferTemplated<is_static> & values, Int receiver,
Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const {
return this->sendImpl(values.storage(), values.size(), receiver, tag, mode);
}
template <typename T>
inline void
send(const T & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto,
std::enable_if_t<std::is_arithmetic<T>::value> * = nullptr) const {
return this->sendImpl(&values, 1, receiver, tag, mode);
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline CommunicationRequest
asyncSend(const Array<T> & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const {
return this->asyncSendImpl(values.storage(),
values.size() * values.getNbComponent(),
receiver, tag, mode);
}
template <typename T>
inline CommunicationRequest
asyncSend(const std::vector<T> & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const {
return this->asyncSendImpl(values.data(), values.size(), receiver, tag,
mode);
}
template <typename Tensor>
inline CommunicationRequest
asyncSend(const Tensor & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto,
std::enable_if_t<is_tensor<Tensor>::value> * = nullptr) const {
return this->asyncSendImpl(values.storage(), values.size(), receiver, tag,
mode);
}
template <bool is_static>
inline CommunicationRequest
asyncSend(const CommunicationBufferTemplated<is_static> & values,
Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const {
return this->asyncSendImpl(values.storage(), values.size(), receiver, tag,
mode);
}
template <typename T>
inline CommunicationRequest
asyncSend(const T & values, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto,
std::enable_if_t<std::is_arithmetic<T>::value> * = nullptr) const {
return this->asyncSendImpl(&values, 1, receiver, tag, mode);
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline CommunicationRequest asyncReceive(Array<T> & values, Int sender,
Int tag) const {
return this->asyncReceiveImpl(
values.storage(), values.size() * values.getNbComponent(), sender, tag);
}
template <typename T>
inline CommunicationRequest asyncReceive(std::vector<T> & values, Int sender,
Int tag) const {
return this->asyncReceiveImpl(values.data(), values.size(), sender, tag);
}
template <typename Tensor,
typename = std::enable_if_t<is_tensor<Tensor>::value>>
inline CommunicationRequest asyncReceive(Tensor & values, Int sender,
Int tag) const {
return this->asyncReceiveImpl(values.storage(), values.size(), sender, tag);
}
template <bool is_static>
inline CommunicationRequest
asyncReceive(CommunicationBufferTemplated<is_static> & values, Int sender,
Int tag) const {
return this->asyncReceiveImpl(values.storage(), values.size(), sender, tag);
}
/* ------------------------------------------------------------------------ */
template <typename T>
void probe(Int sender, Int tag, CommunicationStatus & status) const;
template <typename T>
bool asyncProbe(Int sender, Int tag, CommunicationStatus & status) const;
/* ------------------------------------------------------------------------ */
/* Collectives */
/* ------------------------------------------------------------------------ */
template <typename T>
inline void allReduce(Array<T> & values,
const SynchronizerOperation & op) const {
this->allReduceImpl(values.storage(),
values.size() * values.getNbComponent(), op);
}
template <typename Tensor>
inline void
allReduce(Tensor & values, const SynchronizerOperation & op,
std::enable_if_t<is_tensor<Tensor>::value> * = nullptr) const {
this->allReduceImpl(values.storage(), values.size(), op);
}
template <typename T>
inline void
allReduce(T & values, const SynchronizerOperation & op,
std::enable_if_t<std::is_arithmetic<T>::value> * = nullptr) const {
this->allReduceImpl(&values, 1, op);
}
/* ------------------------------------------------------------------------ */
template <typename T> inline void allGather(Array<T> & values) const {
AKANTU_DEBUG_ASSERT(UInt(psize) == values.size(),
"The array size is not correct");
this->allGatherImpl(values.storage(), values.getNbComponent());
}
template <typename Tensor,
typename = std::enable_if_t<is_tensor<Tensor>::value>>
inline void allGather(Tensor & values) const {
AKANTU_DEBUG_ASSERT(values.size() / UInt(psize) > 0,
"The vector size is not correct");
this->allGatherImpl(values.storage(), values.size() / UInt(psize));
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void allGatherV(Array<T> & values, const Array<Int> & sizes) const {
this->allGatherVImpl(values.storage(), sizes.storage());
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void reduce(Array<T> & values, const SynchronizerOperation & op,
int root = 0) const {
this->reduceImpl(values.storage(), values.size() * values.getNbComponent(),
op, root);
}
/* ------------------------------------------------------------------------ */
template <typename Tensor>
inline void
gather(Tensor & values, int root = 0,
std::enable_if_t<is_tensor<Tensor>::value> * = nullptr) const {
this->gatherImpl(values.storage(), values.getNbComponent(), root);
}
template <typename T>
inline void
gather(T values, int root = 0,
std::enable_if_t<std::is_arithmetic<T>::value> * = nullptr) const {
this->gatherImpl(&values, 1, root);
}
/* ------------------------------------------------------------------------ */
template <typename Tensor, typename T>
inline void
gather(Tensor & values, Array<T> & gathered,
std::enable_if_t<is_tensor<Tensor>::value> * = nullptr) const {
AKANTU_DEBUG_ASSERT(values.size() == gathered.getNbComponent(),
"The array size is not correct");
gathered.resize(psize);
this->gatherImpl(values.data(), values.size(), gathered.storage(),
gathered.getNbComponent());
}
template <typename T>
inline void
gather(T values, Array<T> & gathered,
std::enable_if_t<std::is_arithmetic<T>::value> * = nullptr) const {
this->gatherImpl(&values, 1, gathered.storage(), 1);
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void gatherV(Array<T> & values, const Array<Int> & sizes,
int root = 0) const {
this->gatherVImpl(values.storage(), sizes.storage(), root);
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void broadcast(Array<T> & values, int root = 0) const {
this->broadcastImpl(values.storage(),
values.size() * values.getNbComponent(), root);
}
template <bool is_static>
inline void broadcast(CommunicationBufferTemplated<is_static> & values,
int root = 0) const {
this->broadcastImpl(values.storage(), values.size(), root);
}
template <typename T> inline void broadcast(T & values, int root = 0) const {
this->broadcastImpl(&values, 1, root);
}
/* ------------------------------------------------------------------------ */
void barrier() const;
CommunicationRequest asyncBarrier() const;
/* ------------------------------------------------------------------------ */
/* Request handling */
/* ------------------------------------------------------------------------ */
bool test(CommunicationRequest & request) const;
bool testAll(std::vector<CommunicationRequest> & request) const;
void wait(CommunicationRequest & request) const;
void waitAll(std::vector<CommunicationRequest> & requests) const;
UInt waitAny(std::vector<CommunicationRequest> & requests) const;
inline void freeCommunicationRequest(CommunicationRequest & request) const;
inline void
freeCommunicationRequest(std::vector<CommunicationRequest> & requests) const;
template <typename T, typename MsgProcessor>
inline void
receiveAnyNumber(std::vector<CommunicationRequest> & send_requests,
MsgProcessor && processor, Int tag) const;
protected:
template <typename T>
void
sendImpl(const T * buffer, Int size, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const;
template <typename T>
void receiveImpl(T * buffer, Int size, Int sender, Int tag) const;
template <typename T>
CommunicationRequest asyncSendImpl(
const T * buffer, Int size, Int receiver, Int tag,
const CommunicationMode & mode = CommunicationMode::_auto) const;
template <typename T>
CommunicationRequest asyncReceiveImpl(T * buffer, Int size, Int sender,
Int tag) const;
template <typename T>
void allReduceImpl(T * values, int nb_values,
const SynchronizerOperation & op) const;
template <typename T> void allGatherImpl(T * values, int nb_values) const;
template <typename T> void allGatherVImpl(T * values, int * nb_values) const;
template <typename T>
void reduceImpl(T * values, int nb_values, const SynchronizerOperation & op,
int root = 0) const;
template <typename T>
void gatherImpl(T * values, int nb_values, int root = 0) const;
template <typename T>
void gatherImpl(T * values, int nb_values, T * gathered,
int nb_gathered = 0) const;
template <typename T>
void gatherVImpl(T * values, int * nb_values, int root = 0) const;
template <typename T>
void broadcastImpl(T * values, int nb_values, int root = 0) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
Int getNbProc() const { return psize; };
Int whoAmI() const { return prank; };
static Communicator & getStaticCommunicator();
static Communicator & getStaticCommunicator(int & argc, char **& argv);
int getMaxTag() const;
int getMinTag() const;
AKANTU_GET_MACRO(CommunicatorData, (*communicator_data), decltype(auto));
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
static std::unique_ptr<Communicator> static_communicator;
protected:
Int prank{0};
Int psize{1};
std::unique_ptr<CommunicatorInternalData> communicator_data;
};
inline std::ostream & operator<<(std::ostream & stream,
const CommunicationRequest & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "communicator_inline_impl.hh"
#endif /* __AKANTU_STATIC_COMMUNICATOR_HH__ */
diff --git a/src/synchronizer/communicator_dummy_inline_impl.cc b/src/synchronizer/communicator_dummy_inline_impl.cc
index 79aead7f4..3e69923b3 100644
--- a/src/synchronizer/communicator_dummy_inline_impl.cc
+++ b/src/synchronizer/communicator_dummy_inline_impl.cc
@@ -1,117 +1,120 @@
/**
* @file communicator_dummy_inline_impl.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Nov 07 2017
* @date last modification: Fri Nov 10 2017
*
* @brief Dummy communicator to make everything work im sequential
*
* @section LICENSE
*
* Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
/* -------------------------------------------------------------------------- */
#include <cstring>
#include <type_traits>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace akantu {
-Communicator::Communicator(int & /*argc*/, char **& /*argv*/,
- const private_member & /*unused*/) {}
+struct Communicator::private_member {};
+
+Communicator::Communicator(int & /*argc*/, char **& /*argv*/) {}
+
+Communicator::Communicator(const private_member &) {}
template <typename T>
void Communicator::sendImpl(const T *, Int, Int, Int,
const CommunicationMode &) const {}
template <typename T>
void Communicator::receiveImpl(T *, Int, Int, Int) const {}
template <typename T>
CommunicationRequest
Communicator::asyncSendImpl(const T *, Int, Int, Int,
const CommunicationMode &) const {
return std::shared_ptr<InternalCommunicationRequest>(
new InternalCommunicationRequest(0, 0));
}
template <typename T>
CommunicationRequest Communicator::asyncReceiveImpl(T *, Int, Int, Int) const {
return std::shared_ptr<InternalCommunicationRequest>(
new InternalCommunicationRequest(0, 0));
}
template <typename T>
void Communicator::probe(Int, Int, CommunicationStatus &) const {}
template <typename T>
bool Communicator::asyncProbe(Int, Int, CommunicationStatus &) const {
return true;
}
bool Communicator::test(CommunicationRequest &) const { return true; }
bool Communicator::testAll(std::vector<CommunicationRequest> &) const {
return true;
}
void Communicator::wait(CommunicationRequest &) const {}
void Communicator::waitAll(std::vector<CommunicationRequest> &) const {}
UInt Communicator::waitAny(std::vector<CommunicationRequest> &) const {
return UInt(-1);
}
void Communicator::barrier() const {}
CommunicationRequest Communicator::asyncBarrier() const {
return std::shared_ptr<InternalCommunicationRequest>(
new InternalCommunicationRequest(0, 0));
}
template <typename T>
void Communicator::reduceImpl(T *, int, const SynchronizerOperation &,
int) const {}
template <typename T>
void Communicator::allReduceImpl(T *, int,
const SynchronizerOperation &) const {}
template <typename T> inline void Communicator::allGatherImpl(T *, int) const {}
template <typename T>
inline void Communicator::allGatherVImpl(T *, int *) const {}
template <typename T>
inline void Communicator::gatherImpl(T *, int, int) const {}
template <typename T>
void Communicator::gatherImpl(T * values, int nb_values, T * gathered,
int) const {
static_assert(std::is_trivially_copyable<T>{},
"Cannot send this type of data");
std::memcpy(gathered, values, nb_values);
}
template <typename T>
inline void Communicator::gatherVImpl(T *, int *, int) const {}
template <typename T>
inline void Communicator::broadcastImpl(T *, int, int) const {}
int Communicator::getMaxTag() const { return std::numeric_limits<int>::max(); }
int Communicator::getMinTag() const { return 0; }
} // namespace akantu
diff --git a/src/synchronizer/communicator_mpi_inline_impl.cc b/src/synchronizer/communicator_mpi_inline_impl.cc
index fda048047..3bd3e9bcb 100644
--- a/src/synchronizer/communicator_mpi_inline_impl.cc
+++ b/src/synchronizer/communicator_mpi_inline_impl.cc
@@ -1,459 +1,470 @@
/**
* @file communicator_mpi_inline_impl.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Nov 07 2017
* @date last modification: Mon Dec 18 2017
*
* @brief StaticCommunicatorMPI implementation
*
* @section LICENSE
*
* Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_iterators.hh"
#include "communicator.hh"
#include "mpi_communicator_data.hh"
/* -------------------------------------------------------------------------- */
#include <memory>
#include <type_traits>
#include <unordered_map>
#include <vector>
/* -------------------------------------------------------------------------- */
#include <mpi.h>
/* -------------------------------------------------------------------------- */
#if (defined(__GNUC__) || defined(__GNUG__))
#if AKA_GCC_VERSION < 60000
namespace std {
template <> struct hash<akantu::SynchronizerOperation> {
using argument_type = akantu::SynchronizerOperation;
size_t operator()(const argument_type & e) const noexcept {
auto ue = underlying_type_t<argument_type>(e);
return uh(ue);
}
private:
const hash<underlying_type_t<argument_type>> uh{};
};
} // namespace std
#endif
#endif
namespace akantu {
class CommunicationRequestMPI : public InternalCommunicationRequest {
public:
CommunicationRequestMPI(UInt source, UInt dest)
: InternalCommunicationRequest(source, dest),
request(std::make_unique<MPI_Request>()) {}
MPI_Request & getMPIRequest() { return *request; };
private:
std::unique_ptr<MPI_Request> request;
};
namespace {
template <typename T> inline MPI_Datatype getMPIDatatype();
MPI_Op getMPISynchronizerOperation(const SynchronizerOperation & op) {
std::unordered_map<SynchronizerOperation, MPI_Op> _operations{
{SynchronizerOperation::_sum, MPI_SUM},
{SynchronizerOperation::_min, MPI_MIN},
{SynchronizerOperation::_max, MPI_MAX},
{SynchronizerOperation::_prod, MPI_PROD},
{SynchronizerOperation::_land, MPI_LAND},
{SynchronizerOperation::_band, MPI_BAND},
{SynchronizerOperation::_lor, MPI_LOR},
{SynchronizerOperation::_bor, MPI_BOR},
{SynchronizerOperation::_lxor, MPI_LXOR},
{SynchronizerOperation::_bxor, MPI_BXOR},
{SynchronizerOperation::_min_loc, MPI_MINLOC},
{SynchronizerOperation::_max_loc, MPI_MAXLOC},
{SynchronizerOperation::_null, MPI_OP_NULL}};
return _operations[op];
}
template <typename T> MPI_Datatype inline getMPIDatatype() {
return MPI_DATATYPE_NULL;
}
#define SPECIALIZE_MPI_DATATYPE(type, mpi_type) \
template <> MPI_Datatype inline getMPIDatatype<type>() { return mpi_type; }
#define COMMA ,
SPECIALIZE_MPI_DATATYPE(char, MPI_CHAR)
SPECIALIZE_MPI_DATATYPE(float, MPI_FLOAT)
SPECIALIZE_MPI_DATATYPE(double, MPI_DOUBLE)
SPECIALIZE_MPI_DATATYPE(long double, MPI_LONG_DOUBLE)
SPECIALIZE_MPI_DATATYPE(signed int, MPI_INT)
SPECIALIZE_MPI_DATATYPE(NodeType, getMPIDatatype<Int>())
SPECIALIZE_MPI_DATATYPE(unsigned int, MPI_UNSIGNED)
SPECIALIZE_MPI_DATATYPE(signed long int, MPI_LONG)
SPECIALIZE_MPI_DATATYPE(unsigned long int, MPI_UNSIGNED_LONG)
SPECIALIZE_MPI_DATATYPE(signed long long int, MPI_LONG_LONG)
SPECIALIZE_MPI_DATATYPE(unsigned long long int, MPI_UNSIGNED_LONG_LONG)
SPECIALIZE_MPI_DATATYPE(SCMinMaxLoc<double COMMA int>, MPI_DOUBLE_INT)
SPECIALIZE_MPI_DATATYPE(SCMinMaxLoc<float COMMA int>, MPI_FLOAT_INT)
SPECIALIZE_MPI_DATATYPE(bool, MPI_CXX_BOOL)
inline int getMPISource(int src) {
if (src == _any_source)
return MPI_ANY_SOURCE;
return src;
}
decltype(auto) convertRequests(std::vector<CommunicationRequest> & requests) {
std::vector<MPI_Request> mpi_requests(requests.size());
for (auto && request_pair : zip(requests, mpi_requests)) {
auto && req = std::get<0>(request_pair);
auto && mpi_req = std::get<1>(request_pair);
mpi_req = dynamic_cast<CommunicationRequestMPI &>(req.getInternal())
.getMPIRequest();
}
return mpi_requests;
}
} // namespace
// this is ugly but shorten the code a lot
#define MPIDATA \
(*reinterpret_cast<MPICommunicatorData *>(communicator_data.get()))
/* -------------------------------------------------------------------------- */
/* Implementation */
/* -------------------------------------------------------------------------- */
+struct Communicator::private_member {
+ private_member(MPI_Comm comm) : comm(std::move(comm)) {}
+ MPI_Comm comm;
+};
+
/* -------------------------------------------------------------------------- */
-Communicator::Communicator(int & /*argc*/, char **& /*argv*/,
- const private_member & /*unused*/)
+Communicator::Communicator(int & /*argc*/, char **& /*argv*/)
: communicator_data(std::make_unique<MPICommunicatorData>()) {
prank = MPIDATA.rank();
psize = MPIDATA.size();
}
+/* -------------------------------------------------------------------------- */
+Communicator::Communicator(const private_member & p)
+ : communicator_data(std::make_unique<MPICommunicatorData>(p.comm)) {
+ prank = MPIDATA.rank();
+ psize = MPIDATA.size();
+}
+
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::sendImpl(const T * buffer, Int size, Int receiver, Int tag,
const CommunicationMode & mode) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Datatype type = getMPIDatatype<T>();
switch (mode) {
case CommunicationMode::_auto:
MPI_Send(buffer, size, type, receiver, tag, communicator);
break;
case CommunicationMode::_synchronous:
MPI_Ssend(buffer, size, type, receiver, tag, communicator);
break;
case CommunicationMode::_ready:
MPI_Rsend(buffer, size, type, receiver, tag, communicator);
break;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::receiveImpl(T * buffer, Int size, Int sender,
Int tag) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Status status;
MPI_Datatype type = getMPIDatatype<T>();
MPI_Recv(buffer, size, type, getMPISource(sender), tag, communicator,
&status);
}
/* -------------------------------------------------------------------------- */
template <typename T>
CommunicationRequest
Communicator::asyncSendImpl(const T * buffer, Int size, Int receiver, Int tag,
const CommunicationMode & mode) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
auto * request = new CommunicationRequestMPI(prank, receiver);
MPI_Request & req = request->getMPIRequest();
MPI_Datatype type = getMPIDatatype<T>();
switch (mode) {
case CommunicationMode::_auto:
MPI_Isend(buffer, size, type, receiver, tag, communicator, &req);
break;
case CommunicationMode::_synchronous:
MPI_Issend(buffer, size, type, receiver, tag, communicator, &req);
break;
case CommunicationMode::_ready:
MPI_Irsend(buffer, size, type, receiver, tag, communicator, &req);
break;
}
return std::shared_ptr<InternalCommunicationRequest>(request);
}
/* -------------------------------------------------------------------------- */
template <typename T>
CommunicationRequest Communicator::asyncReceiveImpl(T * buffer, Int size,
Int sender, Int tag) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
auto * request = new CommunicationRequestMPI(sender, prank);
MPI_Datatype type = getMPIDatatype<T>();
MPI_Request & req = request->getMPIRequest();
MPI_Irecv(buffer, size, type, getMPISource(sender), tag, communicator, &req);
return std::shared_ptr<InternalCommunicationRequest>(request);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::probe(Int sender, Int tag,
CommunicationStatus & status) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Status mpi_status;
MPI_Probe(getMPISource(sender), tag, communicator, &mpi_status);
MPI_Datatype type = getMPIDatatype<T>();
int count;
MPI_Get_count(&mpi_status, type, &count);
status.setSource(mpi_status.MPI_SOURCE);
status.setTag(mpi_status.MPI_TAG);
status.setSize(count);
}
/* -------------------------------------------------------------------------- */
template <typename T>
bool Communicator::asyncProbe(Int sender, Int tag,
CommunicationStatus & status) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Status mpi_status;
int test;
MPI_Iprobe(getMPISource(sender), tag, communicator, &test, &mpi_status);
if (not test)
return false;
MPI_Datatype type = getMPIDatatype<T>();
int count;
MPI_Get_count(&mpi_status, type, &count);
status.setSource(mpi_status.MPI_SOURCE);
status.setTag(mpi_status.MPI_TAG);
status.setSize(count);
return true;
}
/* -------------------------------------------------------------------------- */
bool Communicator::test(CommunicationRequest & request) const {
MPI_Status status;
int flag;
auto & req_mpi =
dynamic_cast<CommunicationRequestMPI &>(request.getInternal());
MPI_Request & req = req_mpi.getMPIRequest();
#if !defined(AKANTU_NDEBUG)
int ret =
#endif
MPI_Test(&req, &flag, &status);
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Test.");
return (flag != 0);
}
/* -------------------------------------------------------------------------- */
bool Communicator::testAll(std::vector<CommunicationRequest> & requests) const {
int are_finished;
auto && mpi_requests = convertRequests(requests);
MPI_Testall(mpi_requests.size(), mpi_requests.data(), &are_finished,
MPI_STATUSES_IGNORE);
return are_finished != 0 ? true : false;
}
/* -------------------------------------------------------------------------- */
void Communicator::wait(CommunicationRequest & request) const {
MPI_Status status;
auto & req_mpi =
dynamic_cast<CommunicationRequestMPI &>(request.getInternal());
MPI_Request & req = req_mpi.getMPIRequest();
MPI_Wait(&req, &status);
}
/* -------------------------------------------------------------------------- */
void Communicator::waitAll(std::vector<CommunicationRequest> & requests) const {
auto && mpi_requests = convertRequests(requests);
MPI_Waitall(mpi_requests.size(), mpi_requests.data(), MPI_STATUSES_IGNORE);
}
/* -------------------------------------------------------------------------- */
UInt Communicator::waitAny(std::vector<CommunicationRequest> & requests) const {
auto && mpi_requests = convertRequests(requests);
int pos;
MPI_Waitany(mpi_requests.size(), mpi_requests.data(), &pos,
MPI_STATUSES_IGNORE);
if (pos != MPI_UNDEFINED) {
return pos;
} else {
return UInt(-1);
}
}
/* -------------------------------------------------------------------------- */
void Communicator::barrier() const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Barrier(communicator);
}
/* -------------------------------------------------------------------------- */
CommunicationRequest Communicator::asyncBarrier() const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
auto * request = new CommunicationRequestMPI(0, 0);
MPI_Request & req = request->getMPIRequest();
MPI_Ibarrier(communicator, &req);
return std::shared_ptr<InternalCommunicationRequest>(request);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::reduceImpl(T * values, int nb_values,
const SynchronizerOperation & op,
int root) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Datatype type = getMPIDatatype<T>();
MPI_Reduce(MPI_IN_PLACE, values, nb_values, type,
getMPISynchronizerOperation(op), root, communicator);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::allReduceImpl(T * values, int nb_values,
const SynchronizerOperation & op) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Datatype type = getMPIDatatype<T>();
MPI_Allreduce(MPI_IN_PLACE, values, nb_values, type,
getMPISynchronizerOperation(op), communicator);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::allGatherImpl(T * values, int nb_values) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Datatype type = getMPIDatatype<T>();
MPI_Allgather(MPI_IN_PLACE, nb_values, type, values, nb_values, type,
communicator);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::allGatherVImpl(T * values, int * nb_values) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
std::vector<int> displs(psize);
displs[0] = 0;
for (int i = 1; i < psize; ++i) {
displs[i] = displs[i - 1] + nb_values[i - 1];
}
MPI_Datatype type = getMPIDatatype<T>();
MPI_Allgatherv(MPI_IN_PLACE, *nb_values, type, values, nb_values,
displs.data(), type, communicator);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::gatherImpl(T * values, int nb_values, int root) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
T *send_buf = nullptr, *recv_buf = nullptr;
if (prank == root) {
send_buf = (T *)MPI_IN_PLACE;
recv_buf = values;
} else {
send_buf = values;
}
MPI_Datatype type = getMPIDatatype<T>();
MPI_Gather(send_buf, nb_values, type, recv_buf, nb_values, type, root,
communicator);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::gatherImpl(T * values, int nb_values, T * gathered,
int nb_gathered) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
T * send_buf = values;
T * recv_buf = gathered;
if (nb_gathered == 0)
nb_gathered = nb_values;
MPI_Datatype type = getMPIDatatype<T>();
MPI_Gather(send_buf, nb_values, type, recv_buf, nb_gathered, type,
this->prank, communicator);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::gatherVImpl(T * values, int * nb_values, int root) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
int * displs = nullptr;
if (prank == root) {
displs = new int[psize];
displs[0] = 0;
for (int i = 1; i < psize; ++i) {
displs[i] = displs[i - 1] + nb_values[i - 1];
}
}
T *send_buf = nullptr, *recv_buf = nullptr;
if (prank == root) {
send_buf = (T *)MPI_IN_PLACE;
recv_buf = values;
} else
send_buf = values;
MPI_Datatype type = getMPIDatatype<T>();
MPI_Gatherv(send_buf, *nb_values, type, recv_buf, nb_values, displs, type,
root, communicator);
if (prank == root) {
delete[] displs;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Communicator::broadcastImpl(T * values, int nb_values, int root) const {
MPI_Comm communicator = MPIDATA.getMPICommunicator();
MPI_Datatype type = getMPIDatatype<T>();
MPI_Bcast(values, nb_values, type, root, communicator);
}
/* -------------------------------------------------------------------------- */
int Communicator::getMaxTag() const { return MPIDATA.getMaxTag(); }
int Communicator::getMinTag() const { return 0; }
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/synchronizer/mpi_communicator_data.hh b/src/synchronizer/mpi_communicator_data.hh
index f85d30c8b..1191fc288 100644
--- a/src/synchronizer/mpi_communicator_data.hh
+++ b/src/synchronizer/mpi_communicator_data.hh
@@ -1,134 +1,139 @@
/**
* @file mpi_communicator_data.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Mon Feb 05 2018
*
* @brief Wrapper on MPI types to have a better separation between libraries
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#if defined(__INTEL_COMPILER)
//#pragma warning ( disable : 383 )
#elif defined(__clang__) // test clang to be sure that when we test for gnu it
// is only gnu
#elif (defined(__GNUC__) || defined(__GNUG__))
#if __cplusplus > 199711L
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wliteral-suffix"
#endif
#endif
#include <mpi.h>
#if defined(__INTEL_COMPILER)
//#pragma warning ( disable : 383 )
#elif defined(__clang__) // test clang to be sure that when we test for gnu it
// is only gnu
#elif (defined(__GNUC__) || defined(__GNUG__))
#if __cplusplus > 199711L
#pragma GCC diagnostic pop
#endif
#endif
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
/* -------------------------------------------------------------------------- */
#include <unordered_map>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MPI_TYPE_WRAPPER_HH__
#define __AKANTU_MPI_TYPE_WRAPPER_HH__
namespace akantu {
class MPICommunicatorData : public CommunicatorInternalData {
public:
MPICommunicatorData() {
MPI_Initialized(&is_externaly_initialized);
if (!is_externaly_initialized) {
MPI_Init(nullptr, nullptr); // valid according to the spec
}
MPI_Comm_create_errhandler(MPICommunicatorData::errorHandler,
&error_handler);
MPI_Comm_set_errhandler(MPI_COMM_WORLD, error_handler);
setMPICommunicator(MPI_COMM_WORLD);
}
+ MPICommunicatorData(MPI_Comm comm) { communicator = comm; }
+
~MPICommunicatorData() override {
if (not is_externaly_initialized) {
MPI_Finalize();
}
}
inline void setMPICommunicator(MPI_Comm comm) {
MPI_Comm_set_errhandler(communicator, save_error_handler);
communicator = comm;
MPI_Comm_get_errhandler(comm, &save_error_handler);
MPI_Comm_set_errhandler(comm, error_handler);
}
inline int rank() const {
int prank;
MPI_Comm_rank(communicator, &prank);
return prank;
}
inline int size() const {
int psize;
MPI_Comm_size(communicator, &psize);
return psize;
}
inline MPI_Comm getMPICommunicator() const { return communicator; }
inline int getMaxTag() const {
int flag;
int * value;
MPI_Comm_get_attr(communicator, MPI_TAG_UB, &value, &flag);
+ if (!flag) {
+ MPI_Comm_get_attr(MPI_COMM_WORLD, MPI_TAG_UB, &value, &flag);
+ }
AKANTU_DEBUG_ASSERT(flag, "No attribute MPI_TAG_UB.");
return *value;
}
private:
MPI_Comm communicator{MPI_COMM_WORLD};
MPI_Errhandler save_error_handler{MPI_ERRORS_ARE_FATAL};
static int is_externaly_initialized;
/* ------------------------------------------------------------------------ */
MPI_Errhandler error_handler;
static void errorHandler(MPI_Comm * /*comm*/, int * error_code, ...) {
char error_string[MPI_MAX_ERROR_STRING];
int str_len;
MPI_Error_string(*error_code, error_string, &str_len);
AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(),
"MPI failed with the error code "
<< *error_code << ": \"" << error_string
<< "\"");
}
};
} // namespace akantu
#endif /* __AKANTU_MPI_TYPE_WRAPPER_HH__ */