diff --git a/src/common/aka_enum_macros.hh b/src/common/aka_enum_macros.hh
index 5f38371fd..abfa27b5b 100644
--- a/src/common/aka_enum_macros.hh
+++ b/src/common/aka_enum_macros.hh
@@ -1,116 +1,133 @@
/**
* @file aka_enum_macros.hh
*
* @author Nicolas Richart
*
* @date creation Wed Oct 31 2018
*
* @brief A Documented file.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
+#include
+#include
+/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ENUM_MACROS_HH__
#define __AKANTU_AKA_ENUM_MACROS_HH__
#define AKANTU_PP_ENUM(s, data, i, elem) \
BOOST_PP_TUPLE_REM() \
elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data)))
#if (defined(__GNUC__) || defined(__GNUG__))
#define AKA_GCC_VERSION \
(__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
#if AKA_GCC_VERSION < 60000
#define AKANTU_ENUM_HASH(type_name) \
namespace std { \
template <> struct hash<::akantu::type_name> { \
using argument_type = ::akantu::type_name; \
size_t operator()(const argument_type & e) const noexcept { \
auto ue = underlying_type_t(e); \
return uh(ue); \
} \
\
private: \
const hash> uh{}; \
}; \
}
#else
#define AKANTU_ENUM_HASH(type_name)
#endif // AKA_GCC_VERSION
#endif // GNU
#define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem)
#define AKANTU_PP_TYPE_TO_STR(s, data, elem) \
({BOOST_PP_CAT(data, elem), BOOST_PP_STRINGIZE(elem)})
#define AKANTU_PP_STR_TO_TYPE(s, data, elem) \
({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data, elem)})
#define AKANTU_CLASS_ENUM_DECLARE(type_name, list) \
enum class type_name { \
BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \
};
#define AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, prefix) \
} \
AKANTU_ENUM_HASH(type_name) \
namespace std { \
inline string to_string(const ::akantu::type_name & type) { \
using namespace akantu; \
static unordered_map<::akantu::type_name, string> convert{ \
BOOST_PP_SEQ_FOR_EACH_I( \
AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \
BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, prefix, list))}; \
return convert.at(type); \
} \
} \
namespace akantu { \
inline std::ostream & operator<<(std::ostream & stream, \
const type_name & type) { \
stream << std::to_string(type); \
return stream; \
}
#define AKANTU_ENUM_INPUT_STREAM_(type_name, list, prefix) \
inline std::istream & operator>>(std::istream & stream, type_name & type) { \
std::string str; \
stream >> str; \
static std::unordered_map convert{ \
BOOST_PP_SEQ_FOR_EACH_I( \
AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \
BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, prefix, list))}; \
- type = convert.at(str); \
+ try { \
+ type = convert.at(str); \
+ } catch (std::out_of_range &) { \
+ std::ostringstream values; \
+ std::for_each(convert.begin(), convert.end(), [&values](auto && pair) { \
+ static bool first = true; \
+ if (not first) \
+ values << ", "; \
+ values << "\"" << pair.first << "\""; \
+ first = false; \
+ }); \
+ AKANTU_EXCEPTION("The value " << str << " is not a valid " \
+ << BOOST_PP_STRINGIZE(type_name) \
+ << " valid values are " << values.str()); \
+ } \
return stream; \
}
#define AKANTU_CLASS_ENUM_OUTPUT_STREAM(type_name, list) \
AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, type_name::_)
#define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \
AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, )
#define AKANTU_CLASS_ENUM_INPUT_STREAM(type_name, list) \
AKANTU_ENUM_INPUT_STREAM_(type_name, list, type_name::_)
#define AKANTU_ENUM_INPUT_STREAM(type_name, list) \
AKANTU_ENUM_INPUT_STREAM_(type_name, list, )
#endif /* __AKANTU_AKA_ENUM_MACROS_HH__ */
diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc
index b92532bbf..ed51811d6 100644
--- a/src/mesh_utils/mesh_utils.cc
+++ b/src/mesh_utils/mesh_utils.cc
@@ -1,1810 +1,1817 @@
/**
* @file mesh_utils.cc
*
* @author Guillaume Anciaux
* @author Dana Christen
* @author David Simon Kammer
* @author Nicolas Richart
* @author Leonardo Snozzi
* @author Marco Vocialta
*
* @date creation: Fri Aug 20 2010
* @date last modification: Wed Feb 21 2018
*
* @brief All mesh utils necessary for various tasks
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh_utils.hh"
#include "element_synchronizer.hh"
#include "fe_engine.hh"
#include "mesh_accessor.hh"
#include "mesh_iterators.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
void MeshUtils::buildNode2Elements(const Mesh & mesh,
CSR & node_to_elem,
UInt spatial_dimension) {
AKANTU_DEBUG_IN();
if (spatial_dimension == _all_dimensions)
spatial_dimension = mesh.getSpatialDimension();
/// count number of occurrence of each node
UInt nb_nodes = mesh.getNbNodes();
/// array for the node-element list
node_to_elem.resizeRows(nb_nodes);
node_to_elem.clearRows();
for_each_element(mesh,
[&](auto && element) {
Vector conn = mesh.getConnectivity(element);
for (auto && node : conn) {
++node_to_elem.rowOffset(node);
}
},
_spatial_dimension = spatial_dimension,
_element_kind = _ek_not_defined);
node_to_elem.countToCSR();
node_to_elem.resizeCols();
/// rearrange element to get the node-element list
// Element e;
node_to_elem.beginInsertions();
for_each_element(mesh,
[&](auto && element) {
Vector conn = mesh.getConnectivity(element);
for (auto && node : conn) {
node_to_elem.insertInRow(node, element);
}
},
_spatial_dimension = spatial_dimension,
_element_kind = _ek_not_defined);
node_to_elem.endInsertions();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh,
CSR & node_to_elem,
const ElementType & type,
const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_elements = mesh.getConnectivity(type, ghost_type).size();
UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage();
/// array for the node-element list
node_to_elem.resizeRows(nb_nodes);
node_to_elem.clearRows();
/// count number of occurrence of each node
for (UInt el = 0; el < nb_elements; ++el) {
UInt el_offset = el * nb_nodes_per_element;
for (UInt n = 0; n < nb_nodes_per_element; ++n)
++node_to_elem.rowOffset(conn_val[el_offset + n]);
}
/// convert the occurrence array in a csr one
node_to_elem.countToCSR();
node_to_elem.resizeCols();
node_to_elem.beginInsertions();
/// save the element index in the node-element list
for (UInt el = 0; el < nb_elements; ++el) {
UInt el_offset = el * nb_nodes_per_element;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
node_to_elem.insertInRow(conn_val[el_offset + n], el);
}
}
node_to_elem.endInsertions();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildFacets(Mesh & mesh) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
for (auto ghost_type : ghost_types) {
for (auto & type : mesh.elementTypes(spatial_dimension - 1, ghost_type)) {
mesh.getConnectivity(type, ghost_type).resize(0);
// \todo inform the mesh event handler
}
}
buildFacetsDimension(mesh, mesh, true, spatial_dimension);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets,
UInt to_dimension) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets,
UInt from_dimension, UInt to_dimension) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(
mesh_facets.isMeshFacets(),
"The mesh_facets should be initialized with initMeshFacets");
/// generate facets
buildFacetsDimension(mesh, mesh_facets, false, from_dimension);
/// sort facets and generate sub-facets
for (UInt i = from_dimension - 1; i > to_dimension; --i) {
buildFacetsDimension(mesh_facets, mesh_facets, false, i);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets,
bool boundary_only, UInt dimension) {
AKANTU_DEBUG_IN();
// save the current parent of mesh_facets and set it temporarly to mesh since
// mesh is the one containing the elements for which mesh_facets has the
// sub-elements
// example: if the function is called with mesh = mesh_facets
const Mesh * mesh_facets_parent = nullptr;
try {
mesh_facets_parent = &mesh_facets.getMeshParent();
} catch (...) {
}
mesh_facets.defineMeshParent(mesh);
MeshAccessor mesh_accessor(mesh_facets);
UInt spatial_dimension = mesh.getSpatialDimension();
const Array & mesh_facets_nodes = mesh_facets.getNodes();
const auto mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension);
CSR node_to_elem;
buildNode2Elements(mesh, node_to_elem, dimension);
Array counter;
std::vector connected_elements;
// init the SubelementToElement data to improve performance
for (auto && ghost_type : ghost_types) {
for (auto && type : mesh.elementTypes(dimension, ghost_type)) {
mesh_accessor.getSubelementToElement(type, ghost_type);
auto facet_types = mesh.getAllFacetTypes(type);
for (auto && ft : arange(facet_types.size())) {
auto facet_type = facet_types(ft);
mesh_accessor.getElementToSubelement(facet_type, ghost_type);
mesh_accessor.getConnectivity(facet_type, ghost_type);
}
}
}
const ElementSynchronizer * synchronizer = nullptr;
if (mesh.isDistributed()) {
synchronizer = &(mesh.getElementSynchronizer());
}
Element current_element;
for (auto && ghost_type : ghost_types) {
GhostType facet_ghost_type = ghost_type;
current_element.ghost_type = ghost_type;
for (auto && type : mesh.elementTypes(dimension, ghost_type)) {
auto facet_types = mesh.getAllFacetTypes(type);
current_element.type = type;
for (auto && ft : arange(facet_types.size())) {
auto facet_type = facet_types(ft);
auto nb_element = mesh.getNbElement(type, ghost_type);
auto element_to_subelement =
&mesh_facets.getElementToSubelement(facet_type, ghost_type);
auto connectivity_facets =
&mesh_facets.getConnectivity(facet_type, ghost_type);
auto nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft);
const auto & element_connectivity =
mesh.getConnectivity(type, ghost_type);
Matrix facet_local_connectivity(
mesh.getFacetLocalConnectivity(type, ft));
auto nb_nodes_per_facet = connectivity_facets->getNbComponent();
Vector facet(nb_nodes_per_facet);
for (UInt el = 0; el < nb_element; ++el) {
current_element.element = el;
for (UInt f = 0; f < nb_facet_per_element; ++f) {
for (UInt n = 0; n < nb_nodes_per_facet; ++n)
facet(n) =
element_connectivity(el, facet_local_connectivity(f, n));
UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0));
counter.resize(first_node_nb_elements);
counter.clear();
// loop over the other nodes to search intersecting elements,
// which are the elements that share another node with the
// starting element after first_node
UInt local_el = 0;
auto first_node_elements = node_to_elem.begin(facet(0));
auto first_node_elements_end = node_to_elem.end(facet(0));
for (; first_node_elements != first_node_elements_end;
++first_node_elements, ++local_el) {
for (UInt n = 1; n < nb_nodes_per_facet; ++n) {
auto node_elements_begin = node_to_elem.begin(facet(n));
auto node_elements_end = node_to_elem.end(facet(n));
counter(local_el) +=
std::count(node_elements_begin, node_elements_end,
*first_node_elements);
}
}
// counting the number of elements connected to the facets and
// taking the minimum element number, because the facet should
// be inserted just once
UInt nb_element_connected_to_facet = 0;
Element minimum_el = ElementNull;
connected_elements.clear();
for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) {
Element real_el = node_to_elem(facet(0), el_f);
if (not(counter(el_f) == nb_nodes_per_facet - 1))
continue;
++nb_element_connected_to_facet;
minimum_el = std::min(minimum_el, real_el);
connected_elements.push_back(real_el);
}
if (minimum_el != current_element)
continue;
bool full_ghost_facet = false;
UInt n = 0;
while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n)))
++n;
if (n == nb_nodes_per_facet)
full_ghost_facet = true;
if (full_ghost_facet)
continue;
if (boundary_only and nb_element_connected_to_facet != 1)
continue;
std::vector elements;
// build elements_on_facets: linearized_el must come first
// in order to store the facet in the correct direction
// and avoid to invert the sign in the normal computation
elements.push_back(current_element);
if (nb_element_connected_to_facet == 1) { /// boundary facet
elements.push_back(ElementNull);
} else if (nb_element_connected_to_facet == 2) { /// internal facet
elements.push_back(connected_elements[1]);
/// check if facet is in between ghost and normal
/// elements: if it's the case, the facet is either
/// ghost or not ghost. The criterion to decide this
/// is arbitrary. It was chosen to check the processor
/// id (prank) of the two neighboring elements. If
/// prank of the ghost element is lower than prank of
/// the normal one, the facet is not ghost, otherwise
/// it's ghost
GhostType gt[2] = {_not_ghost, _not_ghost};
for (UInt el = 0; el < connected_elements.size(); ++el)
gt[el] = connected_elements[el].ghost_type;
if ((gt[0] == _not_ghost) xor (gt[1] == _not_ghost)) {
UInt prank[2];
for (UInt el = 0; el < 2; ++el) {
prank[el] = synchronizer->getRank(connected_elements[el]);
}
// ugly trick from Marco detected :P
bool ghost_one = (gt[0] != _ghost);
if (prank[ghost_one] > prank[!ghost_one])
facet_ghost_type = _not_ghost;
else
facet_ghost_type = _ghost;
connectivity_facets =
&mesh_facets.getConnectivity(facet_type, facet_ghost_type);
element_to_subelement = &mesh_facets.getElementToSubelement(
facet_type, facet_ghost_type);
}
} else { /// facet of facet
for (UInt i = 1; i < nb_element_connected_to_facet; ++i) {
elements.push_back(connected_elements[i]);
}
}
element_to_subelement->push_back(elements);
connectivity_facets->push_back(facet);
/// current facet index
UInt current_facet = connectivity_facets->size() - 1;
/// loop on every element connected to current facet and
/// insert current facet in the first free spot of the
/// subelement_to_element vector
for (UInt elem = 0; elem < elements.size(); ++elem) {
Element loc_el = elements[elem];
if (loc_el.type == _not_defined)
continue;
Array & subelement_to_element =
mesh_facets.getSubelementToElement(loc_el.type,
loc_el.ghost_type);
UInt nb_facet_per_loc_element =
subelement_to_element.getNbComponent();
for (UInt f_in = 0; f_in < nb_facet_per_loc_element; ++f_in) {
auto & el = subelement_to_element(loc_el.element, f_in);
if (el.type != _not_defined)
continue;
el.type = facet_type;
el.element = current_facet;
el.ghost_type = facet_ghost_type;
break;
}
}
/// reset connectivity in case a facet was found in
/// between ghost and normal elements
if (facet_ghost_type != ghost_type) {
facet_ghost_type = ghost_type;
connectivity_facets =
&mesh_accessor.getConnectivity(facet_type, facet_ghost_type);
element_to_subelement = &mesh_accessor.getElementToSubelement(
facet_type, facet_ghost_type);
}
}
}
}
}
}
// restore the parent of mesh_facet
if (mesh_facets_parent)
mesh_facets.defineMeshParent(*mesh_facets_parent);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::renumberMeshNodes(Mesh & mesh,
Array & local_connectivities,
UInt nb_local_element, UInt nb_ghost_element,
ElementType type,
Array & old_nodes_numbers) {
AKANTU_DEBUG_IN();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
std::map renumbering_map;
for (UInt i = 0; i < old_nodes_numbers.size(); ++i) {
renumbering_map[old_nodes_numbers(i)] = i;
}
/// renumber the nodes
renumberNodesInConnectivity(local_connectivities,
(nb_local_element + nb_ghost_element) *
nb_nodes_per_element,
renumbering_map);
old_nodes_numbers.resize(renumbering_map.size());
for (auto & renumber_pair : renumbering_map) {
old_nodes_numbers(renumber_pair.second) = renumber_pair.first;
}
renumbering_map.clear();
MeshAccessor mesh_accessor(mesh);
/// copy the renumbered connectivity to the right place
auto & local_conn = mesh_accessor.getConnectivity(type);
local_conn.resize(nb_local_element);
- memcpy(local_conn.storage(), local_connectivities.storage(),
- nb_local_element * nb_nodes_per_element * sizeof(UInt));
+
+ if(nb_local_element > 0) {
+ memcpy(local_conn.storage(), local_connectivities.storage(),
+ nb_local_element * nb_nodes_per_element * sizeof(UInt));
+ }
auto & ghost_conn = mesh_accessor.getConnectivity(type, _ghost);
ghost_conn.resize(nb_ghost_element);
- std::memcpy(ghost_conn.storage(),
- local_connectivities.storage() +
- nb_local_element * nb_nodes_per_element,
- nb_ghost_element * nb_nodes_per_element * sizeof(UInt));
-
+
+ if(nb_ghost_element > 0) {
+ std::memcpy(ghost_conn.storage(),
+ local_connectivities.storage() +
+ nb_local_element * nb_nodes_per_element,
+ nb_ghost_element * nb_nodes_per_element * sizeof(UInt));
+ }
+
auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost);
ghost_counter.resize(nb_ghost_element, 1);
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::renumberNodesInConnectivity(
Array & list_nodes, UInt nb_nodes,
std::map & renumbering_map) {
AKANTU_DEBUG_IN();
UInt * connectivity = list_nodes.storage();
UInt new_node_num = renumbering_map.size();
for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) {
UInt & node = *connectivity;
auto it = renumbering_map.find(node);
if (it == renumbering_map.end()) {
UInt old_node = node;
renumbering_map[old_node] = new_node_num;
node = new_node_num;
++new_node_num;
} else {
node = it->second;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::purifyMesh(Mesh & mesh) {
AKANTU_DEBUG_IN();
std::map renumbering_map;
RemovedNodesEvent remove_nodes(mesh);
Array & nodes_removed = remove_nodes.getList();
for (auto ghost_type : ghost_types) {
for (auto type :
mesh.elementTypes(_all_dimensions, ghost_type, _ek_not_defined)) {
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
Array & connectivity = mesh.getConnectivity(type, ghost_type);
UInt nb_element(connectivity.size());
renumberNodesInConnectivity(
connectivity, nb_element * nb_nodes_per_element, renumbering_map);
}
}
Array & new_numbering = remove_nodes.getNewNumbering();
std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1));
for (auto && pair : renumbering_map) {
new_numbering(std::get<0>(pair)) = std::get<1>(pair);
}
for (UInt i = 0; i < new_numbering.size(); ++i) {
if (new_numbering(i) == UInt(-1))
nodes_removed.push_back(i);
}
mesh.sendEvent(remove_nodes);
AKANTU_DEBUG_OUT();
}
#if defined(AKANTU_COHESIVE_ELEMENT)
/* -------------------------------------------------------------------------- */
UInt MeshUtils::insertCohesiveElements(
Mesh & mesh, Mesh & mesh_facets,
const ElementTypeMapArray & facet_insertion,
Array & doubled_nodes, Array & new_elements,
bool only_double_facets) {
UInt spatial_dimension = mesh.getSpatialDimension();
UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion);
if (elements_to_insert > 0) {
if (spatial_dimension == 1) {
doublePointFacet(mesh, mesh_facets, doubled_nodes);
} else {
doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes,
true);
findSubfacetToDouble(mesh_facets);
if (spatial_dimension == 2) {
doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes);
} else if (spatial_dimension == 3) {
doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false);
findSubfacetToDouble(mesh_facets);
doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes);
}
}
if (!only_double_facets)
updateCohesiveData(mesh, mesh_facets, new_elements);
}
return elements_to_insert;
}
#endif
/* -------------------------------------------------------------------------- */
void MeshUtils::doubleNodes(Mesh & mesh, const std::vector & old_nodes,
Array & doubled_nodes) {
AKANTU_DEBUG_IN();
Array & position = mesh.getNodes();
UInt spatial_dimension = mesh.getSpatialDimension();
UInt old_nb_nodes = position.size();
UInt new_nb_nodes = old_nb_nodes + old_nodes.size();
UInt old_nb_doubled_nodes = doubled_nodes.size();
UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size();
position.resize(new_nb_nodes);
doubled_nodes.resize(new_nb_doubled_nodes);
Array::iterator> position_begin =
position.begin(spatial_dimension);
for (UInt n = 0; n < old_nodes.size(); ++n) {
UInt new_node = old_nb_nodes + n;
/// store doubled nodes
doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n];
doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node;
/// update position
std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1,
position_begin + new_node);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets,
UInt facet_dimension, Array & doubled_nodes,
bool facet_mode) {
AKANTU_DEBUG_IN();
for (auto gt_facet : ghost_types) {
for (auto && type_facet :
mesh_facets.elementTypes(facet_dimension, gt_facet)) {
auto & facets_to_double =
mesh_facets.getData("facet_to_double", type_facet, gt_facet);
auto nb_facet_to_double = facets_to_double.size();
if (nb_facet_to_double == 0)
continue;
// this while fail if multiple facet types
// \TODO handle multiple sub-facet types
auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet);
auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet);
auto nb_nodes_per_facet = conn_facet.getNbComponent();
auto old_nb_facet = conn_facet.size();
auto new_nb_facet = old_nb_facet + nb_facet_to_double;
#ifndef AKANTU_NDEBUG
// memory initialization are slow but help debug
conn_facet.resize(new_nb_facet, UInt(-1));
#else
conn_facet.resize(new_nb_facet);
#endif
auto conn_facet_begin = conn_facet.begin(nb_nodes_per_facet);
auto & subfacet_to_facet =
mesh_facets.getSubelementToElement(type_facet, gt_facet);
#ifndef AKANTU_NDEBUG
subfacet_to_facet.resize(new_nb_facet, ElementNull);
#else
subfacet_to_facet.resize(new_nb_facet);
#endif
auto subfacet_to_facet_begin =
subfacet_to_facet.begin(nb_subfacet_per_facet);
Element new_facet{type_facet, old_nb_facet, gt_facet};
auto conn_facet_new_it = conn_facet_begin + new_facet.element;
auto subfacet_to_facet_new_it =
subfacet_to_facet_begin + new_facet.element;
for (UInt facet = 0; facet < nb_facet_to_double; ++facet,
++new_facet.element, ++conn_facet_new_it,
++subfacet_to_facet_new_it) {
UInt old_facet = facets_to_double(facet);
/// adding a new facet by copying original one
/// copy connectivity in new facet
*conn_facet_new_it = conn_facet_begin[old_facet];
/// update subfacet_to_facet
*subfacet_to_facet_new_it = subfacet_to_facet_begin[old_facet];
/// loop on every subfacet
for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) {
Element & subfacet = subfacet_to_facet(old_facet, sf);
if (subfacet == ElementNull)
continue;
/// update facet_to_subfacet array
mesh_facets.getElementToSubelement(subfacet).push_back(new_facet);
}
}
/// update facet_to_subfacet and _segment_3 facets if any
if (not facet_mode) {
updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true);
updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true);
updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet,
doubled_nodes);
} else
updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet,
doubled_nodes);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt MeshUtils::updateFacetToDouble(
Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh_facets.getSpatialDimension();
UInt nb_facets_to_double = 0.;
for (auto gt_facet : ghost_types) {
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
const auto & f_insertion = facet_insertion(type_facet, gt_facet);
auto & f_to_double =
mesh_facets.getData("facet_to_double", type_facet, gt_facet);
auto & element_to_facet =
mesh_facets.getElementToSubelement(type_facet, gt_facet);
Element old_facet_el{type_facet, 0, gt_facet};
UInt nb_facets = mesh_facets.getNbElement(type_facet, gt_facet);
for (UInt f = 0; f < f_insertion.size(); ++f) {
if (f_insertion(f) == false)
continue;
++nb_facets_to_double;
if (element_to_facet(f)[1].type == _not_defined
#if defined(AKANTU_COHESIVE_ELEMENT)
|| element_to_facet(f)[1].kind() == _ek_cohesive
#endif
) {
AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary");
continue;
}
f_to_double.push_back(f);
UInt new_facet = nb_facets + f_to_double.size() - 1;
old_facet_el.element = f;
/// update facet_to_element vector
auto & elem_to_update = element_to_facet(f)[1];
UInt el = elem_to_update.element;
auto & facet_to_element = mesh_facets.getSubelementToElement(
elem_to_update.type, elem_to_update.ghost_type);
auto el_facets = Vector(
make_view(facet_to_element, facet_to_element.getNbComponent())
.begin()[el]);
auto f_update =
std::find(el_facets.begin(), el_facets.end(), old_facet_el);
AKANTU_DEBUG_ASSERT(f_update != el_facets.end(), "Facet not found");
f_update->element = new_facet;
/// update elements connected to facet
const auto & first_facet_list = element_to_facet(f);
element_to_facet.push_back(first_facet_list);
/// set new and original facets as boundary facets
element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1];
element_to_facet(new_facet)[1] = ElementNull;
element_to_facet(f)[1] = ElementNull;
}
}
}
AKANTU_DEBUG_OUT();
return nb_facets_to_double;
}
/* -------------------------------------------------------------------------- */
void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) {
AKANTU_DEBUG_IN();
for (auto gt : ghost_types) {
for (auto type : mesh_facets.elementTypes(_all_dimensions, gt)) {
mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false);
mesh_facets.getDataPointer>(
"facets_to_subfacet_double", type, gt, 1, false);
mesh_facets.getDataPointer>(
"elements_to_subfacet_double", type, gt, 1, false);
mesh_facets.getDataPointer>(
"subfacets_to_subsubfacet_double", type, gt, 1, false);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MeshUtils::findSubfacetToDouble(Mesh & mesh_facets) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh_facets.getSpatialDimension();
if (spatial_dimension == 1)
return;
for (auto gt_facet : ghost_types) {
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
auto & facets_to_double =
mesh_facets.getData("facet_to_double", type_facet, gt_facet);
auto nb_facet_to_double = facets_to_double.size();
if (nb_facet_to_double == 0)
continue;
ElementType type_subfacet = Mesh::getFacetType(type_facet);
GhostType gt_subfacet = _casper;
ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet);
GhostType gt_subsubfacet = _casper;
Array * conn_subfacet = nullptr;
Array * sf_to_double = nullptr;
Array> * sf_to_subfacet_double = nullptr;
Array> * f_to_subfacet_double = nullptr;
Array> * el_to_subfacet_double = nullptr;
UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet);
UInt nb_subsubfacet;
UInt nb_nodes_per_sf_el;
if (subsubfacet_mode) {
nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet);
nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet);
} else
nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet);
Array & subfacet_to_facet =
mesh_facets.getSubelementToElement(type_facet, gt_facet);
Array> & element_to_facet =
mesh_facets.getElementToSubelement(type_facet, gt_facet);
Array * subsubfacet_to_subfacet = nullptr;
UInt old_nb_facet = subfacet_to_facet.size() - nb_facet_to_double;
Element current_facet{type_facet, 0, gt_facet};
std::vector element_list;
std::vector facet_list;
std::vector * subfacet_list;
if (subsubfacet_mode)
subfacet_list = new std::vector;
/// map to filter subfacets
Array> * facet_to_subfacet = nullptr;
/// this is used to make sure that both new and old facets are
/// checked
UInt facets[2];
/// loop on every facet
for (UInt f_index = 0; f_index < 2; ++f_index) {
for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
facets[bool(f_index)] = facets_to_double(facet);
facets[!bool(f_index)] = old_nb_facet + facet;
UInt old_facet = facets[0];
UInt new_facet = facets[1];
Element & starting_element = element_to_facet(new_facet)[0];
current_facet.element = old_facet;
/// loop on every subfacet
for (UInt sf = 0; sf < nb_subfacet; ++sf) {
Element & subfacet = subfacet_to_facet(old_facet, sf);
if (subfacet == ElementNull)
continue;
if (gt_subfacet != subfacet.ghost_type) {
gt_subfacet = subfacet.ghost_type;
if (subsubfacet_mode) {
subsubfacet_to_subfacet = &mesh_facets.getSubelementToElement(
type_subfacet, gt_subfacet);
} else {
conn_subfacet =
&mesh_facets.getConnectivity(type_subfacet, gt_subfacet);
sf_to_double = &mesh_facets.getData(
"facet_to_double", type_subfacet, gt_subfacet);
f_to_subfacet_double =
&mesh_facets.getData>(
"facets_to_subfacet_double", type_subfacet,
gt_subfacet);
el_to_subfacet_double =
&mesh_facets.getData>(
"elements_to_subfacet_double", type_subfacet,
gt_subfacet);
facet_to_subfacet = &mesh_facets.getElementToSubelement(
type_subfacet, gt_subfacet);
}
}
if (subsubfacet_mode) {
/// loop on every subsubfacet
for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) {
Element & subsubfacet =
(*subsubfacet_to_subfacet)(subfacet.element, ssf);
if (subsubfacet == ElementNull)
continue;
if (gt_subsubfacet != subsubfacet.ghost_type) {
gt_subsubfacet = subsubfacet.ghost_type;
conn_subfacet = &mesh_facets.getConnectivity(type_subsubfacet,
gt_subsubfacet);
sf_to_double = &mesh_facets.getData(
"facet_to_double", type_subsubfacet, gt_subsubfacet);
sf_to_subfacet_double =
&mesh_facets.getData>(
"subfacets_to_subsubfacet_double", type_subsubfacet,
gt_subsubfacet);
f_to_subfacet_double =
&mesh_facets.getData>(
"facets_to_subfacet_double", type_subsubfacet,
gt_subsubfacet);
el_to_subfacet_double =
&mesh_facets.getData>(
"elements_to_subfacet_double", type_subsubfacet,
gt_subsubfacet);
facet_to_subfacet = &mesh_facets.getElementToSubelement(
type_subsubfacet, gt_subsubfacet);
}
UInt global_ssf = subsubfacet.element;
Vector subsubfacet_connectivity(
conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el,
nb_nodes_per_sf_el);
/// check if subsubfacet is to be doubled
if (findElementsAroundSubfacet(
mesh_facets, starting_element, current_facet,
subsubfacet_connectivity, element_list, facet_list,
subfacet_list) == false &&
removeElementsInVector(*subfacet_list,
(*facet_to_subfacet)(global_ssf)) ==
false) {
sf_to_double->push_back(global_ssf);
sf_to_subfacet_double->push_back(*subfacet_list);
f_to_subfacet_double->push_back(facet_list);
el_to_subfacet_double->push_back(element_list);
}
}
} else {
const UInt global_sf = subfacet.element;
Vector subfacet_connectivity(
conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el,
nb_nodes_per_sf_el);
/// check if subfacet is to be doubled
if (findElementsAroundSubfacet(
mesh_facets, starting_element, current_facet,
subfacet_connectivity, element_list,
facet_list) == false &&
removeElementsInVector(
facet_list, (*facet_to_subfacet)(global_sf)) == false) {
sf_to_double->push_back(global_sf);
f_to_subfacet_double->push_back(facet_list);
el_to_subfacet_double->push_back(element_list);
}
}
}
}
}
if (subsubfacet_mode)
delete subfacet_list;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
#if defined(AKANTU_COHESIVE_ELEMENT)
void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets,
Array & new_elements) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
bool third_dimension = spatial_dimension == 3;
MeshAccessor mesh_facets_accessor(mesh_facets);
for (auto gt_facet : ghost_types) {
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
Array & f_to_double =
mesh_facets.getData("facet_to_double", type_facet, gt_facet);
UInt nb_facet_to_double = f_to_double.size();
if (nb_facet_to_double == 0)
continue;
ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet);
auto & facet_to_coh_element =
mesh_facets_accessor.getSubelementToElement(type_cohesive, gt_facet);
auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet);
auto & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet);
UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet);
Array> & element_to_facet =
mesh_facets.getElementToSubelement(type_facet, gt_facet);
UInt old_nb_cohesive_elements = conn_cohesive.size();
UInt new_nb_cohesive_elements = conn_cohesive.size() + nb_facet_to_double;
UInt old_nb_facet = element_to_facet.size() - nb_facet_to_double;
facet_to_coh_element.resize(new_nb_cohesive_elements);
conn_cohesive.resize(new_nb_cohesive_elements);
UInt new_elements_old_size = new_elements.size();
new_elements.resize(new_elements_old_size + nb_facet_to_double);
Element c_element{type_cohesive, 0, gt_facet};
Element f_element{type_facet, 0, gt_facet};
UInt facets[2];
for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
/// (in 3D cohesive elements connectivity is inverted)
facets[third_dimension ? 1 : 0] = f_to_double(facet);
facets[third_dimension ? 0 : 1] = old_nb_facet + facet;
UInt cohesive_element = old_nb_cohesive_elements + facet;
/// store doubled facets
f_element.element = facets[0];
facet_to_coh_element(cohesive_element, 0) = f_element;
f_element.element = facets[1];
facet_to_coh_element(cohesive_element, 1) = f_element;
/// modify cohesive elements' connectivity
for (UInt n = 0; n < nb_nodes_per_facet; ++n) {
conn_cohesive(cohesive_element, n) = conn_facet(facets[0], n);
conn_cohesive(cohesive_element, n + nb_nodes_per_facet) =
conn_facet(facets[1], n);
}
/// update element_to_facet vectors
c_element.element = cohesive_element;
element_to_facet(facets[0])[1] = c_element;
element_to_facet(facets[1])[1] = c_element;
/// add cohesive element to the element event list
new_elements(new_elements_old_size + facet) = c_element;
}
}
}
AKANTU_DEBUG_OUT();
}
#endif
/* -------------------------------------------------------------------------- */
void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets,
Array & doubled_nodes) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
if (spatial_dimension != 1)
return;
auto & position = mesh.getNodes();
for (auto gt_facet : ghost_types) {
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet);
auto & element_to_facet =
mesh_facets.getElementToSubelement(type_facet, gt_facet);
const auto & facets_to_double =
mesh_facets.getData("facet_to_double", type_facet, gt_facet);
auto nb_facet_to_double = facets_to_double.size();
auto new_nb_facet = element_to_facet.size();
auto old_nb_facet = element_to_facet.size() - nb_facet_to_double;
auto old_nb_nodes = position.size();
auto new_nb_nodes = old_nb_nodes + nb_facet_to_double;
position.resize(new_nb_nodes);
conn_facet.resize(new_nb_facet);
auto old_nb_doubled_nodes = doubled_nodes.size();
doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double);
for (auto && data_facet : enumerate(facets_to_double)) {
const auto & old_facet = std::get<1>(data_facet);
auto facet = std::get<0>(data_facet);
auto new_facet = old_nb_facet + facet;
auto el = element_to_facet(new_facet)[0];
auto old_node = conn_facet(old_facet);
auto new_node = old_nb_nodes + facet;
/// update position
position(new_node) = position(old_node);
conn_facet(new_facet) = new_node;
Vector conn_segment = mesh.getConnectivity(el);
/// update facet connectivity
auto it = std::find(conn_segment.begin(), conn_segment.end(), old_node);
*it = new_node;
doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node;
doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node;
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets,
ElementType type_facet,
GhostType gt_facet,
Array & doubled_nodes) {
AKANTU_DEBUG_IN();
if (type_facet != _segment_3)
return;
Array & f_to_double =
mesh_facets.getData("facet_to_double", type_facet, gt_facet);
UInt nb_facet_to_double = f_to_double.size();
UInt old_nb_facet =
mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double;
Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet);
Array> & element_to_facet =
mesh_facets.getElementToSubelement(type_facet, gt_facet);
/// this ones matter only for segments in 3D
Array> * el_to_subfacet_double = nullptr;
Array> * f_to_subfacet_double = nullptr;
if (third_dim_segments) {
el_to_subfacet_double = &mesh_facets.getData>(
"elements_to_subfacet_double", type_facet, gt_facet);
f_to_subfacet_double = &mesh_facets.getData>(
"facets_to_subfacet_double", type_facet, gt_facet);
}
std::vector middle_nodes;
for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
UInt old_facet = f_to_double(facet);
UInt node = conn_facet(old_facet, 2);
if (!mesh.isPureGhostNode(node))
middle_nodes.push_back(node);
}
UInt n = doubled_nodes.size();
doubleNodes(mesh, middle_nodes, doubled_nodes);
for (UInt facet = 0; facet < nb_facet_to_double; ++facet) {
UInt old_facet = f_to_double(facet);
UInt old_node = conn_facet(old_facet, 2);
if (mesh.isPureGhostNode(old_node))
continue;
UInt new_node = doubled_nodes(n, 1);
UInt new_facet = old_nb_facet + facet;
conn_facet(new_facet, 2) = new_node;
if (third_dim_segments) {
updateElementalConnectivity(mesh_facets, old_node, new_node,
element_to_facet(new_facet));
updateElementalConnectivity(mesh, old_node, new_node,
(*el_to_subfacet_double)(facet),
&(*f_to_subfacet_double)(facet));
} else {
updateElementalConnectivity(mesh, old_node, new_node,
element_to_facet(new_facet));
}
++n;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets,
ElementType type_subfacet,
GhostType gt_subfacet, bool facet_mode) {
AKANTU_DEBUG_IN();
Array & sf_to_double =
mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet);
UInt nb_subfacet_to_double = sf_to_double.size();
/// update subfacet_to_facet vector
ElementType type_facet = _not_defined;
GhostType gt_facet = _casper;
Array * subfacet_to_facet = nullptr;
UInt nb_subfacet_per_facet = 0;
UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) -
nb_subfacet_to_double;
Array> * facet_list = nullptr;
if (facet_mode)
facet_list = &mesh_facets.getData>(
"facets_to_subfacet_double", type_subfacet, gt_subfacet);
else
facet_list = &mesh_facets.getData>(
"subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet);
Element old_subfacet_el{type_subfacet, 0, gt_subfacet};
Element new_subfacet_el{type_subfacet, 0, gt_subfacet};
for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
old_subfacet_el.element = sf_to_double(sf);
new_subfacet_el.element = old_nb_subfacet + sf;
for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) {
Element & facet = (*facet_list)(sf)[f];
if (facet.type != type_facet || facet.ghost_type != gt_facet) {
type_facet = facet.type;
gt_facet = facet.ghost_type;
subfacet_to_facet =
&mesh_facets.getSubelementToElement(type_facet, gt_facet);
nb_subfacet_per_facet = subfacet_to_facet->getNbComponent();
}
Element * sf_update = std::find(
subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet,
subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet +
nb_subfacet_per_facet,
old_subfacet_el);
AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() +
facet.element * nb_subfacet_per_facet !=
subfacet_to_facet->storage() +
facet.element * nb_subfacet_per_facet +
nb_subfacet_per_facet,
"Subfacet not found");
*sf_update = new_subfacet_el;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets,
ElementType type_subfacet,
GhostType gt_subfacet, bool facet_mode) {
AKANTU_DEBUG_IN();
Array & sf_to_double =
mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet);
UInt nb_subfacet_to_double = sf_to_double.size();
Array> & facet_to_subfacet =
mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet);
Array> * facet_to_subfacet_double = nullptr;
if (facet_mode) {
facet_to_subfacet_double = &mesh_facets.getData>(
"facets_to_subfacet_double", type_subfacet, gt_subfacet);
} else {
facet_to_subfacet_double = &mesh_facets.getData>(
"subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet);
}
UInt old_nb_subfacet = facet_to_subfacet.size();
facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double);
for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf)
facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets,
Array & doubled_nodes) {
AKANTU_DEBUG_IN();
if (spatial_dimension == 1)
return;
for (auto gt_subfacet : ghost_types) {
for (auto type_subfacet : mesh_facets.elementTypes(0, gt_subfacet)) {
auto & sf_to_double = mesh_facets.getData(
"facet_to_double", type_subfacet, gt_subfacet);
UInt nb_subfacet_to_double = sf_to_double.size();
if (nb_subfacet_to_double == 0)
continue;
AKANTU_DEBUG_ASSERT(
type_subfacet == _point_1,
"Only _point_1 subfacet doubling is supported at the moment");
auto & conn_subfacet =
mesh_facets.getConnectivity(type_subfacet, gt_subfacet);
UInt old_nb_subfacet = conn_subfacet.size();
UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double;
conn_subfacet.resize(new_nb_subfacet);
std::vector nodes_to_double;
UInt old_nb_doubled_nodes = doubled_nodes.size();
/// double nodes
for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
UInt old_subfacet = sf_to_double(sf);
nodes_to_double.push_back(conn_subfacet(old_subfacet));
}
doubleNodes(mesh, nodes_to_double, doubled_nodes);
/// add new nodes in connectivity
for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
UInt new_subfacet = old_nb_subfacet + sf;
UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1);
conn_subfacet(new_subfacet) = new_node;
}
/// update facet and element connectivity
Array> & f_to_subfacet_double =
mesh_facets.getData>("facets_to_subfacet_double",
type_subfacet, gt_subfacet);
Array> & el_to_subfacet_double =
mesh_facets.getData>(
"elements_to_subfacet_double", type_subfacet, gt_subfacet);
Array> * sf_to_subfacet_double = nullptr;
if (spatial_dimension == 3)
sf_to_subfacet_double = &mesh_facets.getData>(
"subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet);
for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) {
UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0);
UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1);
updateElementalConnectivity(mesh, old_node, new_node,
el_to_subfacet_double(sf),
&f_to_subfacet_double(sf));
updateElementalConnectivity(mesh_facets, old_node, new_node,
f_to_subfacet_double(sf));
if (spatial_dimension == 3)
updateElementalConnectivity(mesh_facets, old_node, new_node,
(*sf_to_subfacet_double)(sf));
}
if (spatial_dimension == 2) {
updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true);
updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true);
} else if (spatial_dimension == 3) {
updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false);
updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::flipFacets(
Mesh & mesh_facets,
const ElementTypeMapArray & remote_global_connectivities,
GhostType gt_facet) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh_facets.getSpatialDimension();
/// get global connectivity for local mesh
ElementTypeMapArray local_global_connectivities(
"local_global_connectivity", mesh_facets.getID(),
mesh_facets.getMemoryID());
local_global_connectivities.initialize(
mesh_facets, _spatial_dimension = spatial_dimension - 1,
_ghost_type = gt_facet, _with_nb_nodes_per_element = true,
_with_nb_element = true);
mesh_facets.getGlobalConnectivity(local_global_connectivities);
/// loop on every facet
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
auto & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet);
auto & local_global_connectivity =
local_global_connectivities(type_facet, gt_facet);
const auto & remote_global_connectivity =
remote_global_connectivities(type_facet, gt_facet);
auto & element_per_facet =
mesh_facets.getElementToSubelement(type_facet, gt_facet);
auto & subfacet_to_facet =
mesh_facets.getSubelementToElement(type_facet, gt_facet);
auto nb_nodes_per_facet = connectivity.getNbComponent();
auto nb_nodes_per_P1_facet =
Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet));
for (auto && data :
zip(make_view(connectivity, nb_nodes_per_facet),
make_view(local_global_connectivity, nb_nodes_per_facet),
make_view(remote_global_connectivity, nb_nodes_per_facet),
make_view(subfacet_to_facet, subfacet_to_facet.getNbComponent()),
make_view(element_per_facet))) {
auto & conn = std::get<0>(data);
auto & local_gconn = std::get<1>(data);
const auto & remote_gconn = std::get<2>(data);
/// skip facet if connectivities are the same
if (local_gconn == remote_gconn)
continue;
/// re-arrange connectivity
auto conn_tmp = conn;
auto begin = local_gconn.begin();
auto end = local_gconn.end();
std::transform(remote_gconn.begin(), remote_gconn.end(), conn.begin(),
[&](auto && gnode) {
auto it = std::find(begin, end, gnode);
AKANTU_DEBUG_ASSERT(it != end, "Node not found");
return conn_tmp(it - begin);
});
/// if 3D, check if facets are just rotated
if (spatial_dimension == 3) {
auto begin = remote_gconn.storage();
/// find first node
auto it = std::find(begin, begin + remote_gconn.size(), local_gconn(0));
UInt n, start = it - begin;
/// count how many nodes in the received connectivity follow
/// the same order of those in the local connectivity
for (n = 1; n < nb_nodes_per_P1_facet &&
local_gconn(n) ==
remote_gconn((start + n) % nb_nodes_per_P1_facet);
++n)
;
/// skip the facet inversion if facet is just rotated
if (n == nb_nodes_per_P1_facet)
continue;
}
/// update data to invert facet
auto & element_per_facet = std::get<4>(data);
std::swap(element_per_facet[0], element_per_facet[1]);
auto & subfacets_of_facet = std::get<3>(data);
std::swap(subfacets_of_facet(0), subfacets_of_facet(1));
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::fillElementToSubElementsData(Mesh & mesh) {
AKANTU_DEBUG_IN();
if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) {
AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call "
"the buildFacet method.");
return;
}
UInt spatial_dimension = mesh.getSpatialDimension();
ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID(),
mesh.getMemoryID());
barycenters.initialize(mesh, _nb_component = spatial_dimension,
_spatial_dimension = _all_dimensions);
// mesh.initElementTypeMapArray(barycenters, spatial_dimension,
// _all_dimensions);
Element element;
for (auto ghost_type : ghost_types) {
element.ghost_type = ghost_type;
for (auto & type : mesh.elementTypes(_all_dimensions, ghost_type)) {
element.type = type;
UInt nb_element = mesh.getNbElement(type, ghost_type);
Array & barycenters_arr = barycenters(type, ghost_type);
barycenters_arr.resize(nb_element);
auto bary = barycenters_arr.begin(spatial_dimension);
auto bary_end = barycenters_arr.end(spatial_dimension);
for (UInt el = 0; bary != bary_end; ++bary, ++el) {
element.element = el;
mesh.getBarycenter(element, *bary);
}
}
}
MeshAccessor mesh_accessor(mesh);
for (Int sp(spatial_dimension); sp >= 1; --sp) {
if (mesh.getNbElement(sp) == 0)
continue;
for (auto ghost_type : ghost_types) {
for (auto & type : mesh.elementTypes(sp, ghost_type)) {
mesh_accessor.getSubelementToElement(type, ghost_type)
.resize(mesh.getNbElement(type, ghost_type));
mesh_accessor.getSubelementToElement(type, ghost_type).set(ElementNull);
}
for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) {
mesh_accessor.getElementToSubelement(type, ghost_type)
.resize(mesh.getNbElement(type, ghost_type));
mesh.getElementToSubelement(type, ghost_type).clear();
}
}
CSR nodes_to_elements;
buildNode2Elements(mesh, nodes_to_elements, sp);
Element facet_element;
for (auto ghost_type : ghost_types) {
facet_element.ghost_type = ghost_type;
for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) {
facet_element.type = type;
auto & element_to_subelement =
mesh.getElementToSubelement(type, ghost_type);
const auto & connectivity = mesh.getConnectivity(type, ghost_type);
for (auto && data : enumerate(
make_view(connectivity, mesh.getNbNodesPerElement(type)))) {
const auto & facet = std::get<1>(data);
facet_element.element = std::get<0>(data);
std::map element_seen_counter;
auto nb_nodes_per_facet =
mesh.getNbNodesPerElement(Mesh::getP1ElementType(type));
// count the number of node in common between the facet and the other
// element connected to the nodes of the facet
for (auto node : arange(nb_nodes_per_facet)) {
for (auto & elem : nodes_to_elements.getRow(facet(node))) {
auto cit = element_seen_counter.find(elem);
if (cit != element_seen_counter.end()) {
cit->second++;
} else {
element_seen_counter[elem] = 1;
}
}
}
// check which are the connected elements
std::vector connected_elements;
for (auto && cit : element_seen_counter) {
if (cit.second == nb_nodes_per_facet)
connected_elements.push_back(cit.first);
}
// add the connected elements as sub-elements
for (auto & connected_element : connected_elements) {
element_to_subelement(facet_element.element)
.push_back(connected_element);
}
// add the element as sub-element to the connected elements
for (auto & connected_element : connected_elements) {
Vector subelements_to_element =
mesh.getSubelementToElement(connected_element);
// find the position where to insert the element
auto it = std::find(subelements_to_element.begin(),
subelements_to_element.end(), ElementNull);
AKANTU_DEBUG_ASSERT(
it != subelements_to_element.end(),
"The element "
<< connected_element << " seems to have too many facets!! ("
<< (it - subelements_to_element.begin()) << " < "
<< mesh.getNbFacetsPerElement(connected_element.type)
<< ")");
*it = facet_element;
}
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
bool MeshUtils::findElementsAroundSubfacet(
const Mesh & mesh_facets, const Element & starting_element,
const Element & end_facet, const Vector & subfacet_connectivity,
std::vector & element_list, std::vector & facet_list,
std::vector * subfacet_list) {
AKANTU_DEBUG_IN();
bool facet_matched = false;
element_list.clear();
facet_list.clear();
if (third_dim_points) {
subfacet_list->clear();
}
element_list.push_back(starting_element);
std::queue elements_to_check;
elements_to_check.push(starting_element);
/// keep going as long as there are elements to check
while (not elements_to_check.empty()) {
/// check current element
Element & current_element = elements_to_check.front();
const Vector facets_to_element =
mesh_facets.getSubelementToElement(current_element);
// for every facet of the element
for (auto & current_facet : facets_to_element) {
if (current_facet == ElementNull)
continue;
if (current_facet == end_facet)
facet_matched = true;
// facet already listed
if (std::find(facet_list.begin(), facet_list.end(), current_facet) !=
facet_list.end())
continue;
// subfacet_connectivity is not in the connectivity of current_facet;
if ((std::find(facet_list.begin(), facet_list.end(), current_facet) !=
facet_list.end()) or
not hasElement(mesh_facets.getConnectivity(current_facet),
subfacet_connectivity))
continue;
facet_list.push_back(current_facet);
if (third_dim_points) {
const Vector subfacets_of_facet =
mesh_facets.getSubelementToElement(current_facet);
/// check subfacets
for (const auto & current_subfacet : subfacets_of_facet) {
if (current_subfacet == ElementNull)
continue;
if ((std::find(subfacet_list->begin(), subfacet_list->end(),
current_subfacet) == subfacet_list->end()) and
hasElement(mesh_facets.getConnectivity(current_subfacet),
subfacet_connectivity))
subfacet_list->push_back(current_subfacet);
}
}
/// consider opposing element
const auto & elements_to_facet =
mesh_facets.getElementToSubelement(current_facet);
UInt opposing = 0;
if (elements_to_facet[0] == current_element)
opposing = 1;
auto & opposing_element = elements_to_facet[opposing];
/// skip null elements since they are on a boundary
if (opposing_element == ElementNull)
continue;
/// skip this element if already added
if (std::find(element_list.begin(), element_list.end(),
opposing_element) != element_list.end())
continue;
/// only regular elements have to be checked
if (opposing_element.kind() == _ek_regular)
elements_to_check.push(opposing_element);
element_list.push_back(opposing_element);
AKANTU_DEBUG_ASSERT(
hasElement(
mesh_facets.getMeshParent().getConnectivity(opposing_element),
subfacet_connectivity),
"Subfacet doesn't belong to this element");
}
/// erased checked element from the list
elements_to_check.pop();
}
AKANTU_DEBUG_OUT();
return facet_matched;
}
/* -------------------------------------------------------------------------- */
void MeshUtils::updateElementalConnectivity(
Mesh & mesh, UInt old_node, UInt new_node,
const std::vector & element_list,
const std::vector *
#if defined(AKANTU_COHESIVE_ELEMENT)
facet_list
#endif
) {
AKANTU_DEBUG_IN();
for (auto & element : element_list) {
if (element.type == _not_defined)
continue;
Vector connectivity = mesh.getConnectivity(element);
#if defined(AKANTU_COHESIVE_ELEMENT)
if (element.kind() == _ek_cohesive) {
AKANTU_DEBUG_ASSERT(
facet_list != nullptr,
"Provide a facet list in order to update cohesive elements");
const Vector facets =
mesh.getMeshFacets().getSubelementToElement(element);
auto facet_nb_nodes = connectivity.size() / 2;
/// loop over cohesive element's facets
for (const auto & facet : enumerate(facets)) {
/// skip facets if not present in the list
if (std::find(facet_list->begin(), facet_list->end(),
std::get<1>(facet)) == facet_list->end()) {
continue;
}
auto n = std::get<0>(facet);
auto begin = connectivity.begin() + n * facet_nb_nodes;
auto end = begin + facet_nb_nodes;
auto it = std::find(begin, end, old_node);
AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element");
*it = new_node;
}
} else
#endif
{
auto it = std::find(connectivity.begin(), connectivity.end(), old_node);
AKANTU_DEBUG_ASSERT(it != connectivity.end(),
"Node not found in current element");
/// update connectivity
*it = new_node;
}
}
AKANTU_DEBUG_OUT();
}
} // namespace akantu
diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc
index 9a1289a88..0961702e6 100644
--- a/src/solver/sparse_matrix_aij.cc
+++ b/src/solver/sparse_matrix_aij.cc
@@ -1,291 +1,291 @@
/**
* @file sparse_matrix_aij.cc
*
* @author Nicolas Richart
*
* @date creation: Fri Aug 21 2015
* @date last modification: Mon Dec 04 2017
*
* @brief Implementation of the AIJ sparse matrix
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "sparse_matrix_aij.hh"
#include "aka_iterators.hh"
#include "dof_manager_default.hh"
#include "dof_synchronizer.hh"
#include "solver_vector_default.hh"
#include "terms_to_assemble.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager,
const MatrixType & matrix_type, const ID & id)
: SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager),
irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a") {}
/* -------------------------------------------------------------------------- */
SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id)
: SparseMatrix(matrix, id), dof_manager(matrix.dof_manager),
irn(matrix.irn, id + ":irn"), jcn(matrix.jcn, id + ":jcn"),
a(matrix.a, id + ":a") {}
/* -------------------------------------------------------------------------- */
SparseMatrixAIJ::~SparseMatrixAIJ() = default;
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::applyBoundary(Real block_val) {
AKANTU_DEBUG_IN();
const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs();
auto begin = blocked_dofs.begin();
auto end = blocked_dofs.end();
auto is_blocked = [&](auto && i) -> bool {
auto il = this->dof_manager.globalToLocalEquationNumber(i);
return std::binary_search(begin, end, il);
};
for (auto && ij_a : zip(irn, jcn, a)) {
UInt ni = std::get<0>(ij_a) - 1;
UInt nj = std::get<1>(ij_a) - 1;
if (is_blocked(ni) or is_blocked(nj)) {
// clang-format off
std::get<2>(ij_a) =
std::get<0>(ij_a) != std::get<1>(ij_a) ? 0.
: this->dof_manager.isLocalOrMasterDOF(ni) ? block_val
: 0.;
// clang-format on
}
}
this->value_release++;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::saveProfile(const std::string & filename) const {
AKANTU_DEBUG_IN();
std::ofstream outfile;
outfile.open(filename.c_str());
UInt m = this->size_;
auto & comm = dof_manager.getCommunicator();
// write header
if (comm.whoAmI() == 0) {
outfile << "%%MatrixMarket matrix coordinate pattern";
if (this->matrix_type == _symmetric)
outfile << " symmetric";
else
outfile << " general";
outfile << std::endl;
outfile << m << " " << m << " " << this->nb_non_zero << std::endl;
}
for (auto p : arange(comm.getNbProc())) {
// write content
if (comm.whoAmI() == p) {
for (UInt i = 0; i < this->nb_non_zero; ++i) {
outfile << this->irn.storage()[i] << " " << this->jcn.storage()[i]
<< " 1" << std::endl;
}
}
comm.barrier();
}
outfile.close();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::saveMatrix(const std::string & filename) const {
AKANTU_DEBUG_IN();
auto & comm = dof_manager.getCommunicator();
// open and set the properties of the stream
std::ofstream outfile;
if (0 == comm.whoAmI()) {
outfile.open(filename.c_str());
} else {
outfile.open(filename.c_str(), std::ios_base::app);
}
outfile.precision(std::numeric_limits::digits10);
// write header
decltype(nb_non_zero) nnz = this->nb_non_zero;
comm.allReduce(nnz);
if (comm.whoAmI() == 0) {
outfile << "%%MatrixMarket matrix coordinate real";
if (this->matrix_type == _symmetric)
outfile << " symmetric";
else
outfile << " general";
outfile << std::endl;
outfile << this->size_ << " " << this->size_ << " " << nnz << std::endl;
}
for (auto p : arange(comm.getNbProc())) {
// write content
if (comm.whoAmI() == p) {
for (UInt i = 0; i < this->nb_non_zero; ++i) {
outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i)
<< std::endl;
}
}
comm.barrier();
}
// time to end
outfile.close();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::matVecMul(const Array & x, Array & y,
Real alpha, Real beta) const {
AKANTU_DEBUG_IN();
y *= beta;
auto i_it = this->irn.begin();
auto j_it = this->jcn.begin();
auto a_it = this->a.begin();
auto a_end = this->a.end();
auto x_it = x.begin_reinterpret(x.size() * x.getNbComponent());
auto y_it = y.begin_reinterpret(x.size() * x.getNbComponent());
for (; a_it != a_end; ++i_it, ++j_it, ++a_it) {
Int i = this->dof_manager.globalToLocalEquationNumber(*i_it - 1);
Int j = this->dof_manager.globalToLocalEquationNumber(*j_it - 1);
const Real & A = *a_it;
y_it[i] += alpha * A * x_it[j];
if ((this->matrix_type == _symmetric) && (i != j))
y_it[j] += alpha * A * x_it[i];
}
if (this->dof_manager.hasSynchronizer())
this->dof_manager.getSynchronizer().reduceSynchronize(y);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::matVecMul(const SolverVector & _x, SolverVector & _y,
Real alpha, Real beta) const {
AKANTU_DEBUG_IN();
auto && x = aka::as_type(_x).getVector();
auto && y = aka::as_type(_y).getVector();
this->matVecMul(x, y, alpha, beta);
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) {
AKANTU_DEBUG_IN();
const auto & mat = aka::as_type(matrix);
AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(),
"The to matrix don't have the same profiles");
memcpy(a.storage(), mat.getA().storage(), nb_non_zero * sizeof(Real));
this->value_release++;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::copyProfile(const SparseMatrix & other) {
auto & A = aka::as_type(other);
SparseMatrix::clearProfile();
- this->irn = A.irn;
- this->jcn = A.jcn;
+ this->irn.copy(A.irn);
+ this->jcn.copy(A.jcn);
this->irn_jcn_k.clear();
UInt i, j, k;
for(auto && data : enumerate(irn, jcn)) {
std::tie(k, i, j) = data;
this->irn_jcn_k[this->key(i - 1, j - 1)] = k;
}
this->nb_non_zero = this->irn.size();
this->a.resize(this->nb_non_zero);
this->a.set(0.);
this->size_ = A.size_;
- this->profile_release++;
+ this->profile_release = A.profile_release;
this->value_release++;
}
/* -------------------------------------------------------------------------- */
template
void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const {
UInt i, j;
Real A_ij;
for (auto && tuple : zip(irn, jcn, a)) {
std::tie(i, j, A_ij) = tuple;
B.add(i - 1, j - 1, alpha * A_ij);
}
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const {
if (aka::is_of_type(B)) {
this->addMeToTemplated(aka::as_type(B), alpha);
} else {
// this->addMeToTemplated(*this, alpha);
}
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::mul(Real alpha) {
this->a *= alpha;
this->value_release++;
}
/* -------------------------------------------------------------------------- */
void SparseMatrixAIJ::clear() {
a.set(0.);
this->value_release++;
}
} // namespace akantu
diff --git a/src/synchronizer/node_info_per_processor.cc b/src/synchronizer/node_info_per_processor.cc
index ddae37c59..f81bcae63 100644
--- a/src/synchronizer/node_info_per_processor.cc
+++ b/src/synchronizer/node_info_per_processor.cc
@@ -1,851 +1,849 @@
/**
* @file node_info_per_processor.cc
*
* @author Nicolas Richart
*
* @date creation: Wed Mar 16 2016
* @date last modification: Wed Nov 08 2017
*
* @brief Please type the brief for file: Helper classes to create the
* distributed synchronizer and distribute a mesh
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "node_info_per_processor.hh"
#include "communicator.hh"
#include "node_group.hh"
#include "node_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
NodeInfoPerProc::NodeInfoPerProc(NodeSynchronizer & synchronizer,
UInt message_cnt, UInt root)
: MeshAccessor(synchronizer.getMesh()), synchronizer(synchronizer),
comm(synchronizer.getCommunicator()), rank(comm.whoAmI()),
nb_proc(comm.getNbProc()), root(root), mesh(synchronizer.getMesh()),
spatial_dimension(synchronizer.getMesh().getSpatialDimension()),
message_count(message_cnt) {}
/* -------------------------------------------------------------------------- */
void NodeInfoPerProc::synchronize() {
synchronizeNodes();
synchronizeTypes();
synchronizeGroups();
synchronizePeriodicity();
synchronizeTags();
}
/* -------------------------------------------------------------------------- */
template
void NodeInfoPerProc::fillNodeGroupsFromBuffer(CommunicationBuffer & buffer) {
AKANTU_DEBUG_IN();
std::vector> node_to_group;
buffer >> node_to_group;
AKANTU_DEBUG_ASSERT(node_to_group.size() == mesh.getNbGlobalNodes(),
"Not the good amount of nodes where transmitted");
const auto & global_nodes = mesh.getGlobalNodesIds();
for (auto && data : enumerate(global_nodes)) {
for (const auto & node : node_to_group[std::get<1>(data)]) {
mesh.getNodeGroup(node).add(std::get<0>(data), false);
}
}
for (auto && ng_data : mesh.iterateNodeGroups()) {
ng_data.optimize();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void NodeInfoPerProc::fillNodesType() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
auto & nodes_flags = this->getNodesFlags();
Array nodes_set(nb_nodes);
nodes_set.set(0);
enum NodeSet {
NORMAL_SET = 1,
GHOST_SET = 2,
};
Array already_seen(nb_nodes, 1, false);
for (auto gt : ghost_types) {
UInt set = NORMAL_SET;
if (gt == _ghost)
set = GHOST_SET;
already_seen.set(false);
for (auto && type :
mesh.elementTypes(_all_dimensions, gt, _ek_not_defined)) {
const auto & connectivity = mesh.getConnectivity(type, gt);
for (auto & conn :
make_view(connectivity, connectivity.getNbComponent())) {
for (UInt n = 0; n < conn.size(); ++n) {
AKANTU_DEBUG_ASSERT(conn(n) < nb_nodes,
"Node " << conn(n)
<< " bigger than number of nodes "
<< nb_nodes);
if (!already_seen(conn(n))) {
nodes_set(conn(n)) += set;
already_seen(conn(n)) = true;
}
}
}
}
}
nodes_flags.resize(nb_nodes);
for (UInt i = 0; i < nb_nodes; ++i) {
if (nodes_set(i) == NORMAL_SET)
nodes_flags(i) = NodeFlag::_normal;
else if (nodes_set(i) == GHOST_SET)
nodes_flags(i) = NodeFlag::_pure_ghost;
else if (nodes_set(i) == (GHOST_SET + NORMAL_SET))
nodes_flags(i) = NodeFlag::_master;
else {
AKANTU_EXCEPTION("Gni ?");
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void NodeInfoPerProc::fillCommunicationScheme(const Array & master_info) {
AKANTU_DEBUG_IN();
Communications & communications =
this->synchronizer.getCommunications();
{ // send schemes
- auto it = master_info.begin_reinterpret(2, master_info.size() / 2);
- auto end = master_info.end_reinterpret(2, master_info.size() / 2);
-
std::map> send_array_per_proc;
- for (; it != end; ++it) {
- const Vector & send_info = *it;
-
+ for (const auto & send_info : make_view(master_info, 2)) {
send_array_per_proc[send_info(0)].push_back(send_info(1));
}
for (auto & send_schemes : send_array_per_proc) {
auto & scheme = communications.createSendScheme(send_schemes.first);
auto & sends = send_schemes.second;
std::sort(sends.begin(), sends.end());
std::transform(sends.begin(), sends.end(), sends.begin(),
[this](UInt g) -> UInt { return mesh.getNodeLocalId(g); });
scheme.copy(sends);
+ std::cout << "Proc " << rank << " sends " << sends.size()
+ << " nodes to proc " << send_schemes.first << std::endl;
}
}
{ // receive schemes
std::map> recv_array_per_proc;
for (auto node : arange(mesh.getNbNodes())) {
if (mesh.isSlaveNode(node)) {
recv_array_per_proc[mesh.getNodePrank(node)].push_back(
mesh.getNodeGlobalId(node));
}
}
for (auto & recv_schemes : recv_array_per_proc) {
auto & scheme = communications.createRecvScheme(recv_schemes.first);
auto & recvs = recv_schemes.second;
std::sort(recvs.begin(), recvs.end());
std::transform(recvs.begin(), recvs.end(), recvs.begin(),
[this](UInt g) -> UInt { return mesh.getNodeLocalId(g); });
scheme.copy(recvs);
+ std::cout << "Proc " << rank << " receives " << recvs.size()
+ << " nodes to proc " << recv_schemes.first << std::endl;
+
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void NodeInfoPerProc::fillPeriodicPairs(const Array & global_pairs,
std::vector & missing_nodes) {
this->wipePeriodicInfo();
auto & nodes_flags = this->getNodesFlags();
auto checkIsLocal = [&](auto && global_node) {
auto && node = mesh.getNodeLocalId(global_node);
if (node == UInt(-1)) {
auto & global_nodes = this->getNodesGlobalIds();
node = global_nodes.size();
global_nodes.push_back(global_node);
nodes_flags.push_back(NodeFlag::_pure_ghost);
missing_nodes.push_back(global_node);
std::cout << "Missing node " << node << std::endl;
}
return node;
};
for (auto && pairs : make_view(global_pairs, 2)) {
UInt slave = checkIsLocal(pairs(0));
UInt master = checkIsLocal(pairs(1));
this->addPeriodicSlave(slave, master);
}
this->markMeshPeriodic();
}
/* -------------------------------------------------------------------------- */
void NodeInfoPerProc::receiveMissingPeriodic(
DynamicCommunicationBuffer & buffer) {
auto & nodes = this->getNodes();
Communications & communications =
this->synchronizer.getCommunications();
std::size_t nb_nodes;
buffer >> nb_nodes;
- for (auto _[[gnu::unused]] : arange(nb_nodes)) {
+ for (auto _ [[gnu::unused]] : arange(nb_nodes)) {
Vector pos(spatial_dimension);
Int prank;
buffer >> pos;
buffer >> prank;
UInt node = nodes.size();
this->setNodePrank(node, prank);
nodes.push_back(pos);
auto & scheme = communications.createRecvScheme(prank);
scheme.push_back(node);
}
while (buffer.getLeftToUnpack() != 0) {
Int prank;
UInt gnode;
buffer >> gnode;
buffer >> prank;
auto node = mesh.getNodeLocalId(gnode);
AKANTU_DEBUG_ASSERT(node != UInt(-1),
"I cannot send the node "
<< gnode << " to proc " << prank
<< " because it is note a local node");
auto & scheme = communications.createSendScheme(prank);
scheme.push_back(node);
}
}
/* -------------------------------------------------------------------------- */
void NodeInfoPerProc::fillNodalData(DynamicCommunicationBuffer & buffer,
std::string tag_name) {
#define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, _, elem) \
case MeshDataTypeCode::BOOST_PP_TUPLE_ELEM(2, 0, elem): { \
auto & nodal_data = \
mesh.getNodalData(tag_name); \
nodal_data.resize(mesh.getNbNodes()); \
for (auto && data : make_view(nodal_data)) { \
buffer >> data; \
} \
break; \
}
MeshDataTypeCode data_type_code =
mesh.getTypeCode(tag_name, MeshDataType::_nodal);
switch (data_type_code) {
BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, ,
AKANTU_MESH_DATA_TYPES)
default:
AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!");
break;
}
#undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
MasterNodeInfoPerProc::MasterNodeInfoPerProc(NodeSynchronizer & synchronizer,
UInt message_cnt, UInt root)
: NodeInfoPerProc(synchronizer, message_cnt, root),
all_nodes(0, synchronizer.getMesh().getSpatialDimension()) {
UInt nb_global_nodes = this->mesh.getNbGlobalNodes();
this->comm.broadcast(nb_global_nodes, this->root);
}
/* -------------------------------------------------------------------------- */
void MasterNodeInfoPerProc::synchronizeNodes() {
this->nodes_per_proc.resize(nb_proc);
this->nb_nodes_per_proc.resize(nb_proc);
Array local_nodes(0, spatial_dimension);
Array & nodes = this->getNodes();
all_nodes.copy(nodes);
nodes_pranks.resize(nodes.size(), UInt(-1));
for (UInt p = 0; p < nb_proc; ++p) {
UInt nb_nodes = 0;
// UInt * buffer;
Array * nodes_to_send;
Array & nodespp = nodes_per_proc[p];
if (p != root) {
nodes_to_send = new Array(0, spatial_dimension);
AKANTU_DEBUG_INFO("Receiving number of nodes from proc "
<< p << " " << Tag::genTag(p, 0, Tag::_NB_NODES));
comm.receive(nb_nodes, p, Tag::genTag(p, 0, Tag::_NB_NODES));
nodespp.resize(nb_nodes);
this->nb_nodes_per_proc(p) = nb_nodes;
AKANTU_DEBUG_INFO("Receiving list of nodes from proc "
<< p << " " << Tag::genTag(p, 0, Tag::_NODES));
comm.receive(nodespp, p, Tag::genTag(p, 0, Tag::_NODES));
} else {
Array & local_ids = this->getNodesGlobalIds();
this->nb_nodes_per_proc(p) = local_ids.size();
nodespp.copy(local_ids);
nodes_to_send = &local_nodes;
}
- Array::const_scalar_iterator it = nodespp.begin();
- Array::const_scalar_iterator end = nodespp.end();
/// get the coordinates for the selected nodes
- for (; it != end; ++it) {
- Vector coord(nodes.storage() + spatial_dimension * *it,
+ for (const auto & node : nodespp) {
+ Vector coord(nodes.storage() + spatial_dimension * node,
spatial_dimension);
nodes_to_send->push_back(coord);
}
if (p != root) { /// send them for distant processors
AKANTU_DEBUG_INFO("Sending coordinates to proc "
<< p << " "
<< Tag::genTag(this->rank, 0, Tag::_COORDINATES));
comm.send(*nodes_to_send, p,
Tag::genTag(this->rank, 0, Tag::_COORDINATES));
delete nodes_to_send;
}
}
/// construct the local nodes coordinates
nodes.copy(local_nodes);
}
/* -------------------------------------------------------------------------- */
void MasterNodeInfoPerProc::synchronizeTypes() {
// >
std::multimap> nodes_to_proc;
std::vector> nodes_flags_per_proc(nb_proc);
std::vector> nodes_prank_per_proc(nb_proc);
if (mesh.isPeriodic())
all_periodic_flags.copy(this->getNodesFlags());
// arrays containing pairs of (proc, node)
std::vector> nodes_to_send_per_proc(nb_proc);
for (UInt p = 0; p < nb_proc; ++p) {
nodes_flags_per_proc[p].resize(nb_nodes_per_proc(p), NodeFlag(0xFF));
nodes_prank_per_proc[p].resize(nb_nodes_per_proc(p), -1);
}
this->fillNodesType();
auto is_master = [](auto && flag) {
return (flag & NodeFlag::_shared_mask) == NodeFlag::_master;
};
auto is_local = [](auto && flag) {
return (flag & NodeFlag::_shared_mask) == NodeFlag::_normal;
};
for (auto p : arange(nb_proc)) {
auto & nodes_flags = nodes_flags_per_proc[p];
if (p != root) {
AKANTU_DEBUG_INFO(
"Receiving first nodes types from proc "
<< p << " "
<< Tag::genTag(this->rank, this->message_count, Tag::_NODES_TYPE));
comm.receive(nodes_flags, p, Tag::genTag(p, 0, Tag::_NODES_TYPE));
} else {
nodes_flags.copy(this->getNodesFlags());
}
// stack all processors claiming to be master for a node
for (auto local_node : arange(nb_nodes_per_proc(p))) {
auto global_node = nodes_per_proc[p](local_node);
if (is_master(nodes_flags(local_node))) {
nodes_to_proc.insert(
std::make_pair(global_node, std::make_pair(p, local_node)));
} else if (is_local(nodes_flags(local_node))) {
nodes_pranks[global_node] = p;
}
}
}
for (auto i : arange(mesh.getNbGlobalNodes())) {
auto it_range = nodes_to_proc.equal_range(i);
if (it_range.first == nodes_to_proc.end() || it_range.first->first != i)
continue;
// pick the first processor out of the multi-map as the actual master
UInt master_proc = (it_range.first)->second.first;
nodes_pranks[i] = master_proc;
for (auto && data : range(it_range.first, it_range.second)) {
auto proc = data.second.first;
auto node = data.second.second;
if (proc != master_proc) {
// store the info on all the slaves for a given master
nodes_flags_per_proc[proc](node) = NodeFlag::_slave;
nodes_to_send_per_proc[master_proc].push_back(proc);
nodes_to_send_per_proc[master_proc].push_back(i);
}
}
}
/// Fills the nodes prank per proc
for (auto && data : zip(arange(nb_proc), nodes_per_proc, nodes_prank_per_proc,
nodes_flags_per_proc)) {
for (auto && node_data :
zip(std::get<1>(data), std::get<2>(data), std::get<3>(data))) {
if (std::get<2>(node_data) == NodeFlag::_normal) {
std::get<1>(node_data) = std::get<0>(data);
} else {
std::get<1>(node_data) = nodes_pranks(std::get<0>(node_data));
}
}
}
std::vector requests_send_type;
std::vector requests_send_master_info;
for (UInt p = 0; p < nb_proc; ++p) {
if (p != root) {
AKANTU_DEBUG_INFO("Sending nodes types to proc "
<< p << " "
<< Tag::genTag(this->rank, 0, Tag::_NODES_TYPE));
requests_send_type.push_back(
comm.asyncSend(nodes_flags_per_proc[p], p,
Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)));
requests_send_type.push_back(
comm.asyncSend(nodes_prank_per_proc[p], p,
Tag::genTag(this->rank, 2, Tag::_NODES_TYPE)));
auto & nodes_to_send = nodes_to_send_per_proc[p];
AKANTU_DEBUG_INFO("Sending nodes master info to proc "
<< p << " "
<< Tag::genTag(this->rank, 1, Tag::_NODES_TYPE));
requests_send_master_info.push_back(comm.asyncSend(
nodes_to_send, p, Tag::genTag(this->rank, 1, Tag::_NODES_TYPE)));
} else {
this->getNodesFlags().copy(nodes_flags_per_proc[p]);
for (auto && data : enumerate(nodes_prank_per_proc[p])) {
auto node = std::get<0>(data);
if (not(mesh.isMasterNode(node) or mesh.isLocalNode(node))) {
this->setNodePrank(node, std::get<1>(data));
}
}
this->fillCommunicationScheme(nodes_to_send_per_proc[root]);
}
}
comm.waitAll(requests_send_type);
comm.freeCommunicationRequest(requests_send_type);
requests_send_type.clear();
comm.waitAll(requests_send_master_info);
comm.freeCommunicationRequest(requests_send_master_info);
}
/* -------------------------------------------------------------------------- */
void MasterNodeInfoPerProc::synchronizeGroups() {
AKANTU_DEBUG_IN();
UInt nb_total_nodes = mesh.getNbGlobalNodes();
DynamicCommunicationBuffer buffer;
using NodeToGroup = std::vector>;
NodeToGroup node_to_group;
node_to_group.resize(nb_total_nodes);
GroupManager::const_node_group_iterator ngi = mesh.node_group_begin();
GroupManager::const_node_group_iterator nge = mesh.node_group_end();
for (; ngi != nge; ++ngi) {
NodeGroup & ng = *(ngi->second);
std::string name = ngi->first;
NodeGroup::const_node_iterator nit = ng.begin();
NodeGroup::const_node_iterator nend = ng.end();
for (; nit != nend; ++nit) {
node_to_group[*nit].push_back(name);
}
nit = ng.begin();
if (nit != nend)
ng.empty();
}
buffer << node_to_group;
std::vector requests;
for (UInt p = 0; p < nb_proc; ++p) {
if (p == this->rank)
continue;
AKANTU_DEBUG_INFO("Sending node groups to proc "
<< p << " "
<< Tag::genTag(this->rank, p, Tag::_NODE_GROUP));
requests.push_back(comm.asyncSend(
buffer, p, Tag::genTag(this->rank, p, Tag::_NODE_GROUP)));
}
this->fillNodeGroupsFromBuffer(buffer);
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
requests.clear();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MasterNodeInfoPerProc::synchronizePeriodicity() {
bool is_periodic = mesh.isPeriodic();
comm.broadcast(is_periodic, root);
if (not is_periodic)
return;
std::vector requests;
std::vector> periodic_info_to_send_per_proc;
for (auto p : arange(nb_proc)) {
periodic_info_to_send_per_proc.emplace_back(0, 2);
auto && periodic_info = periodic_info_to_send_per_proc.back();
for (UInt proc_local_node : arange(nb_nodes_per_proc(p))) {
UInt global_node = nodes_per_proc[p](proc_local_node);
if ((all_periodic_flags[global_node] & NodeFlag::_periodic_mask) ==
NodeFlag::_periodic_slave) {
periodic_info.push_back(
Vector{global_node, mesh.getPeriodicMaster(global_node)});
}
}
if (p == root)
continue;
auto && tag = Tag::genTag(this->rank, p, Tag::_PERIODIC_SLAVES);
AKANTU_DEBUG_INFO("Sending periodic info to proc " << p << " " << tag);
requests.push_back(comm.asyncSend(periodic_info, p, tag));
}
CommunicationStatus status;
std::vector buffers(nb_proc);
std::vector> proc_missings(nb_proc);
auto nodes_it = all_nodes.begin(spatial_dimension);
for (UInt p = 0; p < nb_proc; ++p) {
auto & proc_missing = proc_missings[p];
if (p != root) {
auto && tag = Tag::genTag(p, 0, Tag::_PERIODIC_NODES);
comm.probe(p, tag, status);
proc_missing.resize(status.size());
comm.receive(proc_missing, p, tag);
} else {
fillPeriodicPairs(periodic_info_to_send_per_proc[root], proc_missing);
}
auto & buffer = buffers[p];
buffer.reserve((spatial_dimension * sizeof(Real) + sizeof(Int)) *
proc_missing.size());
buffer << proc_missing.size();
for (auto && node : proc_missing) {
buffer << *(nodes_it + node);
buffer << nodes_pranks(node);
}
}
for (UInt p = 0; p < nb_proc; ++p) {
for (auto && node : proc_missings[p]) {
auto & buffer = buffers[nodes_pranks(node)];
buffer << node;
buffer << p;
}
}
for (UInt p = 0; p < nb_proc; ++p) {
if (p != root) {
auto && tag_send = Tag::genTag(p, 1, Tag::_PERIODIC_NODES);
requests.push_back(comm.asyncSend(buffers[p], p, tag_send));
} else {
receiveMissingPeriodic(buffers[p]);
}
}
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
requests.clear();
}
/* -------------------------------------------------------------------------- */
void MasterNodeInfoPerProc::fillTagBuffers(
std::vector & buffers,
const std::string & tag_name) {
#define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, _, elem) \
case MeshDataTypeCode::BOOST_PP_TUPLE_ELEM(2, 0, elem): { \
auto & nodal_data = \
mesh.getNodalData(tag_name); \
for (auto && data : enumerate(nodes_per_proc)) { \
auto proc = std::get<0>(data); \
auto & nodes = std::get<1>(data); \
auto & buffer = buffers[proc]; \
for (auto & node : nodes) { \
for (auto i : arange(nodal_data.getNbComponent())) { \
buffer << nodal_data(node, i); \
} \
} \
} \
break; \
}
MeshDataTypeCode data_type_code =
mesh.getTypeCode(tag_name, MeshDataType::_nodal);
switch (data_type_code) {
BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, ,
AKANTU_MESH_DATA_TYPES)
default:
AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!");
break;
}
#undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA
} // namespace akantu
/* -------------------------------------------------------------------------- */
void MasterNodeInfoPerProc::synchronizeTags() {
/// tag info
auto tag_names = mesh.getTagNames();
DynamicCommunicationBuffer tags_buffer;
for (auto && tag_name : tag_names) {
tags_buffer << tag_name;
tags_buffer << mesh.getTypeCode(tag_name, MeshDataType::_nodal);
tags_buffer << mesh.getNbComponent(tag_name);
}
AKANTU_DEBUG_INFO(
"Broadcasting the information about the nodes mesh data tags: ("
<< tags_buffer.size() << ").");
comm.broadcast(tags_buffer, root);
for (auto && tag_data : enumerate(tag_names)) {
auto tag_count = std::get<0>(tag_data);
auto & tag_name = std::get<1>(tag_data);
std::vector buffers;
std::vector requests;
buffers.resize(nb_proc);
fillTagBuffers(buffers, tag_name);
for (auto && data : enumerate(buffers)) {
auto && proc = std::get<0>(data);
auto & buffer = std::get<1>(data);
if (proc == root) {
fillNodalData(buffer, tag_name);
} else {
auto && tag = Tag::genTag(this->rank, tag_count, Tag::_MESH_DATA);
requests.push_back(comm.asyncSend(buffer, proc, tag));
}
}
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
}
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
SlaveNodeInfoPerProc::SlaveNodeInfoPerProc(NodeSynchronizer & synchronizer,
UInt message_cnt, UInt root)
: NodeInfoPerProc(synchronizer, message_cnt, root) {
UInt nb_global_nodes = 0;
comm.broadcast(nb_global_nodes, root);
this->setNbGlobalNodes(nb_global_nodes);
}
/* -------------------------------------------------------------------------- */
void SlaveNodeInfoPerProc::synchronizeNodes() {
AKANTU_DEBUG_INFO("Sending list of nodes to proc "
<< root << " " << Tag::genTag(this->rank, 0, Tag::_NB_NODES)
<< " " << Tag::genTag(this->rank, 0, Tag::_NODES));
Array & local_ids = this->getNodesGlobalIds();
Array & nodes = this->getNodes();
UInt nb_nodes = local_ids.size();
comm.send(nb_nodes, root, Tag::genTag(this->rank, 0, Tag::_NB_NODES));
comm.send(local_ids, root, Tag::genTag(this->rank, 0, Tag::_NODES));
/* --------<<<<-COORDINATES---------------------------------------------- */
nodes.resize(nb_nodes);
AKANTU_DEBUG_INFO("Receiving coordinates from proc "
<< root << " " << Tag::genTag(root, 0, Tag::_COORDINATES));
comm.receive(nodes, root, Tag::genTag(root, 0, Tag::_COORDINATES));
}
/* -------------------------------------------------------------------------- */
void SlaveNodeInfoPerProc::synchronizeTypes() {
this->fillNodesType();
auto & nodes_types = this->getNodesFlags();
AKANTU_DEBUG_INFO("Sending first nodes types to proc "
<< root << ""
<< Tag::genTag(this->rank, 0, Tag::_NODES_TYPE));
comm.send(nodes_types, root, Tag::genTag(this->rank, 0, Tag::_NODES_TYPE));
AKANTU_DEBUG_INFO("Receiving nodes types from proc "
<< root << " " << Tag::genTag(root, 0, Tag::_NODES_TYPE));
comm.receive(nodes_types, root, Tag::genTag(root, 0, Tag::_NODES_TYPE));
Array nodes_prank(nodes_types.size());
comm.receive(nodes_prank, root, Tag::genTag(root, 2, Tag::_NODES_TYPE));
for (auto && data : enumerate(nodes_prank)) {
auto node = std::get<0>(data);
if (not(mesh.isMasterNode(node) or mesh.isLocalNode(node))) {
this->setNodePrank(node, std::get<1>(data));
}
}
AKANTU_DEBUG_INFO("Receiving nodes master info from proc "
<< root << " " << Tag::genTag(root, 1, Tag::_NODES_TYPE));
CommunicationStatus status;
comm.probe(root, Tag::genTag(root, 1, Tag::_NODES_TYPE), status);
Array nodes_master_info(status.size());
if (nodes_master_info.size() > 0)
comm.receive(nodes_master_info, root,
Tag::genTag(root, 1, Tag::_NODES_TYPE));
this->fillCommunicationScheme(nodes_master_info);
}
/* -------------------------------------------------------------------------- */
void SlaveNodeInfoPerProc::synchronizeGroups() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Receiving node groups from proc "
<< root << " "
<< Tag::genTag(root, this->rank, Tag::_NODE_GROUP));
DynamicCommunicationBuffer buffer;
comm.receive(buffer, root, Tag::genTag(root, this->rank, Tag::_NODE_GROUP));
this->fillNodeGroupsFromBuffer(buffer);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SlaveNodeInfoPerProc::synchronizePeriodicity() {
bool is_periodic;
comm.broadcast(is_periodic, root);
if (not is_periodic)
return;
CommunicationStatus status;
auto && tag = Tag::genTag(root, this->rank, Tag::_PERIODIC_SLAVES);
comm.probe(root, tag, status);
Array periodic_info(status.size() / 2, 2);
comm.receive(periodic_info, root, tag);
std::vector proc_missing;
fillPeriodicPairs(periodic_info, proc_missing);
auto && tag_missing_request =
Tag::genTag(this->rank, 0, Tag::_PERIODIC_NODES);
comm.send(proc_missing, root, tag_missing_request);
DynamicCommunicationBuffer buffer;
auto && tag_missing = Tag::genTag(this->rank, 1, Tag::_PERIODIC_NODES);
comm.receive(buffer, root, tag_missing);
receiveMissingPeriodic(buffer);
}
/* -------------------------------------------------------------------------- */
void SlaveNodeInfoPerProc::synchronizeTags() {
DynamicCommunicationBuffer tags_buffer;
comm.broadcast(tags_buffer, root);
std::vector tag_names;
while (tags_buffer.getLeftToUnpack() > 0) {
std::string name;
MeshDataTypeCode code;
UInt nb_components;
tags_buffer >> name;
tags_buffer >> code;
tags_buffer >> nb_components;
mesh.registerNodalData(name, nb_components, code);
tag_names.push_back(name);
}
for (auto && tag_data : enumerate(tag_names)) {
auto tag_count = std::get<0>(tag_data);
auto & tag_name = std::get<1>(tag_data);
DynamicCommunicationBuffer buffer;
auto && tag = Tag::genTag(this->root, tag_count, Tag::_MESH_DATA);
comm.receive(buffer, this->root, tag);
fillNodalData(buffer, tag_name);
}
}
} // namespace akantu
diff --git a/test/test_gtest_utils.hh b/test/test_gtest_utils.hh
index 98551f437..f136abf9b 100644
--- a/test/test_gtest_utils.hh
+++ b/test/test_gtest_utils.hh
@@ -1,232 +1,242 @@
/**
* @file test_gtest_utils.hh
*
* @author Nicolas Richart
*
* @date creation: Tue Nov 14 2017
* @date last modification: Wed Feb 21 2018
*
* @brief Utils to help write tests
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_iterators.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_TEST_GTEST_UTILS_HH__
#define __AKANTU_TEST_GTEST_UTILS_HH__
namespace {
/* -------------------------------------------------------------------------- */
template <::akantu::ElementType t>
using element_type_t = std::integral_constant<::akantu::ElementType, t>;
/* -------------------------------------------------------------------------- */
template struct gtest_list {};
template struct gtest_list> {
using type = ::testing::Types;
};
template using gtest_list_t = typename gtest_list::type;
/* -------------------------------------------------------------------------- */
template struct tuple_concat {};
template
struct tuple_concat, std::tuple> {
using type = std::tuple;
};
template
using tuple_concat_t = typename tuple_concat::type;
/* -------------------------------------------------------------------------- */
template class Pred, typename... Ts>
struct tuple_filter {};
template class Pred, typename T>
struct tuple_filter> {
using type = std::conditional_t::value, std::tuple, std::tuple<>>;
};
template class Pred, typename T, typename... Ts>
struct tuple_filter> {
using type =
tuple_concat_t>::type,
typename tuple_filter>::type>;
};
template class Pred, typename... Ts>
using tuple_filter_t = typename tuple_filter::type;
/* -------------------------------------------------------------------------- */
template struct tuple_split {};
template
struct tuple_split> {
protected:
using split = tuple_split>;
public:
using type = tuple_concat_t, typename split::type>;
using type_tail = typename split::type_tail;
};
template
struct tuple_split<1, std::tuple> {
using type = std::tuple;
using type_tail = std::tuple;
};
template
using tuple_split_t = typename tuple_split::type;
template
using tuple_split_tail_t = typename tuple_split::type_tail;
/* -------------------------------------------------------------------------- */
template struct cross_product {};
template