diff --git a/src/common/aka_iterators.hh b/src/common/aka_iterators.hh
index 33860debb..5329181f4 100644
--- a/src/common/aka_iterators.hh
+++ b/src/common/aka_iterators.hh
@@ -1,309 +1,320 @@
/**
* @file aka_iterators.hh
*
* @author Nicolas Richart
*
* @date creation Wed Jul 19 2017
*
* @brief iterator interfaces
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_error.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
+#include
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ITERATORS_HH__
#define __AKANTU_AKA_ITERATORS_HH__
namespace akantu {
namespace tuple {
/* ------------------------------------------------------------------------ */
namespace details {
template struct Foreach {
template
static inline decltype(auto) transform_forward(F && func,
Tuple && tuple) {
return std::tuple_cat(
Foreach::transform_forward(std::forward(func),
std::forward(tuple)),
std::forward_as_tuple(std::forward(func)(
std::get(std::forward(tuple)))));
}
template
static inline decltype(auto) transform(F && func, Tuple && tuple) {
return std::tuple_cat(
Foreach::transform(std::forward(func),
std::forward(tuple)),
std::make_tuple(std::forward(func)(
std::get(std::forward(tuple)))));
}
template
static inline void foreach (F && func, Tuple && tuple) {
Foreach::foreach (std::forward(func),
std::forward(tuple));
std::forward(func)(std::get(std::forward(tuple)));
}
template static inline bool equal(Tuple && a, Tuple && b) {
if (not(std::get(std::forward(a)) ==
std::get(std::forward(b))))
return false;
return Foreach::equal(std::forward(a),
std::forward(b));
}
};
/* ------------------------------------------------------------------------
*/
template <> struct Foreach<1> {
template
static inline decltype(auto) transform_forward(F && func,
Tuple && tuple) {
return std::forward_as_tuple(
std::forward(func)(std::get<0>(std::forward(tuple))));
}
template
static inline decltype(auto) transform(F && func, Tuple && tuple) {
return std::make_tuple(
std::forward(func)(std::get<0>(std::forward(tuple))));
}
template
static inline void foreach (F && func, Tuple && tuple) {
std::forward(func)(std::get<0>(std::forward(tuple)));
}
template static inline bool equal(Tuple && a, Tuple && b) {
return std::get<0>(std::forward(a)) ==
std::get<0>(std::forward(b));
}
};
} // namespace details
/* ------------------------------------------------------------------------ */
template bool are_equal(Tuple && a, Tuple && b) {
return details::Foreach>::value>::equal(
std::forward(a), std::forward(b));
}
template void foreach (F && func, Tuple && tuple) {
details::Foreach>::value>::foreach (
std::forward(func), std::forward(tuple));
}
template
decltype(auto) transform_forward(F && func, Tuple && tuple) {
return details::Foreach>::value>::
transform_forward(std::forward(func), std::forward(tuple));
}
template
decltype(auto) transform(F && func, Tuple && tuple) {
return details::Foreach>::value>::
transform(std::forward(func), std::forward(tuple));
}
} // namespace tuple
namespace iterators {
namespace details {
struct dereference_iterator {
template
- typename Iter::reference operator()(Iter & it) const {
- return *it;
+ decltype(auto) operator()(Iter & it) const {
+ return std::forward(*it);
}
};
struct increment_iterator {
template void operator()(Iter & it) const { ++it; }
};
struct begin_container {
template
decltype(auto) operator()(Container && cont) const {
return std::forward(cont).begin();
}
};
struct end_container {
template
decltype(auto) operator()(Container && cont) const {
return std::forward(cont).end();
}
};
} // namespace details
/* ------------------------------------------------------------------------ */
template class ZipIterator {
private:
using tuple_t = std::tuple;
-
public:
explicit ZipIterator(tuple_t iterators) : iterators(std::move(iterators)) {}
decltype(auto) operator*() {
return tuple::transform_forward(details::dereference_iterator(),
iterators);
}
ZipIterator & operator++() {
tuple::foreach (details::increment_iterator(), iterators);
return *this;
}
bool operator==(const ZipIterator & other) const {
return tuple::are_equal(iterators, other.iterators);
}
bool operator!=(const ZipIterator & other) const {
return not operator==(other);
}
private:
tuple_t iterators;
};
} // namespace iterators
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template
decltype(auto) zip_iterator(std::tuple && iterators_tuple) {
- auto zip = iterators::ZipIterator(std::move(iterators_tuple));
+ auto zip = iterators::ZipIterator(std::forward(iterators_tuple));
return zip;
}
/* -------------------------------------------------------------------------- */
namespace containers {
template class ZipContainer {
- using containers_t = std::tuple;
+ using containers_t = std::tuple;
public:
explicit ZipContainer(Containers &&... containers)
- : containers(std::forward_as_tuple(containers...)) {}
+ : containers(std::forward(containers)...) {}
decltype(auto) begin() const {
return zip_iterator(
tuple::transform(iterators::details::begin_container(),
std::forward(containers)));
}
decltype(auto) end() const {
return zip_iterator(
tuple::transform(iterators::details::end_container(),
std::forward(containers)));
}
decltype(auto) begin() {
return zip_iterator(
tuple::transform(iterators::details::begin_container(),
std::forward(containers)));
}
decltype(auto) end() {
return zip_iterator(
tuple::transform(iterators::details::end_container(),
std::forward(containers)));
}
private:
containers_t containers;
};
} // namespace containers
/* -------------------------------------------------------------------------- */
template decltype(auto) zip(Containers &&... conts) {
- return containers::ZipContainer(conts...);
+ return containers::ZipContainer(std::forward(conts)...);
}
/* -------------------------------------------------------------------------- */
/* Arange */
/* -------------------------------------------------------------------------- */
namespace iterators {
template class ArangeIterator {
public:
+ using value_type = T;
+ using pointer = T*;
+ using reference = T&;
+ using iterator_category = std::input_iterator_tag;
+
constexpr ArangeIterator(T value, T step) : value(value), step(step) {}
constexpr ArangeIterator(const ArangeIterator &) = default;
constexpr ArangeIterator & operator++() {
value += step;
return *this;
}
- constexpr T operator*() { return value; }
+ constexpr const T & operator*() const { return value; }
- constexpr bool operator==(const ArangeIterator & other) {
- return value == other.value and step == other.step;
+ constexpr bool operator==(const ArangeIterator & other) const {
+ return (value == other.value) and (step == other.step);
}
- constexpr bool operator!=(const ArangeIterator & other) {
+ constexpr bool operator!=(const ArangeIterator & other) const {
return not operator==(other);
}
private:
- T value;
- const T step;
+ T value{0};
+ const T step{1};
};
} // namespace iterators
namespace containers {
template class ArangeContainer {
public:
using iterator = iterators::ArangeIterator;
constexpr ArangeContainer(T start, T stop, T step = 1)
: start(start), stop((stop - start)% step == 0 ? stop : start + (1 + (stop -start)/step)*step), step(step) {}
explicit constexpr ArangeContainer(T stop) : ArangeContainer(0, stop, 1) {}
constexpr T operator[](size_t i) {
T val = start + i * step;
assert(val < stop && "i is out of range");
return val;
}
constexpr size_t size() { return (stop - start) / step; }
constexpr iterator begin() { return iterator(start, step); }
constexpr iterator end() { return iterator(stop, step); }
private:
const T start{0}, stop{0}, step{1};
};
} // namespace containers
template
inline decltype(auto) arange(T stop) {
return containers::ArangeContainer(stop);
}
template
-inline decltype(auto) arange(T start, T stop, T step = 1) {
+inline constexpr decltype(auto) arange(T start, T stop, T step = 1) {
return containers::ArangeContainer(start, stop, step);
}
+template
+inline constexpr decltype(auto) counting(Container && container) {
+ return zip(arange(std::forward(container).size()),
+ std::forward(container));
+}
+
} // namespace akantu
#endif /* __AKANTU_AKA_ITERATORS_HH__ */
diff --git a/src/common/aka_static_if.hh b/src/common/aka_static_if.hh
new file mode 100644
index 000000000..274a05295
--- /dev/null
+++ b/src/common/aka_static_if.hh
@@ -0,0 +1,105 @@
+// Copyright (c) 2016 Vittorio Romeo
+// License: AFL 3.0 | https://opensource.org/licenses/AFL-3.0
+// http://vittorioromeo.info | vittorio.romeo@outlook.com
+
+#ifndef __AKANTU_AKA_STATIC_IF_HH__
+#define __AKANTU_AKA_STATIC_IF_HH__
+
+#include
+
+#define FWD(...) ::std::forward(__VA_ARGS__)
+
+namespace akantu {
+template using bool_ = std::integral_constant;
+
+template constexpr bool_ bool_v{};
+
+template using int_ = std::integral_constant;
+
+template constexpr int_ int_v{};
+
+template using sz_ = std::integral_constant;
+
+template constexpr sz_ sz_v{};
+
+template auto static_if(TPredicate) noexcept;
+
+namespace impl {
+ template struct static_if_impl;
+
+ template struct static_if_result;
+
+ template auto make_static_if_result(TF && f) noexcept;
+
+ template <> struct static_if_impl {
+ template auto & else_(TF &&) noexcept {
+ // Ignore `else_`, as the predicate is true.
+ return *this;
+ }
+
+ template auto & else_if(TPredicate) noexcept {
+ // Ignore `else_if`, as the predicate is true.
+ return *this;
+ }
+
+ template auto then(TF && f) noexcept {
+ // We found a matching branch, just make a result and
+ // ignore everything else.
+ return make_static_if_result(FWD(f));
+ }
+ };
+
+ template <> struct static_if_impl {
+ template auto & then(TF &&) noexcept {
+ // Ignore `then`, as the predicate is false.
+ return *this;
+ }
+
+ template auto else_(TF && f) noexcept {
+ // (Assuming that `else_` is after all `else_if` calls.)
+
+ // We found a matching branch, just make a result and
+ // ignore everything else.
+
+ return make_static_if_result(FWD(f));
+ }
+
+ template auto else_if(TPredicate) noexcept {
+ return static_if(TPredicate{});
+ }
+
+ template auto operator()(Ts &&...) noexcept {
+ // If there are no `else` branches, we must ignore calls
+ // to a failed `static_if` matching.
+ }
+ };
+
+ template
+ struct static_if_result : TFunctionToCall {
+ // Perfect-forward the function in the result instance.
+ template
+ explicit static_if_result(TFFwd && f) noexcept : TFunctionToCall(FWD(f)) {}
+
+ // Ignore everything, we found a result.
+ template auto & then(TF &&) noexcept { return *this; }
+
+ template auto & else_if(TPredicate) noexcept {
+ return *this;
+ }
+
+ template auto & else_(TF &&) noexcept { return *this; }
+ };
+
+ template auto make_static_if_result(TF && f) noexcept {
+ return static_if_result{FWD(f)};
+ }
+} // namespace impl
+
+template auto static_if(TPredicate) noexcept {
+ return impl::static_if_impl{};
+}
+
+#undef FWD
+} // namespace akantu
+
+#endif /* __AKANTU_AKA_STATIC_IF_HH__ */
diff --git a/src/mesh/mesh_iterators.hh b/src/mesh/mesh_iterators.hh
index 1121458ba..7969b01dc 100644
--- a/src/mesh/mesh_iterators.hh
+++ b/src/mesh/mesh_iterators.hh
@@ -1,224 +1,237 @@
/**
* @file mesh_iterators.hh
*
* @author Nicolas Richart
*
* @date creation Wed Aug 09 2017
*
* @brief Set of helper classes to have fun with range based for
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_ITERATORS_HH__
#define __AKANTU_MESH_ITERATORS_HH__
namespace akantu {
class MeshElementsByTypes {
using elements_iterator = Array::scalar_iterator;
public:
explicit MeshElementsByTypes(const Array & elements) {
this->elements.copy(elements);
std::sort(this->elements.begin(), this->elements.end());
}
/* ------------------------------------------------------------------------ */
class MeshElementsRange {
public:
MeshElementsRange() = default;
MeshElementsRange(const elements_iterator & begin,
const elements_iterator & end)
: type((*begin).type), ghost_type((*begin).ghost_type), begin(begin),
end(end) {}
AKANTU_GET_MACRO(Type, type, const ElementType &);
AKANTU_GET_MACRO(GhostType, ghost_type, const GhostType &);
const Array & getElements() {
elements.resize(end - begin);
auto el_it = elements.begin();
for (auto it = begin; it != end; ++it, ++el_it) {
*el_it = it->element;
}
return elements;
}
private:
ElementType type{_not_defined};
GhostType ghost_type{_casper};
elements_iterator begin;
elements_iterator end;
Array elements;
};
/* ------------------------------------------------------------------------ */
class iterator {
struct element_comparator {
bool operator()(const Element & lhs, const Element & rhs) const {
bool res =
((lhs.kind < rhs.kind) ||
((lhs.kind == rhs.kind) && ((lhs.ghost_type < rhs.ghost_type) ||
((lhs.ghost_type == rhs.ghost_type) &&
((lhs.type < rhs.type))))));
return res;
}
};
public:
iterator(const iterator &) = default;
iterator(const elements_iterator & first, const elements_iterator & last)
: range(std::equal_range(first, last, *first, element_comparator())),
first(first), last(last) {}
decltype(auto) operator*() const {
return MeshElementsRange(range.first, range.second);
}
iterator operator++() {
first = range.second;
range = std::equal_range(first, last, *first, element_comparator());
return *this;
}
bool operator==(const iterator & other) const {
return (first == other.first and last == other.last);
}
bool operator!=(const iterator & other) const {
return (not operator==(other));
}
private:
std::pair range;
elements_iterator first;
elements_iterator last;
};
iterator begin() { return iterator(elements.begin(), elements.end()); }
iterator end() { return iterator(elements.end(), elements.end()); }
private:
Array elements;
};
/* -------------------------------------------------------------------------- */
namespace mesh_iterators {
namespace details {
template class delegated_iterator {
public:
+ using value_type = std::remove_pointer_t;
+ using pointer = value_type*;
+ using reference = value_type&;
+ using iterator_category = std::input_iterator_tag;
+
explicit delegated_iterator(internal_iterator it) : it(std::move(it)) {}
decltype(auto) operator*() {
return std::forwardsecond))>(*(it->second));
}
delegated_iterator operator++() {
++it;
return *this;
}
- bool operator==(const delegated_iterator & other) {
+ bool operator==(const delegated_iterator & other) const {
return other.it == it;
}
- bool operator!=(const delegated_iterator & other) {
+ bool operator!=(const delegated_iterator & other) const {
return other.it != it;
}
private:
internal_iterator it;
};
- template class MeshElementGroups {
+ template class ElementGroupsIterable {
public:
- explicit MeshElementGroups(M && mesh) : mesh(std::forward(mesh)) {}
+ explicit ElementGroupsIterable(GroupManager && group_manager)
+ : group_manager(std::forward(group_manager)) {}
+ size_t size() { return group_manager.getNbElementGroups(); }
decltype(auto) begin() {
- return delegated_iterator(
- mesh.element_group_begin());
+ return delegated_iterator(
+ group_manager.element_group_begin());
}
decltype(auto) begin() const {
- return delegated_iterator(
- mesh.element_group_begin());
+ return delegated_iterator(
+ group_manager.element_group_begin());
}
decltype(auto) end() {
- return delegated_iterator(
- mesh.element_group_end());
+ return delegated_iterator(
+ group_manager.element_group_end());
}
decltype(auto) end() const {
- return delegated_iterator(
- mesh.element_group_end());
+ return delegated_iterator(
+ group_manager.element_group_end());
}
private:
- M && mesh;
+ GroupManager && group_manager;
};
- template class MeshNodeGroups {
+ template class NodeGroupsIterable {
public:
- explicit MeshNodeGroups(M && mesh) : mesh(std::forward(mesh)) {}
+ explicit NodeGroupsIterable(GroupManager && group_manager)
+ : group_manager(std::forward(group_manager)) {}
+ size_t size() { return group_manager.getNbNodeGroups(); }
decltype(auto) begin() {
- return delegated_iterator(
- mesh.node_group_begin());
+ return delegated_iterator(
+ group_manager.node_group_begin());
}
decltype(auto) begin() const {
- return delegated_iterator(
- mesh.node_group_begin());
+ return delegated_iterator(
+ group_manager.node_group_begin());
}
decltype(auto) end() {
- return delegated_iterator(
- mesh.node_group_end());
+ return delegated_iterator(
+ group_manager.node_group_end());
}
decltype(auto) end() const {
- return delegated_iterator(
- mesh.node_group_end());
+ return delegated_iterator(
+ group_manager.node_group_end());
}
private:
- M && mesh;
+ GroupManager && group_manager;
};
} // namespace details
} // namespace mesh_iterators
-template decltype(auto) MeshElementGroups(Mesh && mesh) {
- return mesh_iterators::details::MeshElementGroups(mesh);
+template
+decltype(auto) ElementGroupsIterable(GroupManager && group_manager) {
+ return mesh_iterators::details::ElementGroupsIterable(group_manager);
}
-template decltype(auto) MeshNodeGroups(Mesh && mesh) {
- return mesh_iterators::details::MeshNodeGroups(mesh);
+template
+decltype(auto) NodeGroupsIterable(GroupManager && group_manager) {
+ return mesh_iterators::details::NodeGroupsIterable(group_manager);
}
} // namespace akantu
#endif /* __AKANTU_MESH_ITERATORS_HH__ */
diff --git a/src/mesh_utils/mesh_partition.cc b/src/mesh_utils/mesh_partition.cc
index 01ee4574f..fbb583db9 100644
--- a/src/mesh_utils/mesh_partition.cc
+++ b/src/mesh_utils/mesh_partition.cc
@@ -1,410 +1,417 @@
/**
* @file mesh_partition.cc
*
* @author David Simon Kammer
* @author Nicolas Richart
*
* @date creation: Tue Aug 17 2010
* @date last modification: Fri Jan 22 2016
*
* @brief implementation of common part of all partitioner
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh_partition.hh"
+#include "aka_iterators.hh"
#include "aka_types.hh"
#include "mesh_utils.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
MeshPartition::MeshPartition(const Mesh & mesh, UInt spatial_dimension,
const ID & id, const MemoryID & memory_id)
: Memory(id, memory_id), mesh(mesh), spatial_dimension(spatial_dimension),
partitions("partition", id, memory_id),
ghost_partitions("ghost_partition", id, memory_id),
ghost_partitions_offset("ghost_partition_offset", id, memory_id),
saved_connectivity("saved_connectivity", id, memory_id) {
AKANTU_DEBUG_IN();
+ Element elem;
+ elem.ghost_type = _not_ghost;
+
+ for (auto && type :
+ mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) {
+ elem.type = type;
+ for (auto && i : arange(mesh.getNbElement(elem.type))) {
+ elem.element = i;
+ lin_to_element.push_back(elem);
+ }
+ }
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MeshPartition::~MeshPartition() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* TODO this function should most probably be rewritten in a more modern way
* conversion in c++ of the GENDUALMETIS (mesh.c) function wrote by George in
* Metis (University of Minnesota)
*/
void MeshPartition::buildDualGraph(Array & dxadj, Array & dadjncy,
Array & edge_loads,
const EdgeLoadFunctor & edge_load_func) {
AKANTU_DEBUG_IN();
// tweak mesh;
UInt nb_good_types = 0;
std::vector nb_nodes_per_element_p1;
std::vector magic_number;
std::vector nb_element;
std::vector *> conn;
Element elem;
elem.ghost_type = _not_ghost;
UInt spatial_dimension = mesh.getSpatialDimension();
for (auto & type :
mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) {
elem.type = type;
ElementType type_p1 = Mesh::getP1ElementType(type);
conn.push_back(
&const_cast &>(mesh.getConnectivity(type, _not_ghost)));
nb_nodes_per_element_p1.push_back(Mesh::getNbNodesPerElement(type_p1));
nb_element.push_back(conn.back()->getSize());
magic_number.push_back(
Mesh::getNbNodesPerElement(Mesh::getFacetType(type_p1)));
- for (UInt i = 0; i < nb_element.back(); ++i) {
- elem.element = i;
- lin_to_element.push_back(elem);
- }
-
++nb_good_types;
}
CSR node_to_elem;
MeshUtils::buildNode2Elements(mesh, node_to_elem);
UInt nb_total_element = 0;
UInt nb_total_node_element = 0;
for (UInt t = 0; t < nb_good_types; ++t) {
nb_total_element += nb_element[t];
nb_total_node_element += nb_element[t] * nb_nodes_per_element_p1[t];
}
dxadj.resize(nb_total_element + 1);
/// initialize the dxadj array
for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t)
for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el)
dxadj(linerized_el) = nb_nodes_per_element_p1[t];
/// convert the dxadj_val array in a csr one
for (UInt i = 1; i < nb_total_element; ++i)
dxadj(i) += dxadj(i - 1);
for (UInt i = nb_total_element; i > 0; --i)
dxadj(i) = dxadj(i - 1);
dxadj(0) = 0;
dadjncy.resize(2 * dxadj(nb_total_element));
edge_loads.resize(2 * dxadj(nb_total_element));
/// weight map to determine adjacency
std::unordered_map weight_map;
for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) {
for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) {
/// fill the weight map
for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) {
UInt node = (*conn[t])(el, n);
CSR::iterator k;
for (k = node_to_elem.rbegin(node); k != node_to_elem.rend(node); --k) {
Element current_element = *k;
UInt current_el = lin_to_element.find(current_element);
if (current_el <= linerized_el)
break;
auto it_w = weight_map.find(current_el);
if (it_w == weight_map.end()) {
weight_map[current_el] = 1;
} else {
it_w->second++;
}
}
}
/// each element with a weight of the size of a facet are adjacent
for (auto it_w = weight_map.begin(); it_w != weight_map.end(); ++it_w) {
if (it_w->second == magic_number[t]) {
UInt adjacent_el = it_w->first;
#if defined(AKANTU_COHESIVE_ELEMENT)
/// Patch in order to prevent neighboring cohesive elements
/// from detecting each other
ElementKind linearized_el_kind = lin_to_element(linerized_el).kind;
ElementKind adjacent_el_kind = lin_to_element(adjacent_el).kind;
if (linearized_el_kind == adjacent_el_kind &&
linearized_el_kind == _ek_cohesive)
continue;
#endif
UInt index_adj = dxadj(adjacent_el)++;
UInt index_lin = dxadj(linerized_el)++;
dadjncy(index_lin) = adjacent_el;
dadjncy(index_adj) = linerized_el;
}
}
weight_map.clear();
}
}
Int k_start = 0;
for (UInt t = 0, linerized_el = 0, j = 0; t < nb_good_types; ++t)
for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) {
for (Int k = k_start; k < dxadj(linerized_el); ++k, ++j)
dadjncy(j) = dadjncy(k);
dxadj(linerized_el) = j;
k_start += nb_nodes_per_element_p1[t];
}
for (UInt i = nb_total_element; i > 0; --i)
dxadj(i) = dxadj(i - 1);
dxadj(0) = 0;
UInt adj = 0;
for (UInt i = 0; i < nb_total_element; ++i) {
UInt nb_adj = dxadj(i + 1) - dxadj(i);
for (UInt j = 0; j < nb_adj; ++j, ++adj) {
Int el_adj_id = dadjncy(dxadj(i) + j);
Element el = lin_to_element(i);
Element el_adj = lin_to_element(el_adj_id);
Int load = edge_load_func(el, el_adj);
edge_loads(adj) = load;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshPartition::fillPartitionInformation(
const Mesh & mesh, const Int * linearized_partitions) {
AKANTU_DEBUG_IN();
CSR node_to_elem;
MeshUtils::buildNode2Elements(mesh, node_to_elem);
UInt linearized_el = 0;
- for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) {
+ for (auto & type :
+ mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) {
UInt nb_element = mesh.getNbElement(type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
partitions.alloc(nb_element, 1, type, _not_ghost);
auto & ghost_part_csr = ghost_partitions_csr(type, _not_ghost);
ghost_part_csr.resizeRows(nb_element);
ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost);
ghost_partitions.alloc(0, 1, type, _ghost);
const Array & connectivity = mesh.getConnectivity(type, _not_ghost);
for (UInt el = 0; el < nb_element; ++el, ++linearized_el) {
UInt part = linearized_partitions[linearized_el];
partitions(type, _not_ghost)(el) = part;
std::list list_adj_part;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connectivity.storage()[el * nb_nodes_per_element + n];
CSR::iterator ne;
for (ne = node_to_elem.begin(node); ne != node_to_elem.end(node);
++ne) {
const Element & adj_element = *ne;
UInt adj_el = lin_to_element.find(adj_element);
UInt adj_part = linearized_partitions[adj_el];
if (part != adj_part) {
list_adj_part.push_back(adj_part);
}
}
}
list_adj_part.sort();
list_adj_part.unique();
for (auto & adj_part : list_adj_part) {
ghost_part_csr.getRows().push_back(adj_part);
ghost_part_csr.rowOffset(el)++;
ghost_partitions(type, _ghost).push_back(adj_part);
ghost_partitions_offset(type, _ghost)(el)++;
}
}
ghost_part_csr.countToCSR();
/// convert the ghost_partitions_offset array in an offset array
Array & ghost_partitions_offset_ptr =
ghost_partitions_offset(type, _ghost);
for (UInt i = 1; i < nb_element; ++i)
ghost_partitions_offset_ptr(i) += ghost_partitions_offset_ptr(i - 1);
for (UInt i = nb_element; i > 0; --i)
ghost_partitions_offset_ptr(i) = ghost_partitions_offset_ptr(i - 1);
ghost_partitions_offset_ptr(0) = 0;
}
// All Facets
for (Int sp = spatial_dimension - 1; sp >= 0; --sp) {
for (auto & type : mesh.elementTypes(sp, _not_ghost, _ek_not_defined)) {
UInt nb_element = mesh.getNbElement(type);
partitions.alloc(nb_element, 1, type, _not_ghost);
AKANTU_DEBUG_INFO("Allocating partitions for " << type);
auto & ghost_part_csr = ghost_partitions_csr(type, _not_ghost);
ghost_part_csr.resizeRows(nb_element);
ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost);
ghost_partitions.alloc(0, 1, type, _ghost);
AKANTU_DEBUG_INFO("Allocating ghost_partitions for " << type);
const Array> & elem_to_subelem =
mesh.getElementToSubelement(type, _not_ghost);
for (UInt i(0); i < mesh.getNbElement(type, _not_ghost);
++i) { // Facet loop
const std::vector & adjacent_elems = elem_to_subelem(i);
if (!adjacent_elems.empty()) {
Element min_elem;
UInt min_part(std::numeric_limits::max());
std::set adjacent_parts;
for (UInt j(0); j < adjacent_elems.size(); ++j) {
UInt adjacent_elem_id = adjacent_elems[j].element;
UInt adjacent_elem_part =
partitions(adjacent_elems[j].type,
adjacent_elems[j].ghost_type)(adjacent_elem_id);
if (adjacent_elem_part < min_part) {
min_part = adjacent_elem_part;
min_elem = adjacent_elems[j];
}
adjacent_parts.insert(adjacent_elem_part);
}
partitions(type, _not_ghost)(i) = min_part;
- auto git =
- ghost_partitions_csr(min_elem.type, _not_ghost)
- .begin(min_elem.element);
- auto gend =
- ghost_partitions_csr(min_elem.type, _not_ghost)
- .end(min_elem.element);
+ auto git = ghost_partitions_csr(min_elem.type, _not_ghost)
+ .begin(min_elem.element);
+ auto gend = ghost_partitions_csr(min_elem.type, _not_ghost)
+ .end(min_elem.element);
for (; git != gend; ++git) {
adjacent_parts.insert(min_part);
}
adjacent_parts.erase(min_part);
for (auto & part : adjacent_parts) {
ghost_part_csr.getRows().push_back(part);
ghost_part_csr.rowOffset(i)++;
ghost_partitions(type, _ghost).push_back(part);
}
ghost_partitions_offset(type, _ghost)(i + 1) =
ghost_partitions_offset(type, _ghost)(i + 1) +
adjacent_elems.size();
} else {
partitions(type, _not_ghost)(i) = 0;
}
}
ghost_part_csr.countToCSR();
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshPartition::tweakConnectivity(const Array & pairs) {
AKANTU_DEBUG_IN();
if (pairs.getSize() == 0)
return;
Mesh::type_iterator it =
mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined);
Mesh::type_iterator end =
mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined);
for (; it != end; ++it) {
ElementType type = *it;
Array & conn =
const_cast &>(mesh.getConnectivity(type, _not_ghost));
UInt nb_nodes_per_element = conn.getNbComponent();
UInt nb_element = conn.getSize();
Array & saved_conn = saved_connectivity.alloc(
nb_element, nb_nodes_per_element, type, _not_ghost);
saved_conn.copy(conn);
for (UInt i = 0; i < pairs.getSize(); ++i) {
for (UInt el = 0; el < nb_element; ++el) {
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
if (pairs(i, 1) == conn(el, n))
conn(el, n) = pairs(i, 0);
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshPartition::restoreConnectivity() {
AKANTU_DEBUG_IN();
ElementTypeMapArray::type_iterator it = saved_connectivity.firstType(
spatial_dimension, _not_ghost, _ek_not_defined);
ElementTypeMapArray::type_iterator end = saved_connectivity.lastType(
spatial_dimension, _not_ghost, _ek_not_defined);
for (; it != end; ++it) {
ElementType type = *it;
Array & conn =
const_cast &>(mesh.getConnectivity(type, _not_ghost));
Array & saved_conn = saved_connectivity(type, _not_ghost);
conn.copy(saved_conn);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc b/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc
index 7ba52fd3c..f5994e699 100644
--- a/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc
+++ b/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc
@@ -1,478 +1,480 @@
/**
* @file mesh_partition_scotch.cc
*
* @author David Simon Kammer
* @author Nicolas Richart
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Jan 22 2016
*
* @brief implementation of the MeshPartitionScotch class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*
*/
-
+/* -------------------------------------------------------------------------- */
+#include "mesh_partition_scotch.hh"
+#include "aka_common.hh"
+#include "aka_random_generator.hh"
+#include "aka_static_if.hh"
+#include "mesh_utils.hh"
/* -------------------------------------------------------------------------- */
#include
#include
-
/* -------------------------------------------------------------------------- */
-#include "aka_common.hh"
-#include "mesh_partition_scotch.hh"
-#include "mesh_utils.hh"
-/* -------------------------------------------------------------------------- */
#if !defined(AKANTU_USE_PTSCOTCH)
#ifndef AKANTU_SCOTCH_NO_EXTERN
extern "C" {
#endif // AKANTU_SCOTCH_NO_EXTERN
#include
#ifndef AKANTU_SCOTCH_NO_EXTERN
}
#endif // AKANTU_SCOTCH_NO_EXTERN
#else // AKANTU_USE_PTSCOTCH
#include
#endif // AKANTU_USE_PTSCOTCH
namespace akantu {
+namespace {
+ constexpr int scotch_version = int(SCOTCH_VERSION);
+}
+
/* -------------------------------------------------------------------------- */
MeshPartitionScotch::MeshPartitionScotch(const Mesh & mesh,
UInt spatial_dimension, const ID & id,
const MemoryID & memory_id)
: MeshPartition(mesh, spatial_dimension, id, memory_id) {
AKANTU_DEBUG_IN();
-// check if the akantu types and Scotch one are consistent
-#if __cplusplus > 199711L
- static_assert(sizeof(Int) == sizeof(SCOTCH_Num),
- "The integer type of Akantu does not match the one from Scotch");
-#else
- AKANTU_DEBUG_ASSERT(
+ // check if the akantu types and Scotch one are consistent
+ static_assert(
sizeof(Int) == sizeof(SCOTCH_Num),
"The integer type of Akantu does not match the one from Scotch");
-#endif
+
+ static_if(scotch_version >= 6)
+ .then([](auto && y) { SCOTCH_randomSeed(y); })
+ .else_([](auto && y) { srandom(y); })
+ (std::forward(RandomGenerator::seed()));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
static SCOTCH_Mesh * createMesh(const Mesh & mesh) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
UInt nb_nodes = mesh.getNbNodes();
UInt total_nb_element = 0;
UInt nb_edge = 0;
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator end = mesh.lastType(spatial_dimension);
for (; it != end; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
total_nb_element += nb_element;
nb_edge += nb_element * nb_nodes_per_element;
}
SCOTCH_Num vnodbas = 0;
SCOTCH_Num vnodnbr = nb_nodes;
SCOTCH_Num velmbas = vnodnbr;
SCOTCH_Num velmnbr = total_nb_element;
SCOTCH_Num * verttab = new SCOTCH_Num[vnodnbr + velmnbr + 1];
SCOTCH_Num * vendtab = verttab + 1;
SCOTCH_Num * velotab = NULL;
SCOTCH_Num * vnlotab = NULL;
SCOTCH_Num * vlbltab = NULL;
memset(verttab, 0, (vnodnbr + velmnbr + 1) * sizeof(SCOTCH_Num));
it = mesh.firstType(spatial_dimension);
for (; it != end; ++it) {
ElementType type = *it;
if (Mesh::getSpatialDimension(type) != spatial_dimension)
continue;
UInt nb_element = mesh.getNbElement(type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Array & connectivity = mesh.getConnectivity(type, _not_ghost);
/// count number of occurrence of each node
for (UInt el = 0; el < nb_element; ++el) {
UInt * conn_val = connectivity.storage() + el * nb_nodes_per_element;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
verttab[*(conn_val++)]++;
}
}
}
/// convert the occurrence array in a csr one
for (UInt i = 1; i < nb_nodes; ++i)
verttab[i] += verttab[i - 1];
for (UInt i = nb_nodes; i > 0; --i)
verttab[i] = verttab[i - 1];
verttab[0] = 0;
/// rearrange element to get the node-element list
SCOTCH_Num edgenbr = verttab[vnodnbr] + nb_edge;
SCOTCH_Num * edgetab = new SCOTCH_Num[edgenbr];
UInt linearized_el = 0;
it = mesh.firstType(spatial_dimension);
for (; it != end; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Array & connectivity = mesh.getConnectivity(type, _not_ghost);
for (UInt el = 0; el < nb_element; ++el, ++linearized_el) {
UInt * conn_val = connectivity.storage() + el * nb_nodes_per_element;
for (UInt n = 0; n < nb_nodes_per_element; ++n)
edgetab[verttab[*(conn_val++)]++] = linearized_el + velmbas;
}
}
for (UInt i = nb_nodes; i > 0; --i)
verttab[i] = verttab[i - 1];
verttab[0] = 0;
SCOTCH_Num * verttab_tmp = verttab + vnodnbr + 1;
SCOTCH_Num * edgetab_tmp = edgetab + verttab[vnodnbr];
it = mesh.firstType(spatial_dimension);
for (; it != end; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Array & connectivity = mesh.getConnectivity(type, _not_ghost);
for (UInt el = 0; el < nb_element; ++el) {
*verttab_tmp = *(verttab_tmp - 1) + nb_nodes_per_element;
verttab_tmp++;
UInt * conn = connectivity.storage() + el * nb_nodes_per_element;
for (UInt i = 0; i < nb_nodes_per_element; ++i) {
*(edgetab_tmp++) = *(conn++) + vnodbas;
}
}
}
SCOTCH_Mesh * meshptr = new SCOTCH_Mesh;
SCOTCH_meshInit(meshptr);
SCOTCH_meshBuild(meshptr, velmbas, vnodbas, velmnbr, vnodnbr, verttab,
vendtab, velotab, vnlotab, vlbltab, edgenbr, edgetab);
/// Check the mesh
AKANTU_DEBUG_ASSERT(SCOTCH_meshCheck(meshptr) == 0,
"Scotch mesh is not consistent");
#ifndef AKANTU_NDEBUG
if (AKANTU_DEBUG_TEST(dblDump)) {
/// save initial graph
FILE * fmesh = fopen("ScotchMesh.msh", "w");
SCOTCH_meshSave(meshptr, fmesh);
fclose(fmesh);
/// write geometry file
std::ofstream fgeominit;
fgeominit.open("ScotchMesh.xyz");
- fgeominit << spatial_dimension << std::endl
- << nb_nodes << std::endl;
+ fgeominit << spatial_dimension << std::endl << nb_nodes << std::endl;
const Array & nodes = mesh.getNodes();
Real * nodes_val = nodes.storage();
for (UInt i = 0; i < nb_nodes; ++i) {
fgeominit << i << " ";
for (UInt s = 0; s < spatial_dimension; ++s)
fgeominit << *(nodes_val++) << " ";
fgeominit << std::endl;
;
}
fgeominit.close();
}
#endif
AKANTU_DEBUG_OUT();
return meshptr;
}
/* -------------------------------------------------------------------------- */
static void destroyMesh(SCOTCH_Mesh * meshptr) {
AKANTU_DEBUG_IN();
SCOTCH_Num velmbas, vnodbas, vnodnbr, velmnbr, *verttab, *vendtab, *velotab,
*vnlotab, *vlbltab, edgenbr, *edgetab, degrptr;
SCOTCH_meshData(meshptr, &velmbas, &vnodbas, &velmnbr, &vnodnbr, &verttab,
&vendtab, &velotab, &vnlotab, &vlbltab, &edgenbr, &edgetab,
°rptr);
delete[] verttab;
delete[] edgetab;
SCOTCH_meshExit(meshptr);
delete meshptr;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshPartitionScotch::partitionate(UInt nb_part,
const EdgeLoadFunctor & edge_load_func,
const Array & pairs) {
AKANTU_DEBUG_IN();
nb_partitions = nb_part;
tweakConnectivity(pairs);
AKANTU_DEBUG_INFO("Partitioning the mesh " << mesh.getID() << " in "
<< nb_part << " parts.");
Array dxadj;
Array dadjncy;
Array edge_loads;
buildDualGraph(dxadj, dadjncy, edge_loads, edge_load_func);
/// variables that will hold our structures in scotch format
SCOTCH_Graph scotch_graph;
SCOTCH_Strat scotch_strat;
/// description number and arrays for struct mesh for scotch
SCOTCH_Num baseval = 0; // base numbering for element and
// nodes (0 -> C , 1 -> fortran)
SCOTCH_Num vertnbr = dxadj.getSize() - 1; // number of vertexes
SCOTCH_Num * parttab; // array of partitions
SCOTCH_Num edgenbr = dxadj(vertnbr); // twice the number of "edges"
//(an "edge" bounds two nodes)
SCOTCH_Num * verttab = dxadj.storage(); // array of start indices in edgetab
SCOTCH_Num * vendtab = NULL; // array of after-last indices in edgetab
SCOTCH_Num * velotab = NULL; // integer load associated with
// every vertex ( optional )
SCOTCH_Num * edlotab = edge_loads.storage(); // integer load associated with
// every edge ( optional )
SCOTCH_Num * edgetab = dadjncy.storage(); // adjacency array of every vertex
SCOTCH_Num * vlbltab = NULL; // vertex label array (optional)
/// Allocate space for Scotch arrays
parttab = new SCOTCH_Num[vertnbr];
/// Initialize the strategy structure
SCOTCH_stratInit(&scotch_strat);
/// Initialize the graph structure
SCOTCH_graphInit(&scotch_graph);
/// Build the graph from the adjacency arrays
SCOTCH_graphBuild(&scotch_graph, baseval, vertnbr, verttab, vendtab, velotab,
vlbltab, edgenbr, edgetab, edlotab);
#ifndef AKANTU_NDEBUG
if (AKANTU_DEBUG_TEST(dblDump)) {
/// save initial graph
FILE * fgraphinit = fopen("GraphIniFile.grf", "w");
SCOTCH_graphSave(&scotch_graph, fgraphinit);
fclose(fgraphinit);
/// write geometry file
std::ofstream fgeominit;
fgeominit.open("GeomIniFile.xyz");
- fgeominit << spatial_dimension << std::endl
- << vertnbr << std::endl;
+ fgeominit << spatial_dimension << std::endl << vertnbr << std::endl;
const Array & nodes = mesh.getNodes();
Mesh::type_iterator f_it =
mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined);
Mesh::type_iterator f_end =
mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined);
Array::const_vector_iterator nodes_it =
nodes.begin(spatial_dimension);
UInt out_linerized_el = 0;
for (; f_it != f_end; ++f_it) {
ElementType type = *f_it;
UInt nb_element = mesh.getNbElement(*f_it);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Array & connectivity = mesh.getConnectivity(type);
Vector mid(spatial_dimension);
for (UInt el = 0; el < nb_element; ++el) {
mid.set(0.);
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connectivity.storage()[nb_nodes_per_element * el + n];
mid += Vector(nodes_it[node]);
}
mid /= nb_nodes_per_element;
fgeominit << out_linerized_el++ << " ";
for (UInt s = 0; s < spatial_dimension; ++s)
fgeominit << mid[s] << " ";
fgeominit << std::endl;
;
}
}
fgeominit.close();
}
#endif
/// Check the graph
AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0,
"Graph to partition is not consistent");
/// Partition the mesh
SCOTCH_graphPart(&scotch_graph, nb_part, &scotch_strat, parttab);
/// Check the graph
AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0,
"Partitioned graph is not consistent");
#ifndef AKANTU_NDEBUG
if (AKANTU_DEBUG_TEST(dblDump)) {
/// save the partitioned graph
FILE * fgraph = fopen("GraphFile.grf", "w");
SCOTCH_graphSave(&scotch_graph, fgraph);
fclose(fgraph);
/// save the partition map
std::ofstream fmap;
fmap.open("MapFile.map");
fmap << vertnbr << std::endl;
for (Int i = 0; i < vertnbr; i++)
fmap << i << " " << parttab[i] << std::endl;
fmap.close();
}
#endif
/// free the scotch data structures
SCOTCH_stratExit(&scotch_strat);
SCOTCH_graphFree(&scotch_graph);
SCOTCH_graphExit(&scotch_graph);
fillPartitionInformation(mesh, parttab);
delete[] parttab;
restoreConnectivity();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshPartitionScotch::reorder() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Reordering the mesh " << mesh.getID());
SCOTCH_Mesh * scotch_mesh = createMesh(mesh);
UInt nb_nodes = mesh.getNbNodes();
SCOTCH_Strat scotch_strat;
// SCOTCH_Ordering scotch_order;
SCOTCH_Num * permtab = new SCOTCH_Num[nb_nodes];
SCOTCH_Num * peritab = NULL;
SCOTCH_Num cblknbr = 0;
SCOTCH_Num * rangtab = NULL;
SCOTCH_Num * treetab = NULL;
/// Initialize the strategy structure
SCOTCH_stratInit(&scotch_strat);
SCOTCH_Graph scotch_graph;
SCOTCH_graphInit(&scotch_graph);
SCOTCH_meshGraph(scotch_mesh, &scotch_graph);
#ifndef AKANTU_NDEBUG
if (AKANTU_DEBUG_TEST(dblDump)) {
FILE * fgraphinit = fopen("ScotchMesh.grf", "w");
SCOTCH_graphSave(&scotch_graph, fgraphinit);
fclose(fgraphinit);
}
#endif
/// Check the graph
// AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0,
// "Mesh to Graph is not consistent");
SCOTCH_graphOrder(&scotch_graph, &scotch_strat, permtab, peritab, &cblknbr,
rangtab, treetab);
SCOTCH_graphExit(&scotch_graph);
SCOTCH_stratExit(&scotch_strat);
destroyMesh(scotch_mesh);
/// Renumbering
UInt spatial_dimension = mesh.getSpatialDimension();
for (UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType)g;
Mesh::type_iterator it = mesh.firstType(_all_dimensions, gt);
Mesh::type_iterator end = mesh.lastType(_all_dimensions, gt);
for (; it != end; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type, gt);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Array & connectivity = mesh.getConnectivity(type, gt);
UInt * conn = connectivity.storage();
for (UInt el = 0; el < nb_element * nb_nodes_per_element; ++el, ++conn) {
*conn = permtab[*conn];
}
}
}
/// \todo think of a in-place way to do it
Real * new_coordinates = new Real[spatial_dimension * nb_nodes];
Real * old_coordinates = mesh.getNodes().storage();
for (UInt i = 0; i < nb_nodes; ++i) {
memcpy(new_coordinates + permtab[i] * spatial_dimension,
old_coordinates + i * spatial_dimension,
spatial_dimension * sizeof(Real));
}
memcpy(old_coordinates, new_coordinates,
nb_nodes * spatial_dimension * sizeof(Real));
delete[] new_coordinates;
delete[] permtab;
AKANTU_DEBUG_OUT();
}
-} // akantu
+} // namespace akantu
diff --git a/test/test_synchronizer/test_data_distribution.cc b/test/test_synchronizer/test_data_distribution.cc
index 2efd1498e..cc042b24d 100644
--- a/test/test_synchronizer/test_data_distribution.cc
+++ b/test/test_synchronizer/test_data_distribution.cc
@@ -1,197 +1,173 @@
/**
* @file test_data_distribution.cc
*
* @author Nicolas Richart
*
* @date creation: Fri Sep 05 2014
* @date last modification: Sun Oct 19 2014
*
* @brief Test the mesh distribution on creation of a distributed synchonizer
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 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 "aka_random_generator.hh"
#include "element_group.hh"
#include "element_synchronizer.hh"
+#include "mesh_iterators.hh"
#include "mesh_partition_mesh_data.hh"
-#include "aka_random_generator.hh"
-
+/* -------------------------------------------------------------------------- */
+#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize(argc, argv);
const UInt spatial_dimension = 3;
Mesh mesh_group_after(spatial_dimension, "after");
Mesh mesh_group_before(spatial_dimension, "before");
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int psize = comm.getNbProc();
Int prank = comm.whoAmI();
if (prank == 0) {
mesh_group_before.read("data_split.msh");
mesh_group_after.read("data_split.msh");
mesh_group_before.registerData("global_id");
mesh_group_after.registerData("global_id");
- for (Mesh::type_iterator tit = mesh_group_after.firstType(_all_dimensions);
- tit != mesh_group_after.lastType(_all_dimensions); ++tit) {
- Array & gidb =
- mesh_group_before.getDataPointer("global_id", *tit);
- Array & gida =
- mesh_group_after.getDataPointer("global_id", *tit);
-
- Array::scalar_iterator ait = gida.begin();
- Array::scalar_iterator bit = gidb.begin();
- Array::scalar_iterator end = gida.end();
- for (UInt i = 0; ait != end; ++ait, ++i, ++bit) {
- *bit = i;
- *ait = i;
+ for (const auto & type : mesh_group_after.elementTypes(_all_dimensions)) {
+ auto & gidb = mesh_group_before.getDataPointer("global_id", type);
+ auto & gida = mesh_group_after.getDataPointer("global_id", type);
+
+ for (auto && data : zip(arange(gida.getSize()), gida, gidb)) {
+ std::get<1>(data) = std::get<0>(data);
+ std::get<2>(data) = std::get<0>(data);
}
}
}
RandomGenerator::seed(1);
mesh_group_before.distribute();
RandomGenerator::seed(1);
mesh_group_after.distribute();
if (prank == 0)
std::cout << mesh_group_after;
- GroupManager::element_group_iterator grp_ait =
- mesh_group_after.element_group_begin();
- GroupManager::element_group_iterator grp_end =
- mesh_group_after.element_group_end();
- for (; grp_ait != grp_end; ++grp_ait) {
- std::string grp = grp_ait->first;
- const ElementGroup & bgrp = mesh_group_before.getElementGroup(grp);
- const ElementGroup & agrp = *grp_ait->second;
-
- for (ghost_type_t::iterator git = ghost_type_t::begin();
- git != ghost_type_t::end(); ++git) {
- GhostType ghost_type = *git;
-
- for (Mesh::type_iterator tit =
- bgrp.firstType(_all_dimensions, ghost_type);
- tit != bgrp.lastType(_all_dimensions, ghost_type); ++tit) {
+ for (const auto & agrp : ElementGroupsIterable(mesh_group_after)) {
+ const ElementGroup & bgrp =
+ mesh_group_before.getElementGroup(agrp.getName());
+
+ for (auto && ghost_type : ghost_types) {
+ for (const auto & type : bgrp.elementTypes(_all_dimensions, ghost_type)) {
Array & gidb = mesh_group_before.getDataPointer(
- "global_id", *tit, ghost_type);
+ "global_id", type, ghost_type);
Array & gida = mesh_group_after.getDataPointer(
- "global_id", *tit, ghost_type);
+ "global_id", type, ghost_type);
- Array bgelem(bgrp.getElements(*tit, ghost_type));
- Array agelem(agrp.getElements(*tit, ghost_type));
+ Array bgelem(bgrp.getElements(type, ghost_type));
+ Array agelem(agrp.getElements(type, ghost_type));
- Array::scalar_iterator ait = agelem.begin();
- Array