diff --git a/cmake/libakantu/v3/printers.py b/cmake/libakantu/v3/printers.py
index de6356340..a8ca32bee 100755
--- a/cmake/libakantu/v3/printers.py
+++ b/cmake/libakantu/v3/printers.py
@@ -1,234 +1,339 @@
#!/usr/bin/env python3
# encoding: utf-8
#
# Inspired from boost's pretty printers from
# Rüdiger Sonderfeld <ruediger@c-plusplus.de>
# and from Pretty-printers for libstc++ from Free Software Foundation, Inc.
#
-
import gdb
import re
-#import sys
-# import libstdcxx.v6.printers as std
+import sys
__use_gdb_pp__ = True
try:
import gdb.printing
except ImportError:
__use_gdb_pp__ = False
class AkantuPrinter(object):
regex = None
@classmethod
def supports(cls, typename):
#print('{0} ~= {1}'.format(typename, cls.regex), file=sys.stderr)
return cls.regex.search(typename)
@classmethod
def get_basic_type(cls, value):
""" Determines the type associated to a value"""
_type = value.type
# If it points to a reference, get the reference.
if _type.code == gdb.TYPE_CODE_REF:
_type = _type.target()
# Get the unqualified type, stripped of typedefs.
_type = _type.unqualified().strip_typedefs()
return _type.tag
if __use_gdb_pp__:
__akantu_pretty_printers__ = \
gdb.printing.RegexpCollectionPrettyPrinter("libakantu-v3")
else:
class AkantuPrettyPrinters(object):
def __init__(self, name):
super(AkantuPrettyPrinters, self).__init__()
self.name = name
self.printers = {}
@property
def enabled(self):
return True
def add_printer(self, name, regex, printer):
self.printers[name] = printer
def __call__(self, val):
typename = AkantuPrinter.get_basic_type(val)
if not typename:
return None
for name, printer in self.printers.items():
if(printer.supports(typename)):
return printer
return None
__akantu_pretty_printers__ = AkantuPrettyPrinters("libakantu-v3")
def register_pretty_printer(pretty_printer):
"Registers a Pretty Printer"
__akantu_pretty_printers__.add_printer(pretty_printer.name,
pretty_printer.regex,
pretty_printer)
return pretty_printer
@register_pretty_printer
class AkaArrayPrinter(AkantuPrinter):
"""Pretty printer for akantu::Array<T>"""
regex = re.compile('^akantu::Array<(.*?), (true|false)>$')
name = 'akantu::Array'
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.ptr = self.value['values']
self.size = int(self.value['size_'])
self.nb_component = int(self.value['nb_component'])
def display_hint(self):
return 'array'
def to_string(self):
m = self.regex.search(self.typename)
return 'Array<{0}>({1}, {2}) stored at {3}'.format(
m.group(1), self.size, self.nb_component, self.ptr)
def children(self):
_ptr = self.ptr
for i in range(self.size):
_values = ["{0}".format((_ptr + j).dereference())
for j in range(self.nb_component)]
_ptr = _ptr + self.nb_component
yield ('[{0}]'.format(i),
('{0}' if self.nb_component == 1 else '[{0}]').format(
', '.join(_values)))
-# @register_pretty_printer
-# class AkaArrayIteratorPrinter(AkantuPrinter):
-# """Pretty printer for akantu::Array<T>"""
-# regex = re.compile('^akantu::Array<(.*?), (true|false)>::internal_iterator<(.*?), (.*?), (false|true)>$')
-# name = 'akantu::Array::iterator'
-
-# def __init__(self, value):
-# self.typename = self.get_basic_type(value)
-# self.value = value
-# self.ret = self.value['ret'].dereference()
-
-# def to_string(self):
-# m = self.regex.search(self.typename)
-# return 'Array<{0}>::iterator<{3}>'.format(
-# m.group(1), m.group(1))
-
-# def children(self):
-# yield ('[data]', self.ret)
+if sys.version_info[0] > 2:
+ ### Python 3 stuff
+ Iterator = object
+else:
+ ### Python 2 stuff
+ class Iterator:
+ """Compatibility mixin for iterators
+
+ Instead of writing next() methods for iterators, write
+ __next__() methods and use this mixin to make them work in
+ Python 2 as well as Python 3.
+
+ Idea stolen from the "six" documentation:
+ <http://pythonhosted.org/six/#six.Iterator>
+ """
+ def next(self):
+ return self.__next__()
+
+
+class RbtreeIterator(Iterator):
+ """
+ Turn an RB-tree-based container (std::map, std::set etc.) into
+ a Python iterable object.
+ """
+
+ def __init__(self, rbtree):
+ self.size = rbtree['_M_t']['_M_impl']['_M_node_count']
+ self.node = rbtree['_M_t']['_M_impl']['_M_header']['_M_left']
+ self.count = 0
+
+ def __iter__(self):
+ return self
+
+ def __len__(self):
+ return int (self.size)
+
+ def __next__(self):
+ if self.count == self.size:
+ raise StopIteration
+ result = self.node
+ self.count = self.count + 1
+ if self.count < self.size:
+ # Compute the next node.
+ node = self.node
+ if node.dereference()['_M_right']:
+ node = node.dereference()['_M_right']
+ while node.dereference()['_M_left']:
+ node = node.dereference()['_M_left']
+ else:
+ parent = node.dereference()['_M_parent']
+ while node == parent.dereference()['_M_right']:
+ node = parent
+ parent = parent.dereference()['_M_parent']
+ if node.dereference()['_M_right'] != parent:
+ node = parent
+ self.node = node
+ return result
+
+def get_value_from_Rb_tree_node(node):
+ """Returns the value held in an _Rb_tree_node<_Val>"""
+ try:
+ member = node.type.fields()[1].name
+ if member == '_M_value_field':
+ # C++03 implementation, node contains the value as a member
+ return node['_M_value_field']
+ elif member == '_M_storage':
+ # C++11 implementation, node stores value in __aligned_membuf
+ valtype = node.type.template_argument(0)
+ return get_value_from_aligned_membuf(node['_M_storage'], valtype)
+ except:
+ pass
+ raise ValueError("Unsupported implementation for %s" % str(node.type))
+
+@register_pretty_printer
+class AkaElementTypeMapArrayPrinter(AkantuPrinter):
+ """Pretty printer for akantu::ElementTypeMap<Array<T>>"""
+ regex = re.compile('^akantu::ElementTypeMap<akantu::Array<(.*?), (true|false)>\*, akantu::(.*?)>$')
+ name = 'akantu::ElementTypeMapArray'
+
+ # Turn an RbtreeIterator into a pretty-print iterator.
+ class _iter(Iterator):
+ def __init__(self, rbiter, type):
+ self.rbiter = rbiter
+ self.count = 0
+ self.type = type
+
+ def __iter__(self):
+ return self
+
+ def __next__(self):
+ if self.count % 2 == 0:
+ n = next(self.rbiter)
+ n = n.cast(self.type).dereference()
+ n = get_value_from_Rb_tree_node(n)
+ self.pair = n
+ item = n['first']
+ else:
+ item = self.pair['second']
+ result = ('[{0}]'.format(self.count), item.dereference())
+ self.count = self.count + 1
+ return result
+
+ def __init__(self, value):
+ self.typename = self.get_basic_type(value)
+ self.value = value
+ self.data = self.value['data']
+ self.ghost_data = self.value['ghost_data']
+
+ def to_string(self):
+ m = self.regex.search(self.typename)
+ return 'ElementTypMapArray<{0}> with {1} _not_ghost and {2} _ghost'.format(
+ m.group(1), len(RbtreeIterator(self.data)), len(RbtreeIterator(self.ghost_data)))
+
+ def children(self):
+ m = self.regex.search(self.typename)
+ try:
+ _type = gdb.lookup_type("akantu::Array<{0}, {1}>".format(
+ m.group(1), m.group(2))).strip_typedefs().pointer()
+ yield ('[_not_ghost]', self._iter(RbtreeIterator(self.data), _type))
+ yield ('[_ghost]', self._iter(RbtreeIterator(self.ghost_data), _type))
+ except RuntimeError:
+ pass
+
+ def display_hint(self):
+ return 'map'
# @register_pretty_printer
class AkaTensorPrinter(AkantuPrinter):
"""Pretty printer for akantu::Tensor<T>"""
regex = re.compile('^akantu::Tensor<(.*), +(.*), +(.*)>$')
name = 'akantu::Tensor'
value = None
typename = ""
ptr = None
dims = []
ndims = 0
def pretty_print(self):
def ij2str(i, j, m):
return "{0}".format((self.ptr+m*j + i).dereference())
def line(i, m, n):
return "[{0}]".format(", ".join((ij2str(i, j, m) for j in
range(n))))
m = int(self.dims[0])
if (self.ndims == 1):
n = 1
else:
n = int(self.dims[1])
return "[{0}]".format(", ".join(line(i, m, n) for i in range(m)))
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.ptr = self.value['values']
self.dims = self.value['n']
def children(self):
yield ('values', self.pretty_print())
yield ('wrapped', self.value['wrapped'])
@register_pretty_printer
class AkaVectorPrinter(AkaTensorPrinter):
"""Pretty printer for akantu::Vector<T>"""
regex = re.compile('^akantu::Vector<(.*)>$')
name = 'akantu::Vector'
n = 0
ptr = 0x0
def __init__(self, value):
super(AkaVectorPrinter, self).__init__(value)
self.ndims = 1
def to_string(self):
m = self.regex.search(self.typename)
return 'Vector<{0}>({1}) [{2}]'.format(m.group(1), int(self.dims[0]),
str(self.ptr))
@register_pretty_printer
class AkaMatrixPrinter(AkaTensorPrinter):
"""Pretty printer for akantu::Matrix<T>"""
regex = re.compile('^akantu::Matrix<(.*)>$')
name = 'akantu::Matrix'
def __init__(self, value):
super(AkaMatrixPrinter, self).__init__(value)
self.ndims = 2
def to_string(self):
m = self.regex.search(self.typename)
return 'Matrix<%s>(%d, %d) [%s]' % (m.group(1), int(self.dims[0]),
int(self.dims[1]),
str(self.ptr))
@register_pretty_printer
class AkaElementPrinter(AkantuPrinter):
"""Pretty printer for akantu::Element"""
regex = re.compile('^akantu::Element$')
name = 'akantu::Element'
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.element = self.value['element']
self.eltype = self.value['type']
self.ghost_type = self.value['ghost_type']
def to_string(self):
return 'Element({0}, {1}, {2})'.format(self.element, self.eltype,
self.ghost_type)
def register_akantu_printers(obj):
"Register Akantu Pretty Printers."
if __use_gdb_pp__:
gdb.printing.register_pretty_printer(obj, __akantu_pretty_printers__)
else:
if obj is None:
obj = gdb
obj.pretty_printers.append(__akantu_pretty_printers__)
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index a5c46d45c..ac1b92fd9 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1260 +1,1263 @@
/**
* @file aka_array_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Jul 15 2010
* @date last modification: Fri Jan 22 2016
*
* @brief Inline functions of the classes Array<T> and ArrayBase
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* Inline Functions Array<T> */
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <memory>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__
#define __AKANTU_AKA_ARRAY_TMPL_HH__
namespace akantu {
namespace debug {
struct ArrayException : public Exception {};
} // namespace debug
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline T & Array<T, is_scal>::operator()(UInt i, UInt j) {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
"The value at position ["
<< i << "," << j << "] is out of range in array \""
<< id << "\"");
return values[i * nb_component + j];
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline const T & Array<T, is_scal>::operator()(UInt i, UInt j) const {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
"The value at position ["
<< i << "," << j << "] is out of range in array \""
<< id << "\"");
return values[i * nb_component + j];
}
template <class T, bool is_scal>
inline T & Array<T, is_scal>::operator[](UInt i) {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
"The value at position ["
<< i << "] is out of range in array \"" << id
<< "\"");
return values[i];
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline const T & Array<T, is_scal>::operator[](UInt i) const {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
"The value at position ["
<< i << "] is out of range in array \"" << id
<< "\"");
return values[i];
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array with the value value for all components
* @param value the new last tuple or the array will contain nb_component copies
* of value
*/
template <class T, bool is_scal>
inline void Array<T, is_scal>::push_back(const T & value) {
resizeUnitialized(size_ + 1, true, value);
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array
* @param new_elem a C-array containing the values to be copied to the end of
* the array */
// template <class T, bool is_scal>
// inline void Array<T, is_scal>::push_back(const T new_elem[]) {
// UInt pos = size_;
// resizeUnitialized(size_ + 1, false);
// T * tmp = values + nb_component * pos;
// std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
// }
/* -------------------------------------------------------------------------- */
#ifndef SWIG
/**
* append a matrix or a vector to the array
* @param new_elem a reference to a Matrix<T> or Vector<T> */
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline void Array<T, is_scal>::push_back(const C<T> & new_elem) {
AKANTU_DEBUG_ASSERT(
nb_component == new_elem.size(),
"The vector("
<< new_elem.size()
<< ") as not a size compatible with the Array (nb_component="
<< nb_component << ").");
UInt pos = size_;
resizeUnitialized(size_ + 1, false);
T * tmp = values + nb_component * pos;
std::uninitialized_copy(new_elem.storage(), new_elem.storage() + nb_component,
tmp);
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array
* @param it an iterator to the tuple to be copied to the end of the array */
template <class T, bool is_scal>
template <class Ret>
inline void
Array<T, is_scal>::push_back(const Array<T, is_scal>::iterator<Ret> & it) {
UInt pos = size_;
resizeUnitialized(size_ + 1, false);
T * tmp = values + nb_component * pos;
T * new_elem = it.data();
std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
}
#endif
/* -------------------------------------------------------------------------- */
/**
* erase an element. If the erased element is not the last of the array, the
* last element is moved into the hole in order to maintain contiguity. This
* may invalidate existing iterators (For instance an iterator obtained by
* Array::end() is no longer correct) and will change the order of the
* elements.
* @param i index of element to erase
*/
template <class T, bool is_scal> inline void Array<T, is_scal>::erase(UInt i) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT((size_ > 0), "The array is empty");
AKANTU_DEBUG_ASSERT((i < size_),
"The element at position [" << i << "] is out of range ("
<< i << ">=" << size_ << ")");
if (i != (size_ - 1)) {
for (UInt j = 0; j < nb_component; ++j) {
values[i * nb_component + j] = values[(size_ - 1) * nb_component + j];
}
}
resize(size_ - 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Subtract another array entry by entry from this array in place. Both arrays
* must
* have the same size and nb_component. If the arrays have different shapes,
* code compiled in debug mode will throw an expeption and optimised code
* will behave in an unpredicted manner
* @param other array to subtract from this
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator-=(const Array<T, is_scal> & vect) {
AKANTU_DEBUG_ASSERT((size_ == vect.size_) &&
(nb_component == vect.nb_component),
"The too array don't have the same sizes");
T * a = values;
T * b = vect.storage();
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a -= *b;
++a;
++b;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Add another array entry by entry to this array in place. Both arrays must
* have the same size and nb_component. If the arrays have different shapes,
* code compiled in debug mode will throw an expeption and optimised code
* will behave in an unpredicted manner
* @param other array to add to this
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator+=(const Array<T, is_scal> & vect) {
AKANTU_DEBUG_ASSERT((size_ == vect.size) &&
(nb_component == vect.nb_component),
"The too array don't have the same sizes");
T * a = values;
T * b = vect.storage();
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a++ += *b++;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Multiply all entries of this array by a scalar in place
* @param alpha scalar multiplicant
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::operator*=(const T & alpha) {
T * a = values;
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a++ *= alpha;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Compare this array element by element to another.
* @param other array to compare to
* @return true it all element are equal and arrays have the same shape, else
* false
*/
template <class T, bool is_scal>
bool Array<T, is_scal>::operator==(const Array<T, is_scal> & array) const {
bool equal = nb_component == array.nb_component && size_ == array.size_ &&
id == array.id;
if (!equal)
return false;
if (values == array.storage())
return true;
else
return std::equal(values, values + size_ * nb_component, array.storage());
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
bool Array<T, is_scal>::operator!=(const Array<T, is_scal> & array) const {
return !operator==(array);
}
/* -------------------------------------------------------------------------- */
#ifndef SWIG
/**
* set all tuples of the array to a given vector or matrix
* @param vm Matrix or Vector to fill the array with
*/
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline void Array<T, is_scal>::set(const C<T> & vm) {
AKANTU_DEBUG_ASSERT(
nb_component == vm.size(),
"The size of the object does not match the number of components");
for (T * it = values; it < values + nb_component * size_;
it += nb_component) {
std::copy(vm.storage(), vm.storage() + nb_component, it);
}
}
#endif
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::append(const Array<T> & other) {
AKANTU_DEBUG_ASSERT(
nb_component == other.nb_component,
"Cannot append an array with a different number of component");
UInt old_size = this->size_;
this->resizeUnitialized(this->size_ + other.size(), false);
T * tmp = values + nb_component * old_size;
std::uninitialized_copy(other.storage(),
other.storage() + other.size() * nb_component, tmp);
}
/* -------------------------------------------------------------------------- */
/* Functions Array<T, is_scal> */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const ID & id)
: ArrayBase(id), values(nullptr) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
if (!is_scal) {
T val = T();
std::uninitialized_fill(values, values + size * nb_component, val);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const T def_values[],
const ID & id)
: ArrayBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
T * tmp = values;
for (UInt i = 0; i < size; ++i) {
tmp = values + nb_component * i;
std::uninitialized_copy(def_values, def_values + nb_component, tmp);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const T & value,
const ID & id)
: ArrayBase(id), values(nullptr) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
std::uninitialized_fill_n(values, size * nb_component, value);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(const Array<T, is_scal> & vect, bool deep,
const ID & id)
: ArrayBase(vect) {
AKANTU_DEBUG_IN();
this->id = (id == "") ? vect.id : id;
if (deep) {
allocate(vect.size_, vect.nb_component);
T * tmp = values;
std::uninitialized_copy(vect.storage(),
vect.storage() + size_ * nb_component, tmp);
} else {
this->values = vect.storage();
this->size_ = vect.size_;
this->nb_component = vect.nb_component;
this->allocated_size = vect.allocated_size;
this->size_of_type = vect.size_of_type;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
#ifndef SWIG
template <class T, bool is_scal>
Array<T, is_scal>::Array(const std::vector<T> & vect) {
AKANTU_DEBUG_IN();
this->id = "";
allocate(vect.size(), 1);
T * tmp = values;
std::uninitialized_copy(&(vect[0]), &(vect[size_ - 1]), tmp);
AKANTU_DEBUG_OUT();
}
#endif
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> Array<T, is_scal>::~Array() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG(dblAccessory,
"Freeing " << printMemorySize<T>(allocated_size * nb_component)
<< " (" << id << ")");
if (values) {
if (!is_scal)
for (UInt i = 0; i < size_ * nb_component; ++i) {
T * obj = values + i;
obj->~T();
}
free(values);
}
size_ = allocated_size = 0;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* perform the allocation for the constructors
* @param size is the size of the array
* @param nb_component is the number of component of the array
*/
template <class T, bool is_scal>
void Array<T, is_scal>::allocate(UInt size, UInt nb_component) {
AKANTU_DEBUG_IN();
if (size == 0) {
values = nullptr;
} else {
values = static_cast<T *>(malloc(nb_component * size * sizeof(T)));
AKANTU_DEBUG_ASSERT(values != nullptr,
"Cannot allocate "
<< printMemorySize<T>(size * nb_component) << " ("
<< id << ")");
}
if (values == nullptr) {
this->size_ = this->allocated_size = 0;
} else {
AKANTU_DEBUG(dblAccessory,
"Allocated " << printMemorySize<T>(size * nb_component) << " ("
<< id << ")");
this->size_ = this->allocated_size = size;
}
this->size_of_type = sizeof(T);
this->nb_component = nb_component;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::reserve(UInt new_size) {
UInt tmp_size = this->size_;
resizeUnitialized(new_size, false);
this->size_ = tmp_size;
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed. If the
* size increases, the new tuples are filled with zeros
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal> void Array<T, is_scal>::resize(UInt new_size) {
resizeUnitialized(new_size, !is_scal);
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed. If the
* size increases, the new tuples are filled with zeros
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal>
void Array<T, is_scal>::resize(UInt new_size, const T & val) {
this->resizeUnitialized(new_size, true, val);
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed.
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal>
void Array<T, is_scal>::resizeUnitialized(UInt new_size, bool fill,
const T & val) {
// AKANTU_DEBUG_IN();
// free some memory
if (new_size <= allocated_size) {
if (!is_scal) {
T * old_values = values;
if (new_size < size_) {
for (UInt i = new_size * nb_component; i < size_ * nb_component; ++i) {
T * obj = old_values + i;
obj->~T();
}
}
}
if (allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
AKANTU_DEBUG(dblAccessory,
"Freeing " << printMemorySize<T>((allocated_size - size_) *
nb_component)
<< " (" << id << ")");
// Normally there are no allocation problem when reducing an array
if (new_size == 0) {
free(values);
values = nullptr;
} else {
auto * tmp_ptr = static_cast<T *>(
realloc(values, new_size * nb_component * sizeof(T)));
if (tmp_ptr == nullptr) {
AKANTU_EXCEPTION("Cannot free data ("
<< id << ")"
<< " [current allocated size : " << allocated_size
<< " | "
<< "requested size : " << new_size << "]");
}
values = tmp_ptr;
}
allocated_size = new_size;
}
} else {
// allocate more memory
UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION)
? allocated_size + AKANTU_MIN_ALLOCATION
: new_size;
auto * tmp_ptr = static_cast<T *>(
realloc(values, size_to_alloc * nb_component * sizeof(T)));
AKANTU_DEBUG_ASSERT(
tmp_ptr != nullptr,
"Cannot allocate " << printMemorySize<T>(size_to_alloc * nb_component));
if (tmp_ptr == nullptr) {
AKANTU_DEBUG_ERROR("Cannot allocate more data ("
<< id << ")"
<< " [current allocated size : " << allocated_size
<< " | "
<< "requested size : " << new_size << "]");
}
AKANTU_DEBUG(dblAccessory,
"Allocating " << printMemorySize<T>(
(size_to_alloc - allocated_size) * nb_component));
allocated_size = size_to_alloc;
values = tmp_ptr;
}
if (fill && this->size_ < new_size) {
std::uninitialized_fill(values + size_ * nb_component,
values + new_size * nb_component, val);
}
size_ = new_size;
// AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* change the number of components by interlacing data
* @param multiplicator number of interlaced components add
* @param block_size blocks of data in the array
* Examaple for block_size = 2, multiplicator = 2
* array = oo oo oo -> new array = oo nn nn oo nn nn oo nn nn */
template <class T, bool is_scal>
void Array<T, is_scal>::extendComponentsInterlaced(UInt multiplicator,
UInt block_size) {
AKANTU_DEBUG_IN();
if (multiplicator == 1)
return;
AKANTU_DEBUG_ASSERT(multiplicator > 1, "invalid multiplicator");
AKANTU_DEBUG_ASSERT(nb_component % block_size == 0,
"stride must divide actual number of components");
values = static_cast<T *>(
realloc(values, nb_component * multiplicator * size_ * sizeof(T)));
UInt new_component = nb_component / block_size * multiplicator;
for (UInt i = 0, k = size_ - 1; i < size_; ++i, --k) {
for (UInt j = 0; j < new_component; ++j) {
UInt m = new_component - j - 1;
UInt n = m / multiplicator;
for (UInt l = 0, p = block_size - 1; l < block_size; ++l, --p) {
values[k * nb_component * multiplicator + m * block_size + p] =
values[k * nb_component + n * block_size + p];
}
}
}
nb_component = nb_component * multiplicator;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* search elem in the array, return the position of the first occurrence or
* -1 if not found
* @param elem the element to look for
* @return index of the first occurrence of elem or -1 if elem is not present
*/
template <class T, bool is_scal>
UInt Array<T, is_scal>::find(const T & elem) const {
AKANTU_DEBUG_IN();
auto begin = this->begin();
auto end = this->end();
auto it = std::find(begin, end, elem);
AKANTU_DEBUG_OUT();
return (it != end) ? it - begin : UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> UInt Array<T, is_scal>::find(T elem[]) const {
AKANTU_DEBUG_IN();
T * it = values;
UInt i = 0;
for (; i < size_; ++i) {
if (*it == elem[0]) {
T * cit = it;
UInt c = 0;
for (; (c < nb_component) && (*cit == elem[c]); ++c, ++cit)
;
if (c == nb_component) {
AKANTU_DEBUG_OUT();
return i;
}
}
it += nb_component;
}
return UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline UInt Array<T, is_scal>::find(const C<T> & elem) {
AKANTU_DEBUG_ASSERT(elem.size() == nb_component,
"Cannot find an element with a wrong size ("
<< elem.size() << ") != " << nb_component);
return this->find(elem.storage());
}
/* -------------------------------------------------------------------------- */
/**
* copy the content of another array. This overwrites the current content.
* @param other Array to copy into this array. It has to have the same
* nb_component as this. If compiled in debug mode, an incorrect other will
* result in an exception being thrown. Optimised code may result in
* unpredicted behaviour.
*/
template <class T, bool is_scal>
void Array<T, is_scal>::copy(const Array<T, is_scal> & vect,
bool no_sanity_check) {
AKANTU_DEBUG_IN();
if (!no_sanity_check)
if (vect.nb_component != nb_component)
AKANTU_DEBUG_ERROR(
"The two arrays do not have the same number of components");
resize((vect.size_ * vect.nb_component) / nb_component);
T * tmp = values;
std::uninitialized_copy(vect.storage(), vect.storage() + size_ * nb_component,
tmp);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <bool is_scal> class ArrayPrintHelper {
public:
template <typename T>
static void print_content(const Array<T> & vect, std::ostream & stream,
int indent) {
if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << " + values : {";
for (UInt i = 0; i < vect.size(); ++i) {
stream << "{";
for (UInt j = 0; j < vect.getNbComponent(); ++j) {
stream << vect(i, j);
if (j != vect.getNbComponent() - 1)
stream << ", ";
}
stream << "}";
if (i != vect.size() - 1)
stream << ", ";
}
stream << "}" << std::endl;
}
}
};
template <> class ArrayPrintHelper<false> {
public:
template <typename T>
static void print_content(__attribute__((unused)) const Array<T> & vect,
__attribute__((unused)) std::ostream & stream,
__attribute__((unused)) int indent) {}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
std::streamsize prec = stream.precision();
std::ios_base::fmtflags ff = stream.flags();
stream.setf(std::ios_base::showbase);
stream.precision(2);
stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> ["
<< std::endl;
stream << space << " + id : " << this->id << std::endl;
stream << space << " + size : " << this->size_ << std::endl;
stream << space << " + nb_component : " << this->nb_component << std::endl;
stream << space << " + allocated size : " << this->allocated_size
<< std::endl;
stream << space << " + memory size : "
<< printMemorySize<T>(allocated_size * nb_component) << std::endl;
if (!AKANTU_DEBUG_LEVEL_IS_TEST())
stream << space << " + address : " << std::hex << this->values
<< std::dec << std::endl;
stream.precision(prec);
stream.flags(ff);
ArrayPrintHelper<is_scal>::print_content(*this, stream, indent);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
/* Inline Functions ArrayBase */
/* -------------------------------------------------------------------------- */
inline UInt ArrayBase::getMemorySize() const {
return allocated_size * nb_component * size_of_type;
}
inline void ArrayBase::empty() { size_ = 0; }
#ifndef SWIG
/* -------------------------------------------------------------------------- */
/* Iterators */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <class R, class daughter, class IR, bool is_tensor>
class Array<T, is_scal>::iterator_internal {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using proxy = typename R::proxy;
using const_proxy = const typename R::proxy;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
iterator_internal() = default;
iterator_internal(pointer_type data, UInt _offset)
: _offset(_offset), initial(data), ret(nullptr), ret_ptr(data) {
AKANTU_DEBUG_ERROR(
"The constructor should never be called it is just an ugly trick...");
}
iterator_internal(std::unique_ptr<internal_value_type> && wrapped)
: _offset(wrapped->size()), initial(wrapped->storage()),
ret(std::move(wrapped)), ret_ptr(ret->storage()) {}
iterator_internal(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
}
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
if (this->ret)
this->ret->shallowCopy(*it.ret);
else
this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
return *this;
}
UInt getCurrentIndex() {
return (this->ret_ptr - this->initial) / this->_offset;
};
inline reference operator*() {
ret->values = ret_ptr;
return *ret;
};
inline const_reference operator*() const {
ret->values = ret_ptr;
return *ret;
};
inline pointer operator->() {
ret->values = ret_ptr;
return ret.get();
};
inline daughter & operator++() {
ret_ptr += _offset;
return static_cast<daughter &>(*this);
};
inline daughter & operator--() {
ret_ptr -= _offset;
return static_cast<daughter &>(*this);
};
inline daughter & operator+=(const UInt n) {
ret_ptr += _offset * n;
return static_cast<daughter &>(*this);
}
inline daughter & operator-=(const UInt n) {
ret_ptr -= _offset * n;
return static_cast<daughter &>(*this);
}
inline proxy operator[](const UInt n) {
ret->values = ret_ptr + n * _offset;
return proxy(*ret);
}
inline const_proxy operator[](const UInt n) const {
ret->values = ret_ptr + n * _offset;
return const_proxy(*ret);
}
inline bool operator==(const iterator_internal & other) const {
return this->ret_ptr == other.ret_ptr;
}
inline bool operator!=(const iterator_internal & other) const {
return this->ret_ptr != other.ret_ptr;
}
inline bool operator<(const iterator_internal & other) const {
return this->ret_ptr < other.ret_ptr;
}
inline bool operator<=(const iterator_internal & other) const {
return this->ret_ptr <= other.ret_ptr;
}
inline bool operator>(const iterator_internal & other) const {
return this->ret_ptr > other.ret_ptr;
}
inline bool operator>=(const iterator_internal & other) const {
return this->ret_ptr >= other.ret_ptr;
}
inline daughter operator+(difference_type n) {
daughter tmp(static_cast<daughter &>(*this));
tmp += n;
return tmp;
}
inline daughter operator-(difference_type n) {
daughter tmp(static_cast<daughter &>(*this));
tmp -= n;
return tmp;
}
inline difference_type operator-(const iterator_internal & b) {
return (this->ret_ptr - b.ret_ptr) / _offset;
}
inline pointer_type data() const { return ret_ptr; }
inline difference_type offset() const { return _offset; }
protected:
UInt _offset{0};
pointer_type initial{nullptr};
std::unique_ptr<internal_value_type> ret{nullptr};
pointer_type ret_ptr{nullptr};
};
/* -------------------------------------------------------------------------- */
/**
* Specialization for scalar types
*/
template <class T, bool is_scal>
template <class R, class daughter, class IR>
class Array<T, is_scal>::iterator_internal<R, daughter, IR, false> {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
iterator_internal(pointer data = nullptr) : ret(data), initial(data){};
iterator_internal(const iterator_internal & it) = default;
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) = default;
UInt getCurrentIndex() { return (this->ret - this->initial); };
inline reference operator*() { return *ret; };
inline const_reference operator*() const { return *ret; };
inline pointer operator->() { return ret; };
inline daughter & operator++() {
++ret;
return static_cast<daughter &>(*this);
};
inline daughter & operator--() {
--ret;
return static_cast<daughter &>(*this);
};
inline daughter & operator+=(const UInt n) {
ret += n;
return static_cast<daughter &>(*this);
}
inline daughter & operator-=(const UInt n) {
ret -= n;
return static_cast<daughter &>(*this);
}
inline reference operator[](const UInt n) { return ret[n]; }
inline bool operator==(const iterator_internal & other) const {
return ret == other.ret;
}
inline bool operator!=(const iterator_internal & other) const {
return ret != other.ret;
}
inline bool operator<(const iterator_internal & other) const {
return ret < other.ret;
}
inline bool operator<=(const iterator_internal & other) const {
return ret <= other.ret;
}
inline bool operator>(const iterator_internal & other) const {
return ret > other.ret;
}
inline bool operator>=(const iterator_internal & other) const {
return ret >= other.ret;
}
inline daughter operator-(difference_type n) { return daughter(ret - n); }
inline daughter operator+(difference_type n) { return daughter(ret + n); }
inline difference_type operator-(const iterator_internal & b) {
return ret - b.ret;
}
inline pointer data() const { return ret; }
protected:
pointer ret{nullptr};
pointer initial{nullptr};
};
/* -------------------------------------------------------------------------- */
/* Begin/End functions implementation */
/* -------------------------------------------------------------------------- */
namespace detail {
template <class Tuple, size_t... Is>
constexpr auto take_front_impl(Tuple && t, std::index_sequence<Is...>) {
return std::make_tuple(std::get<Is>(std::forward<Tuple>(t))...);
}
template <size_t N, class Tuple> constexpr auto take_front(Tuple && t) {
return take_front_impl(std::forward<Tuple>(t),
std::make_index_sequence<N>{});
}
template <typename... V>
constexpr auto product_all(V &&... v) ->
typename std::common_type<V...>::type {
typename std::common_type<V...>::type result = 1;
(void)std::initializer_list<int>{(result *= v, 0)...};
return result;
}
template <typename... T> std::string to_string_all(T &&... t) {
if (sizeof...(T) == 0)
return "";
std::stringstream ss;
bool noComma = true;
ss << "(";
(void)std::initializer_list<bool>{
(ss << (noComma ? "" : ", ") << t, noComma = false)...};
ss << ")";
return ss.str();
}
template <std::size_t N> struct InstantiationHelper {
template <typename type, typename T, typename... Ns>
static auto instantiate(T && data, Ns... ns) {
return std::make_unique<type>(data, ns...);
}
};
template <> struct InstantiationHelper<0> {
template <typename type, typename T> static auto instantiate(T && data) {
return data;
}
};
template <typename Arr, typename T, typename... Ns>
decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) {
using type = IteratorHelper_t<sizeof...(Ns) - 1, T>;
using array_type = std::decay_t<Arr>;
using iterator =
- std::conditional_t<std::is_const<Arr>::value,
+ std::conditional_t<std::is_const<std::remove_reference_t<Arr>>::value,
typename array_type::template const_iterator<type>,
typename array_type::template iterator<type>>;
+ static_assert(sizeof...(Ns), "You should provide a least one size");
+
if (array.getNbComponent() * array.size() !=
product_all(std::forward<Ns>(ns)...)) {
AKANTU_CUSTOM_EXCEPTION_INFO(
debug::ArrayException(),
- "The iterator on Array "
+ "The iterator on " << debug::demangle(typeid(Arr).name())
<< to_string_all(array.size(), array.getNbComponent())
<< "is not compatible with the type "
<< debug::demangle(typeid(type).name()) << to_string_all(ns...));
}
auto && wrapped = aka::apply(
[&](auto... n) {
return InstantiationHelper<sizeof...(n)>::template instantiate<type>(
data, n...);
},
take_front<sizeof...(Ns) - 1>(std::make_tuple(ns...)));
return iterator(std::move(wrapped));
}
} // namespace detail
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) {
return detail::get_iterator(*this, values + nb_component * size_,
std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) const {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) const {
return detail::get_iterator(*this, values + nb_component * size_,
std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) {
return detail::get_iterator(
*this, values + detail::product_all(std::forward<Ns>(ns)...),
std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) const {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) const {
return detail::get_iterator(
*this, values + detail::product_all(std::forward<Ns>(ns)...),
std::forward<Ns>(ns)...);
}
/* -------------------------------------------------------------------------- */
/* Views */
/* -------------------------------------------------------------------------- */
namespace detail {
template <typename Array, typename... Ns> class ArrayView {
using tuple = std::tuple<Ns...>;
public:
ArrayView(Array && array, Ns... ns)
: array(std::forward<Array>(array)), sizes(std::move(ns)...){};
decltype(auto) begin() {
return aka::apply(
[&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
}
decltype(auto) end() {
return aka::apply(
[&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
}
decltype(auto) size() {
return std::get<std::tuple_size<tuple>::value - 1>(sizes);
}
decltype(auto) dims() {
return std::tuple_size<tuple>::value - 1;
}
private:
Array array;
tuple sizes;
};
} // namespace detail
/* -------------------------------------------------------------------------- */
template <typename Array, typename... Ns>
decltype(auto) make_view(Array && array, Ns... ns) {
+ static_assert(sizeof...(Ns), "You should provide a least one dimension");
auto size = std::forward<decltype(array)>(array).size() *
std::forward<decltype(array)>(array).getNbComponent() /
detail::product_all(ns...);
return detail::ArrayView<Array, std::common_type_t<size_t, Ns>...,
std::common_type_t<size_t, decltype(size)>>(
std::forward<Array>(array), std::move(ns)..., size);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::const_iterator
: public iterator_internal<const R, Array<T, is_scal>::const_iterator<R>,
R> {
public:
using parent = iterator_internal<const R, const_iterator, R>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
const_iterator() : parent(){};
// const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}
// const_iterator(pointer warped) : parent(warped) {}
// const_iterator(const parent & it) : parent(it) {}
const_iterator(const const_iterator & it) = default;
const_iterator(const_iterator && it) = default;
template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
const_iterator(P * data) : parent(data) {}
template <typename UP_P,
typename =
std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
const_iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
const_iterator & operator=(const const_iterator & it) = default;
};
template <class T, class R, bool is_tensor_ = is_tensor<R>::value>
struct ConstConverterIteratorHelper {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
return const_iterator(std::unique_ptr<R>(new R(*it, false)));
}
};
template <class T, class R> struct ConstConverterIteratorHelper<T, R, false> {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
return const_iterator(it.data());
}
};
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::iterator
: public iterator_internal<R, Array<T, is_scal>::iterator<R>> {
public:
using parent = iterator_internal<R, iterator>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
iterator() : parent(){};
iterator(const iterator & it) = default;
iterator(iterator && it) = default;
template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
iterator(P * data) : parent(data) {}
template <typename UP_P,
typename =
std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
iterator & operator=(const iterator & it) = default;
operator const_iterator<R>() {
return ConstConverterIteratorHelper<T, R>::convert(*this);
}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
inline typename Array<T, is_scal>::template iterator<R>
Array<T, is_scal>::erase(const iterator<R> & it) {
T * curr = it.data();
UInt pos = (curr - values) / nb_component;
erase(pos);
iterator<R> rit = it;
return --rit;
}
#endif
} // namespace akantu
#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */
diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in
index 8f7e5057c..ecbec6274 100644
--- a/src/common/aka_config.hh.in
+++ b/src/common/aka_config.hh.in
@@ -1,103 +1,103 @@
/**
* @file aka_config.hh.in
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Sun Sep 26 2010
* @date last modification: Thu Jan 21 2016
*
* @brief Compilation time configuration of Akantu
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_CONFIG_HH__
#define __AKANTU_AKA_CONFIG_HH__
#define AKANTU_VERSION_MAJOR @AKANTU_MAJOR_VERSION@
#define AKANTU_VERSION_MINOR @AKANTU_MINOR_VERSION@
#define AKANTU_VERSION_PATCH @AKANTU_PATCH_VERSION@
-#define AKANTU_VERSION (AKANTU_VERSION_MAJOR * 100000 \
- + AKANTU_VERSION_MINOR * 1000 \
+#define AKANTU_VERSION (AKANTU_VERSION_MAJOR * 10000 \
+ + AKANTU_VERSION_MINOR * 100 \
+ AKANTU_VERSION_PATCH)
@AKANTU_TYPES_EXTRA_INCLUDES@
namespace akantu {
using Real = @AKANTU_FLOAT_TYPE@;
using Int = @AKANTU_SIGNED_INTEGER_TYPE@;
using UInt = @AKANTU_UNSIGNED_INTEGER_TYPE@;
} // akantu
#define AKANTU_INTEGER_SIZE @AKANTU_INTEGER_SIZE@
#define AKANTU_FLOAT_SIZE @AKANTU_FLOAT_SIZE@
#cmakedefine AKANTU_HAS_STRDUP
#cmakedefine AKANTU_USE_BLAS
#cmakedefine AKANTU_USE_LAPACK
#cmakedefine AKANTU_PARALLEL
#cmakedefine AKANTU_USE_MPI
#cmakedefine AKANTU_USE_SCOTCH
#cmakedefine AKANTU_USE_PTSCOTCH
#cmakedefine AKANTU_SCOTCH_NO_EXTERN
#cmakedefine AKANTU_IMPLICIT
#cmakedefine AKANTU_USE_MUMPS
#cmakedefine AKANTU_USE_PETSC
#cmakedefine AKANTU_USE_IOHELPER
#cmakedefine AKANTU_USE_QVIEW
#cmakedefine AKANTU_USE_BLACKDYNAMITE
#cmakedefine AKANTU_USE_PYBIND11
#cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY
#cmakedefine AKANTU_EXTRA_MATERIALS
#cmakedefine AKANTU_STUDENTS_EXTRA_PACKAGE
#cmakedefine AKANTU_DAMAGE_NON_LOCAL
#cmakedefine AKANTU_SOLID_MECHANICS
#cmakedefine AKANTU_STRUCTURAL_MECHANICS
#cmakedefine AKANTU_HEAT_TRANSFER
#cmakedefine AKANTU_PYTHON_INTERFACE
#cmakedefine AKANTU_COHESIVE_ELEMENT
#cmakedefine AKANTU_PARALLEL_COHESIVE_ELEMENT
#cmakedefine AKANTU_IGFEM
#cmakedefine AKANTU_USE_CGAL
#cmakedefine AKANTU_EMBEDDED
// Debug tools
//#cmakedefine AKANTU_NDEBUG
#cmakedefine AKANTU_DEBUG_TOOLS
#cmakedefine READLINK_COMMAND @READLINK_COMMAND@
#cmakedefine ADDR2LINE_COMMAND @ADDR2LINE_COMMAND@
#define __aka_inline__ inline
#endif /* __AKANTU_AKA_CONFIG_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index fcc686edf..3693383f0 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1194 +1,1195 @@
/**
* @file solid_mechanics_model.cc
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Clement Roux <clement.roux@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Tue Jan 19 2016
*
* @brief Implementation of the SolidMechanicsModel 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "solid_mechanics_model_tmpl.hh"
#include "communicator.hh"
#include "element_synchronizer.hh"
#include "sparse_matrix.hh"
#include "synchronizer_registry.hh"
#include "dumpable_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
#include "material_non_local.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
*/
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id,
const MemoryID & memory_id, const ModelType model_type)
: Model(mesh, model_type, dim, id, memory_id),
BoundaryCondition<SolidMechanicsModel>(), f_m2a(1.0),
displacement(nullptr), previous_displacement(nullptr),
displacement_increment(nullptr), mass(nullptr), velocity(nullptr),
acceleration(nullptr), external_force(nullptr), internal_force(nullptr),
blocked_dofs(nullptr), current_position(nullptr),
material_index("material index", id, memory_id),
material_local_numbering("material local numbering", id, memory_id),
increment_flag(false),
are_materials_instantiated(false) { //, pbc_synch(nullptr) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
Model::spatial_dimension);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
_ek_regular);
#endif
material_selector = std::make_shared<DefaultMaterialSelector>(material_index),
this->initDOFManager();
this->registerDataAccessor(*this);
if (this->mesh.isDistributed()) {
auto & synchronizer = this->mesh.getElementSynchronizer();
this->registerSynchronizer(synchronizer, _gst_material_id);
this->registerSynchronizer(synchronizer, _gst_smm_mass);
this->registerSynchronizer(synchronizer, _gst_smm_stress);
this->registerSynchronizer(synchronizer, _gst_smm_boundary);
this->registerSynchronizer(synchronizer, _gst_for_dump);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
for (auto & internal : this->registered_internals) {
delete internal.second;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) {
Model::setTimeStep(time_step, solver_id);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialization */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilities
*/
void SolidMechanicsModel::initFullImpl(const ModelOptions & options) {
material_index.initialize(mesh, _element_kind = _ek_not_defined,
_default_value = UInt(-1), _with_nb_element = true);
material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true);
Model::initFullImpl(options);
// initialize pbc
if (this->pbc_pair.size() != 0)
this->initPBC();
// initialize the materials
if (this->parser.getLastParsedFile() != "") {
this->instantiateMaterials();
}
this->initMaterials();
this->initBC(*this, *displacement, *displacement_increment, *external_force);
}
/* -------------------------------------------------------------------------- */
TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const {
return _tsst_dynamic_lumped;
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case _tsst_dynamic_lumped: {
options.non_linear_solver_type = _nls_lumped;
options.integration_scheme_type["displacement"] = _ist_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case _tsst_static: {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
case _tsst_dynamic: {
if (this->method == _explicit_consistent_mass) {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
} else {
options.non_linear_solver_type = _nls_newton_raphson;
options.integration_scheme_type["displacement"] = _ist_trapezoidal_rule_2;
options.solution_type["displacement"] = IntegrationScheme::_displacement;
}
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
}
return options;
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped", _tsst_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", _tsst_dynamic);
}
case _static: {
return std::make_tuple("static", _tsst_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", _tsst_dynamic);
}
default:
return std::make_tuple("unknown", _tsst_not_defined);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type,
NonLinearSolverType) {
auto & dof_manager = this->getDOFManager();
/* ------------------------------------------------------------------------ */
// for alloc type of solvers
this->allocNodalField(this->displacement, spatial_dimension, "displacement");
this->allocNodalField(this->previous_displacement, spatial_dimension,
"previous_displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
"displacement_increment");
this->allocNodalField(this->internal_force, spatial_dimension,
"internal_force");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs");
this->allocNodalField(this->current_position, spatial_dimension,
"current_position");
// initialize the current positions
this->current_position->copy(this->mesh.getNodes());
/* ------------------------------------------------------------------------ */
if (!dof_manager.hasDOFs("displacement")) {
dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
dof_manager.registerDOFsIncrement("displacement",
*this->displacement_increment);
dof_manager.registerDOFsPrevious("displacement",
*this->previous_displacement);
}
/* ------------------------------------------------------------------------ */
// for dynamic
if (time_step_solver_type == _tsst_dynamic ||
time_step_solver_type == _tsst_dynamic_lumped) {
this->allocNodalField(this->velocity, spatial_dimension, "velocity");
this->allocNodalField(this->acceleration, spatial_dimension,
"acceleration");
if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
dof_manager.registerDOFsDerivative("displacement", 2,
*this->acceleration);
}
}
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
/* ------------------------------------------------------------------------ */
// computes the internal forces
this->assembleInternalForces();
/* ------------------------------------------------------------------------ */
this->getDOFManager().assembleToResidual("displacement",
*this->external_force, 1);
this->getDOFManager().assembleToResidual("displacement",
*this->internal_force, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) {
// \TODO check the materials to know what is the correct answer
if (matrix_id == "C")
return _mt_not_defined;
return _symmetric;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
} else if (matrix_id == "M") {
this->assembleMass();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) {
if (matrix_id == "M") {
this->assembleMassLumped();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::beforeSolveStep() {
for (auto & material : materials)
material->beforeSolveStep();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::afterSolveStep() {
for (auto & material : materials)
material->afterSolveStep();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::predictor() { ++displacement_release; }
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::corrector() { ++displacement_release; }
/* -------------------------------------------------------------------------- */
/**
* This function computes the internal forces as F_{int} = \int_{\Omega} N
* \sigma d\Omega@f$
*/
void SolidMechanicsModel::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
this->internal_force->clear();
// compute the stresses of local elements
AKANTU_DEBUG_INFO("Compute local stresses");
for (auto & material : materials) {
material->computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
if (this->non_local_manager)
this->non_local_manager->computeAllNonLocalStresses();
// communicate the stresses
AKANTU_DEBUG_INFO("Send data for residual assembly");
this->asynchronousSynchronize(_gst_smm_stress);
// assemble the forces due to local stresses
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for (auto & material : materials) {
material->assembleInternalForces(_not_ghost);
}
// finalize communications
AKANTU_DEBUG_INFO("Wait distant stresses");
this->waitEndSynchronize(_gst_smm_stress);
// assemble the stresses due to ghost elements
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for (auto & material : materials) {
material->assembleInternalForces(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
// Check if materials need to recompute the matrix
bool need_to_reassemble = false;
for (auto & material : materials) {
need_to_reassemble |= material->hasStiffnessMatrixChanged();
}
if (need_to_reassemble) {
this->getDOFManager().getMatrix("K").clear();
// call compute stiffness matrix on each local elements
for (auto & material : materials) {
material->assembleStiffnessMatrix(_not_ghost);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
if (this->current_position_release == this->displacement_release)
return;
this->current_position->copy(this->mesh.getNodes());
auto cpos_it = this->current_position->begin(Model::spatial_dimension);
auto cpos_end = this->current_position->end(Model::spatial_dimension);
auto disp_it = this->displacement->begin(Model::spatial_dimension);
for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
*cpos_it += *disp_it;
}
this->current_position_release = this->displacement_release;
}
/* -------------------------------------------------------------------------- */
const Array<Real> & SolidMechanicsModel::getCurrentPosition() {
this->updateCurrentPosition();
return *this->current_position;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateDataForNonLocalCriterion(
ElementTypeMapReal & criterion) {
const ID field_name = criterion.getName();
for (auto & material : materials) {
if (!material->isInternal<Real>(field_name, _ek_regular))
continue;
for (auto ghost_type : ghost_types) {
material->flattenInternal(field_name, criterion, ghost_type, _ek_regular);
}
}
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Real min_dt = std::numeric_limits<Real>::max();
this->updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
for (auto type :
mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) {
elem.type = type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
UInt nb_element = mesh.getNbElement(type);
auto mat_indexes = material_index(type, ghost_type).begin();
auto mat_loc_num = material_local_numbering(type, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position, X, type,
_not_ghost);
auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element;
++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
elem.element = *mat_loc_num;
Real el_h = getFEEngine().getElementInradius(*X_el, type);
Real el_c = this->materials[*mat_indexes]->getCelerity(elem);
Real el_dt = el_h / el_c;
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
if (this->getDOFManager().hasLumpedMatrix("M")) {
auto m_it = this->mass->begin(Model::spatial_dimension);
auto m_end = this->mass->end(Model::spatial_dimension);
auto v_it = this->velocity->begin(Model::spatial_dimension);
for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) {
const Vector<Real> & v = *v_it;
const Vector<Real> & m = *m_it;
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
// bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node; // && is_not_pbc_slave_node;
if (count_node) {
for (UInt i = 0; i < Model::spatial_dimension; ++i) {
if (m(i) > std::numeric_limits<Real>::epsilon())
mv2 += v(i) * v(i) * m(i);
}
}
ekin += mv2;
}
} else if (this->getDOFManager().hasMatrix("M")) {
Array<Real> Mv(nb_nodes, Model::spatial_dimension);
this->getDOFManager().getMatrix("M").matVecMul(*this->velocity, Mv);
auto mv_it = Mv.begin(Model::spatial_dimension);
auto mv_end = Mv.end(Model::spatial_dimension);
auto v_it = this->velocity->begin(Model::spatial_dimension);
for (; mv_it != mv_end; ++mv_it, ++v_it) {
ekin += v_it->dot(*mv_it);
}
} else {
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
}
mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(const ElementType & type,
UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
Model::spatial_dimension, type,
_not_ghost, filter_element);
Array<Real>::vector_iterator vit =
vel_on_quad.begin(Model::spatial_dimension);
Array<Real>::vector_iterator vend = vel_on_quad.end(Model::spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[material_index(type)(index)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5 * getFEEngine().integrate(rho_v2, type, index);
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
auto ext_force_it = external_force->begin(Model::spatial_dimension);
auto int_force_it = internal_force->begin(Model::spatial_dimension);
auto boun_it = blocked_dofs->begin(Model::spatial_dimension);
decltype(ext_force_it) incr_or_velo_it;
if (this->method == _static) {
incr_or_velo_it =
this->displacement_increment->begin(Model::spatial_dimension);
} else {
incr_or_velo_it = this->velocity->begin(Model::spatial_dimension);
}
Real work = 0.;
UInt nb_nodes = this->mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes;
++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) {
const auto & int_force = *int_force_it;
const auto & ext_force = *ext_force_it;
const auto & boun = *boun_it;
const auto & incr_or_velo = *incr_or_velo_it;
bool is_local_node = this->mesh.isLocalOrMasterNode(n);
// bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
bool count_node = is_local_node; // && is_not_pbc_slave_node;
if (count_node) {
for (UInt i = 0; i < Model::spatial_dimension; ++i) {
if (boun(i))
work -= int_force(i) * incr_or_velo(i);
else
work += ext_force(i) * incr_or_velo(i);
}
}
}
mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum);
if (this->method != _static)
work *= this->getTimeStep();
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy();
} else if (energy_id == "external work") {
return getExternalWork();
}
Real energy = 0.;
for (auto & material : materials)
energy += material->getEnergy(energy_id);
/// reduction sum over all processors
mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy(type, index);
}
UInt mat_index = this->material_index(type, _not_ghost)(index);
UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
Real energy =
this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
this->material_index.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true,
_default_value = UInt(-1));
this->material_local_numbering.initialize(
mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
_default_value = UInt(-1));
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID(),
this->getMemoryID());
for (auto & elem : element_list) {
if (!filter.exists(elem.type, elem.ghost_type))
filter.alloc(0, 1, elem.type, elem.ghost_type);
filter(elem.type, elem.ghost_type).push_back(elem.element);
}
this->assignMaterialToElements(&filter);
for (auto & material : materials)
material->onElementsAdded(element_list, event);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(
const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
for (auto & material : materials) {
material->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if (displacement) {
displacement->resize(nb_nodes, 0.);
++displacement_release;
}
if (mass)
mass->resize(nb_nodes, 0.);
if (velocity)
velocity->resize(nb_nodes, 0.);
if (acceleration)
acceleration->resize(nb_nodes, 0.);
if (external_force)
external_force->resize(nb_nodes, 0.);
if (internal_force)
internal_force->resize(nb_nodes, 0.);
if (blocked_dofs)
blocked_dofs->resize(nb_nodes, 0.);
if (current_position)
current_position->resize(nb_nodes, 0.);
if (previous_displacement)
previous_displacement->resize(nb_nodes, 0.);
if (displacement_increment)
displacement_increment->resize(nb_nodes, 0.);
for (auto & material : materials) {
material->onNodesAdded(nodes_list, event);
}
need_to_reassemble_lumped_mass = true;
need_to_reassemble_mass = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/,
const Array<UInt> & new_numbering,
const RemovedNodesEvent & /*event*/) {
if (displacement) {
mesh.removeNodesFromArray(*displacement, new_numbering);
++displacement_release;
}
if (mass)
mesh.removeNodesFromArray(*mass, new_numbering);
if (velocity)
mesh.removeNodesFromArray(*velocity, new_numbering);
if (acceleration)
mesh.removeNodesFromArray(*acceleration, new_numbering);
if (internal_force)
mesh.removeNodesFromArray(*internal_force, new_numbering);
if (external_force)
mesh.removeNodesFromArray(*external_force, new_numbering);
if (blocked_dofs)
mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
// if (increment_acceleration)
// mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
if (displacement_increment)
mesh.removeNodesFromArray(*displacement_increment, new_numbering);
if (previous_displacement)
mesh.removeNodesFromArray(*previous_displacement, new_numbering);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << Model::spatial_dimension
<< std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
if (velocity)
velocity->printself(stream, indent + 2);
if (acceleration)
acceleration->printself(stream, indent + 2);
if (mass)
mass->printself(stream, indent + 2);
external_force->printself(stream, indent + 2);
internal_force->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + materials [" << std::endl;
for (auto & material : materials) {
material->printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeNonLocal() {
this->non_local_manager->synchronize(*this, _gst_material_id);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods(
const GhostType & ghost_type) {
for (auto & mat : materials) {
MaterialNonLocalInterface * mat_non_local;
if ((mat_non_local =
dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr)
continue;
ElementTypeMapArray<Real> quadrature_points_coordinates(
"quadrature_points_coordinates_tmp_nl", this->id, this->memory_id);
quadrature_points_coordinates.initialize(this->getFEEngine(),
_nb_component = spatial_dimension,
_ghost_type = ghost_type);
for (auto & type : quadrature_points_coordinates.elementTypes(
Model::spatial_dimension, ghost_type)) {
this->getFEEngine().computeIntegrationPointsCoordinates(
quadrature_points_coordinates(type, ghost_type), type, ghost_type);
}
mat_non_local->initMaterialNonLocal();
mat_non_local->insertIntegrationPointsInNeighborhoods(
ghost_type, quadrature_points_coordinates);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeNonLocalStresses(
const GhostType & ghost_type) {
for (auto & mat : materials) {
if (dynamic_cast<MaterialNonLocalInterface *>(mat.get()) == nullptr)
continue;
auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
mat_non_local.computeNonLocalStresses(ghost_type);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateLocalInternal(
ElementTypeMapReal & internal_flat, const GhostType & ghost_type,
const ElementKind & kind) {
const ID field_name = internal_flat.getName();
for (auto & material : materials) {
if (material->isInternal<Real>(field_name, kind))
material->flattenInternal(field_name, internal_flat, ghost_type, kind);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateNonLocalInternal(
ElementTypeMapReal & internal_flat, const GhostType & ghost_type,
const ElementKind & kind) {
const ID field_name = internal_flat.getName();
for (auto & mat : materials) {
if (dynamic_cast<MaterialNonLocalInterface *>(mat.get()) == nullptr)
continue;
auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
kind);
}
}
/* -------------------------------------------------------------------------- */
FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) {
return dynamic_cast<FEEngine &>(
getFEEngineClassBoundary<MyFEEngineType>(name));
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::splitElementByMaterial(
const Array<Element> & elements,
std::vector<Array<Element>> & elements_per_mat) const {
ElementType current_element_type = _not_defined;
GhostType current_ghost_type = _casper;
const Array<UInt> * mat_indexes = nullptr;
const Array<UInt> * mat_loc_num = nullptr;
- for (auto && el : elements) {
+ for (const auto & el : elements) {
if (el.type != current_element_type ||
el.ghost_type != current_ghost_type) {
current_element_type = el.type;
current_ghost_type = el.ghost_type;
mat_indexes = &(this->material_index(el.type, el.ghost_type));
mat_loc_num = &(this->material_local_numbering(el.type, el.ghost_type));
}
- UInt old_id = el.element;
- el.element = (*mat_loc_num)(old_id);
- elements_per_mat[(*mat_indexes)(old_id)].push_back(el);
+ Element mat_el = el;
+
+ mat_el.element = (*mat_loc_num)(el.element);
+ elements_per_mat[(*mat_indexes)(el.element)].push_back(mat_el);
}
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModel::getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
UInt nb_nodes_per_element = 0;
for (const Element & el : elements) {
nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
}
switch (tag) {
case _gst_material_id: {
size += elements.size() * sizeof(UInt);
break;
}
case _gst_smm_mass: {
size += nb_nodes_per_element * sizeof(Real) *
Model::spatial_dimension; // mass vector
break;
}
case _gst_smm_for_gradu: {
size += nb_nodes_per_element * Model::spatial_dimension *
sizeof(Real); // displacement
break;
}
case _gst_smm_boundary: {
// force, displacement, boundary
size += nb_nodes_per_element * Model::spatial_dimension *
(2 * sizeof(Real) + sizeof(bool));
break;
}
case _gst_for_dump: {
// displacement, velocity, acceleration, residual, force
size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5;
break;
}
default: {}
}
if (tag != _gst_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
size += mat.getNbData(elements, tag);
});
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
switch (tag) {
case _gst_material_id: {
this->packElementalDataHelper(material_index, buffer, elements, false,
getFEEngine());
break;
}
case _gst_smm_mass: {
packNodalDataHelper(*mass, buffer, elements, mesh);
break;
}
case _gst_smm_for_gradu: {
packNodalDataHelper(*displacement, buffer, elements, mesh);
break;
}
case _gst_for_dump: {
packNodalDataHelper(*displacement, buffer, elements, mesh);
packNodalDataHelper(*velocity, buffer, elements, mesh);
packNodalDataHelper(*acceleration, buffer, elements, mesh);
packNodalDataHelper(*internal_force, buffer, elements, mesh);
packNodalDataHelper(*external_force, buffer, elements, mesh);
break;
}
case _gst_smm_boundary: {
packNodalDataHelper(*external_force, buffer, elements, mesh);
packNodalDataHelper(*velocity, buffer, elements, mesh);
packNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
break;
}
default: {}
}
if (tag != _gst_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
mat.packData(buffer, elements, tag);
});
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
switch (tag) {
case _gst_material_id: {
for (auto && element : elements) {
UInt recv_mat_index;
buffer >> recv_mat_index;
UInt & mat_index = material_index(element);
if (mat_index != UInt(-1))
continue;
// add ghosts element to the correct material
mat_index = recv_mat_index;
UInt index = materials[mat_index]->addElement(element);
material_local_numbering(element) = index;
}
break;
}
case _gst_smm_mass: {
unpackNodalDataHelper(*mass, buffer, elements, mesh);
break;
}
case _gst_smm_for_gradu: {
unpackNodalDataHelper(*displacement, buffer, elements, mesh);
break;
}
case _gst_for_dump: {
unpackNodalDataHelper(*displacement, buffer, elements, mesh);
unpackNodalDataHelper(*velocity, buffer, elements, mesh);
unpackNodalDataHelper(*acceleration, buffer, elements, mesh);
unpackNodalDataHelper(*internal_force, buffer, elements, mesh);
unpackNodalDataHelper(*external_force, buffer, elements, mesh);
break;
}
case _gst_smm_boundary: {
unpackNodalDataHelper(*external_force, buffer, elements, mesh);
unpackNodalDataHelper(*velocity, buffer, elements, mesh);
unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
break;
}
default: {}
}
if (tag != _gst_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
mat.unpackData(buffer, elements, tag);
});
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
// UInt nb_nodes = mesh.getNbNodes();
switch (tag) {
case _gst_smm_uv: {
size += sizeof(Real) * Model::spatial_dimension * 2;
break;
}
case _gst_smm_res: {
size += sizeof(Real) * Model::spatial_dimension;
break;
}
case _gst_smm_mass: {
size += sizeof(Real) * Model::spatial_dimension;
break;
}
case _gst_for_dump: {
size += sizeof(Real) * Model::spatial_dimension * 5;
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size * dofs.size();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const Array<UInt> & dofs,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
switch (tag) {
case _gst_smm_uv: {
packDOFDataHelper(*displacement, buffer, dofs);
packDOFDataHelper(*velocity, buffer, dofs);
break;
}
case _gst_smm_res: {
packDOFDataHelper(*internal_force, buffer, dofs);
break;
}
case _gst_smm_mass: {
packDOFDataHelper(*mass, buffer, dofs);
break;
}
case _gst_for_dump: {
packDOFDataHelper(*displacement, buffer, dofs);
packDOFDataHelper(*velocity, buffer, dofs);
packDOFDataHelper(*acceleration, buffer, dofs);
packDOFDataHelper(*internal_force, buffer, dofs);
packDOFDataHelper(*external_force, buffer, dofs);
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const Array<UInt> & dofs,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
switch (tag) {
case _gst_smm_uv: {
unpackDOFDataHelper(*displacement, buffer, dofs);
unpackDOFDataHelper(*velocity, buffer, dofs);
break;
}
case _gst_smm_res: {
unpackDOFDataHelper(*internal_force, buffer, dofs);
break;
}
case _gst_smm_mass: {
unpackDOFDataHelper(*mass, buffer, dofs);
break;
}
case _gst_for_dump: {
unpackDOFDataHelper(*displacement, buffer, dofs);
unpackDOFDataHelper(*velocity, buffer, dofs);
unpackDOFDataHelper(*acceleration, buffer, dofs);
unpackDOFDataHelper(*internal_force, buffer, dofs);
unpackDOFDataHelper(*external_force, buffer, dofs);
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
} // namespace akantu
diff --git a/src/synchronizer/communication_descriptor.hh b/src/synchronizer/communication_descriptor.hh
index e3f24ad99..609218b86 100644
--- a/src/synchronizer/communication_descriptor.hh
+++ b/src/synchronizer/communication_descriptor.hh
@@ -1,153 +1,153 @@
/**
* @file synchronizer_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Aug 3 13:49:36 2016
*
* @brief Implementation of the helper classes for the synchronizer
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "communication_tag.hh"
#include "data_accessor.hh"
#include "communication_request.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMUNICATION_DESCRIPTOR_HH__
#define __AKANTU_COMMUNICATION_DESCRIPTOR_HH__
namespace akantu {
/* ------------------------------------------------------------------------ */
enum CommunicationSendRecv { _send, _recv, _csr_not_defined };
/* -------------------------------------------------------------------------- */
struct CommunicationSRType {
using type = CommunicationSendRecv;
static const type _begin_ = _send;
static const type _end_ = _csr_not_defined;
};
using send_recv_t = safe_enum<CommunicationSRType>;
namespace {
send_recv_t iterate_send_recv{};
}
/* ------------------------------------------------------------------------ */
class Communication {
public:
explicit Communication(const CommunicationSendRecv & type = _csr_not_defined)
: _type(type) {}
Communication(const Communication &) = delete;
Communication & operator=(const Communication &) = delete;
void resize(UInt size) {
this->_size = size;
this->_buffer.resize(size);
}
inline const CommunicationSendRecv & type() const { return this->_type; }
inline const UInt & size() const { return this->_size; }
inline const CommunicationRequest & request() const { return this->_request; }
inline CommunicationRequest & request() { return this->_request; }
inline const CommunicationBuffer & buffer() const { return this->_buffer; }
inline CommunicationBuffer & buffer() { return this->_buffer; }
private:
UInt _size{0};
CommunicationBuffer _buffer;
CommunicationRequest _request;
CommunicationSendRecv _type;
};
template <class Entity> class Communications;
/* ------------------------------------------------------------------------ */
template <class Entity> class CommunicationDescriptor {
public:
CommunicationDescriptor(Communication & communication, Array<Entity> & scheme,
Communications<Entity> & communications,
const SynchronizationTag & tag, UInt proc);
CommunicationDescriptor(const CommunicationDescriptor &) = default;
CommunicationDescriptor &
operator=(const CommunicationDescriptor &) = default;
/// get the quantity of data in the buffer
UInt getNbData() { return communication.size(); }
/// set the quantity of data in the buffer
void setNbData(UInt size) { communication.resize(size); }
/// get the corresponding tag
const SynchronizationTag & getTag() const { return tag; }
/// get the data buffer
CommunicationBuffer & getBuffer();
/// get the corresponding request
CommunicationRequest & getRequest();
/// get the communication scheme
- Array<Entity> & getScheme();
+ const Array<Entity> & getScheme();
/// reset the buffer before pack or after unpack
void resetBuffer();
/// pack data for entities in the buffer
void packData(const DataAccessor<Entity> & accessor);
/// unpack data for entities from the buffer
void unpackData(DataAccessor<Entity> & accessor);
/// posts asynchronous send requests
void postSend(int hash_id);
/// posts asynchronous receive requests
void postRecv(int hash_id);
/// free the request
void freeRequest();
UInt getProc() { return proc; }
protected:
Communication & communication;
- Array<Entity> & scheme;
+ const Array<Entity> & scheme;
Communications<Entity> & communications;
const SynchronizationTag & tag;
UInt proc;
UInt rank;
UInt counter;
};
/* -------------------------------------------------------------------------- */
} // namespace akantu
#include "communication_descriptor_tmpl.hh"
#endif /* __AKANTU_COMMUNICATION_DESCRIPTOR_HH__ */
diff --git a/src/synchronizer/communication_descriptor_tmpl.hh b/src/synchronizer/communication_descriptor_tmpl.hh
index 45c57628f..8e6982fec 100644
--- a/src/synchronizer/communication_descriptor_tmpl.hh
+++ b/src/synchronizer/communication_descriptor_tmpl.hh
@@ -1,152 +1,152 @@
/**
* @file communication_descriptor_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Sep 6 14:03:16 2016
*
* @brief implementation of CommunicationDescriptor
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "communication_descriptor.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMUNICATION_DESCRIPTOR_TMPL_HH__
#define __AKANTU_COMMUNICATION_DESCRIPTOR_TMPL_HH__
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Implementations */
/* -------------------------------------------------------------------------- */
template <class Entity>
CommunicationDescriptor<Entity>::CommunicationDescriptor(
Communication & communication, Array<Entity> & scheme,
Communications<Entity> & communications, const SynchronizationTag & tag,
UInt proc)
: communication(communication), scheme(scheme),
communications(communications), tag(tag), proc(proc),
rank(communications.getCommunicator().whoAmI()) {
counter = communications.getCounter(tag);
}
/* -------------------------------------------------------------------------- */
template <class Entity>
CommunicationBuffer & CommunicationDescriptor<Entity>::getBuffer() {
return communication.buffer();
}
/* -------------------------------------------------------------------------- */
template <class Entity>
CommunicationRequest & CommunicationDescriptor<Entity>::getRequest() {
return communication.request();
}
/* -------------------------------------------------------------------------- */
template <class Entity> void CommunicationDescriptor<Entity>::freeRequest() {
const auto & comm = communications.getCommunicator();
// comm.test(communication.request());
comm.freeCommunicationRequest(communication.request());
communications.decrementPending(tag, communication.type());
}
/* -------------------------------------------------------------------------- */
template <class Entity>
-Array<Entity> & CommunicationDescriptor<Entity>::getScheme() {
+const Array<Entity> & CommunicationDescriptor<Entity>::getScheme() {
return scheme;
}
template <class Entity>
void CommunicationDescriptor<Entity>::resetBuffer() {
this->communication.buffer().reset();
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void CommunicationDescriptor<Entity>::packData(
const DataAccessor<Entity> & accessor) {
AKANTU_DEBUG_ASSERT(
communication.type() == _send,
"Cannot pack data on communication that is not of type _send");
accessor.packData(communication.buffer(), scheme, tag);
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void CommunicationDescriptor<Entity>::unpackData(
DataAccessor<Entity> & accessor) {
AKANTU_DEBUG_ASSERT(
communication.type() == _recv,
"Cannot unpack data from communication that is not of type _recv");
accessor.unpackData(communication.buffer(), scheme, tag);
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void CommunicationDescriptor<Entity>::postSend(int hash_id) {
AKANTU_DEBUG_ASSERT(communication.type() == _send,
"Cannot send a communication that is not of type _send");
Tag comm_tag = Tag::genTag(rank, counter, tag, hash_id);
AKANTU_DEBUG_ASSERT(communication.buffer().getPackedSize() ==
communication.size(),
"a problem have been introduced with "
<< "false sent sizes declaration "
<< communication.buffer().getPackedSize()
<< " != " << communication.size());
AKANTU_DEBUG_INFO("Posting send to proc " << proc << " (tag: " << tag << " - "
<< communication.size()
<< " data to send) "
<< " [ " << comm_tag << " ]");
CommunicationRequest & request = communication.request();
request = communications.getCommunicator().asyncSend(communication.buffer(),
proc, comm_tag);
communications.incrementPending(tag, communication.type());
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void CommunicationDescriptor<Entity>::postRecv(int hash_id) {
AKANTU_DEBUG_ASSERT(communication.type() == _recv,
"Cannot receive data for communication ("
<< communication.type()
<< ")that is not of type _recv");
Tag comm_tag = Tag::genTag(proc, counter, tag, hash_id);
AKANTU_DEBUG_INFO("Posting receive from proc "
<< proc << " (tag: " << tag << " - " << communication.size()
<< " data to receive) "
<< " [ " << comm_tag << " ]");
CommunicationRequest & request = communication.request();
request = communications.getCommunicator().asyncReceive(
communication.buffer(), proc, comm_tag);
communications.incrementPending(tag, communication.type());
}
} // akantu
#endif /* __AKANTU_COMMUNICATION_DESCRIPTOR_TMPL_HH__ */
diff --git a/src/synchronizer/element_synchronizer.cc b/src/synchronizer/element_synchronizer.cc
index 76a26d983..25a7f77d1 100644
--- a/src/synchronizer/element_synchronizer.cc
+++ b/src/synchronizer/element_synchronizer.cc
@@ -1,320 +1,349 @@
/**
* @file element_synchronizer.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Dana Christen <dana.christen@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Sep 01 2010
* @date last modification: Fri Jan 22 2016
*
* @brief implementation of a communicator using a static_communicator for
* real
* send/receive
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "element_synchronizer.hh"
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_utils.hh"
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include <iostream>
#include <map>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ElementSynchronizer::ElementSynchronizer(Mesh & mesh, const ID & id,
MemoryID memory_id,
bool register_to_event_manager,
EventHandlerPriority event_priority)
: SynchronizerImpl<Element>(mesh.getCommunicator(), id, memory_id),
mesh(mesh), element_to_prank("element_to_prank", id, memory_id) {
AKANTU_DEBUG_IN();
if (register_to_event_manager)
this->mesh.registerEventHandler(*this, event_priority);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
ElementSynchronizer::~ElementSynchronizer() = default;
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::substituteElements(
const std::map<Element, Element> & old_to_new_elements) {
- auto found_element_end =
- old_to_new_elements.end();
+ auto found_element_end = old_to_new_elements.end();
// substitute old elements with new ones
for (auto && sr : iterate_send_recv) {
for (auto && scheme_pair : communications.iterateSchemes(sr)) {
auto & list = scheme_pair.second;
for (auto & el : list) {
auto found_element_it = old_to_new_elements.find(el);
if (found_element_it != found_element_end)
el = found_element_it->second;
}
}
}
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::onElementsChanged(
const Array<Element> & old_elements_list,
const Array<Element> & new_elements_list, const ElementTypeMapArray<UInt> &,
const ChangedElementsEvent &) {
// create a map to link old elements to new ones
std::map<Element, Element> old_to_new_elements;
for (UInt el = 0; el < old_elements_list.size(); ++el) {
AKANTU_DEBUG_ASSERT(old_to_new_elements.find(old_elements_list(el)) ==
old_to_new_elements.end(),
"The same element cannot appear twice in the list");
old_to_new_elements[old_elements_list(el)] = new_elements_list(el);
}
substituteElements(old_to_new_elements);
communications.invalidateSizes();
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::onElementsRemoved(
const Array<Element> & element_to_remove,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent &) {
AKANTU_DEBUG_IN();
this->removeElements(element_to_remove);
this->renumberElements(new_numbering);
communications.invalidateSizes();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::buildElementToPrank() {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
element_to_prank.initialize(mesh, _spatial_dimension = spatial_dimension,
_element_kind = _ek_not_defined,
_with_nb_element = true, _default_value = rank);
/// assign prank to all ghost elements
for (auto && scheme : communications.iterateSchemes(_recv)) {
auto & recv = scheme.second;
auto proc = scheme.first;
for (auto & element : recv) {
element_to_prank(element) = proc;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Int ElementSynchronizer::getRank(const Element & element) const {
- if(not element_to_prank.exists(element.type, element.ghost_type)) {
+ if (not element_to_prank.exists(element.type, element.ghost_type)) {
// Nicolas: Ok This is nasty I know....
const_cast<ElementSynchronizer *>(this)->buildElementToPrank();
}
return element_to_prank(element);
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::reset() {
AKANTU_DEBUG_IN();
communications.resetSchemes();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::removeElements(
const Array<Element> & element_to_remove) {
AKANTU_DEBUG_IN();
std::vector<CommunicationRequest> send_requests;
std::map<UInt, Array<UInt>> list_of_elements_per_proc;
auto filter_list = [](const Array<UInt> & keep, Array<Element> & list) {
Array<Element> new_list;
for (UInt e = 0; e < keep.size() - 1; ++e) {
Element & el = list(keep(e));
new_list.push_back(el);
}
list.copy(new_list);
};
// Handling ghost elements
for (auto && scheme_pair : communications.iterateSchemes(_recv)) {
auto proc = scheme_pair.first;
auto & recv = scheme_pair.second;
Array<UInt> & keep_list = list_of_elements_per_proc[proc];
auto rem_it = element_to_remove.begin();
auto rem_end = element_to_remove.end();
for (auto && pair : enumerate(recv)) {
const auto & el = std::get<1>(pair);
auto pos = std::find(rem_it, rem_end, el);
if (pos == rem_end) {
keep_list.push_back(std::get<0>(pair));
}
}
keep_list.push_back(UInt(-1)); // To no send empty arrays
auto && tag = Tag::genTag(proc, 0, Tag::_MODIFY_SCHEME, this->hash_id);
AKANTU_DEBUG_INFO("Sending a message of size "
<< keep_list.size() << " to proc " << proc << " TAG("
<< tag << ")");
send_requests.push_back(this->communicator.asyncSend(keep_list, proc, tag));
auto old_size = recv.size();
filter_list(keep_list, recv);
AKANTU_DEBUG_INFO("I had " << old_size << " elements to recv from proc "
<< proc << " and " << keep_list.size()
<< " elements to keep. I have " << recv.size()
<< " elements left.");
}
for (auto && scheme_pair : communications.iterateSchemes(_send)) {
auto proc = scheme_pair.first;
auto & send = scheme_pair.second;
CommunicationStatus status;
auto && tag = Tag::genTag(rank, 0, Tag::_MODIFY_SCHEME, this->hash_id);
AKANTU_DEBUG_INFO("Getting number of elements of proc "
<< proc << " to keep - TAG(" << tag << ")");
this->communicator.probe<UInt>(proc, tag, status);
Array<UInt> keep_list(status.size());
AKANTU_DEBUG_INFO("Receiving list of elements ("
<< keep_list.size() << " elements) to keep for proc "
<< proc << " TAG(" << tag << ")");
this->communicator.receive(keep_list, proc, tag);
auto old_size = send.size();
filter_list(keep_list, send);
AKANTU_DEBUG_INFO("I had " << old_size << " elements to send to proc "
<< proc << " and " << keep_list.size()
<< " elements to keep. I have " << send.size()
<< " elements left.");
}
this->communicator.waitAll(send_requests);
this->communicator.freeCommunicationRequest(send_requests);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::renumberElements(
const ElementTypeMapArray<UInt> & new_numbering) {
- for(auto && sr : iterate_send_recv) {
+ for (auto && sr : iterate_send_recv) {
for (auto && scheme_pair : communications.iterateSchemes(sr)) {
auto & list = scheme_pair.second;
for (auto && el : list) {
- if(new_numbering.exists(el.type, el.ghost_type))
+ if (new_numbering.exists(el.type, el.ghost_type))
el.element = new_numbering(el);
}
}
}
}
/* -------------------------------------------------------------------------- */
UInt ElementSynchronizer::sanityCheckDataSize(
const Array<Element> & elements, const SynchronizationTag &) const {
- return (elements.size() * mesh.getSpatialDimension() * sizeof(Real) +
- sizeof(SynchronizationTag));
+ UInt size = 0;
+ size += sizeof(SynchronizationTag); // tag
+ size += sizeof(UInt); // comm_desc.getNbData();
+ size += sizeof(UInt); // comm_desc.getProc();
+ size += sizeof(UInt); // mesh.getCommunicator().whoAmI();
+
+ // barycenters
+ size += (elements.size() * mesh.getSpatialDimension() * sizeof(Real));
+ return size;
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::packSanityCheckData(
CommunicationDescriptor<Element> & comm_desc) const {
auto & buffer = comm_desc.getBuffer();
buffer << comm_desc.getTag();
+ buffer << comm_desc.getNbData();
+ buffer << comm_desc.getProc();
+ buffer << mesh.getCommunicator().whoAmI();
auto & send_element = comm_desc.getScheme();
/// pack barycenters in debug mode
- for(auto && element : send_element) {
+ for (auto && element : send_element) {
Vector<Real> barycenter(mesh.getSpatialDimension());
mesh.getBarycenter(element, barycenter);
buffer << barycenter;
}
}
/* -------------------------------------------------------------------------- */
void ElementSynchronizer::unpackSanityCheckData(
CommunicationDescriptor<Element> & comm_desc) const {
auto & buffer = comm_desc.getBuffer();
const auto & tag = comm_desc.getTag();
+ auto nb_data = comm_desc.getNbData();
+ auto proc = comm_desc.getProc();
+ auto rank = mesh.getCommunicator().whoAmI();
+
+ decltype(nb_data) recv_nb_data;
+ decltype(proc) recv_proc;
+ decltype(rank) recv_rank;
+
SynchronizationTag t;
buffer >> t;
+ buffer >> recv_nb_data;
+ buffer >> recv_proc;
+ buffer >> recv_rank;
AKANTU_DEBUG_ASSERT(
t == tag, "The tag received does not correspond to the tag expected");
+ AKANTU_DEBUG_ASSERT(
+ nb_data == recv_nb_data,
+ "The nb_data received does not correspond to the nb_data expected");
+
+ AKANTU_DEBUG_ASSERT(UInt(recv_rank) == proc,
+ "The rank received does not correspond to the proc");
+
+ AKANTU_DEBUG_ASSERT(recv_proc == UInt(rank),
+ "The proc received does not correspond to the rank");
+
auto & recv_element = comm_desc.getScheme();
auto spatial_dimension = mesh.getSpatialDimension();
- for(auto && element : recv_element) {
+ for (auto && element : recv_element) {
Vector<Real> barycenter_loc(spatial_dimension);
mesh.getBarycenter(element, barycenter_loc);
Vector<Real> barycenter(spatial_dimension);
buffer >> barycenter;
- for (UInt i = 0; i < spatial_dimension; ++i) {
- if (!Math::are_float_equal(barycenter_loc(i), barycenter(i)))
- AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: "
- << element << "(barycenter[" << i
- << "] = " << barycenter_loc(i) << " and buffer[" << i
- << "] = " << barycenter(i) << ") ["
- << std::abs(barycenter(i) - barycenter_loc(i))
- << "] - tag: " << tag);
- }
+
+ auto dist = barycenter_loc.distance(barycenter);
+ if (not Math::are_float_equal(dist, 0.))
+ AKANTU_EXCEPTION("Unpacking an unknown value for the element "
+ << element << "(barycenter " << barycenter_loc
+ << " != buffer " << barycenter << ") [" << dist
+ << "] - tag: " << tag << " comm from " << proc << " to "
+ << rank);
}
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/synchronizer/master_element_info_per_processor.cc b/src/synchronizer/master_element_info_per_processor.cc
index ff2fced01..8a4759d6d 100644
--- a/src/synchronizer/master_element_info_per_processor.cc
+++ b/src/synchronizer/master_element_info_per_processor.cc
@@ -1,462 +1,461 @@
/**
* @file master_element_info_per_processor.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Mar 11 14:57:13 2016
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_iterators.hh"
+#include "communicator.hh"
#include "element_group.hh"
#include "element_info_per_processor.hh"
#include "element_synchronizer.hh"
#include "mesh_iterators.hh"
#include "mesh_utils.hh"
-#include "communicator.hh"
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include <iostream>
#include <map>
#include <tuple>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
MasterElementInfoPerProc::MasterElementInfoPerProc(
ElementSynchronizer & synchronizer, UInt message_cnt, UInt root,
ElementType type, const MeshPartition & partition)
: ElementInfoPerProc(synchronizer, message_cnt, root, type),
partition(partition), all_nb_local_element(nb_proc, 0),
all_nb_ghost_element(nb_proc, 0), all_nb_element_to_send(nb_proc, 0) {
Vector<UInt> size(5);
size(0) = (UInt)type;
if (type != _not_defined) {
nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
nb_element = mesh.getNbElement(type);
const auto & partition_num =
this->partition.getPartition(this->type, _not_ghost);
const auto & ghost_partition =
this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
for (UInt el = 0; el < nb_element; ++el) {
this->all_nb_local_element[partition_num(el)]++;
for (auto part = ghost_partition.begin(el);
part != ghost_partition.end(el); ++part) {
this->all_nb_ghost_element[*part]++;
}
this->all_nb_element_to_send[partition_num(el)] +=
ghost_partition.getNbCols(el) + 1;
}
/// tag info
std::vector<std::string> tag_names;
this->getMeshData().getTagNames(tag_names, type);
this->nb_tags = tag_names.size();
size(4) = nb_tags;
for (UInt p = 0; p < nb_proc; ++p) {
if (p != root) {
size(1) = this->all_nb_local_element[p];
size(2) = this->all_nb_ghost_element[p];
size(3) = this->all_nb_element_to_send[p];
AKANTU_DEBUG_INFO(
"Sending connectivities informations to proc "
<< p << " TAG("
<< Tag::genTag(this->rank, this->message_count, Tag::_SIZES)
<< ")");
comm.send(size, p,
Tag::genTag(this->rank, this->message_count, Tag::_SIZES));
} else {
this->nb_local_element = this->all_nb_local_element[p];
this->nb_ghost_element = this->all_nb_ghost_element[p];
}
}
} else {
for (UInt p = 0; p < this->nb_proc; ++p) {
if (p != this->root) {
AKANTU_DEBUG_INFO(
"Sending empty connectivities informations to proc "
<< p << " TAG("
<< Tag::genTag(this->rank, this->message_count, Tag::_SIZES)
<< ")");
comm.send(size, p,
Tag::genTag(this->rank, this->message_count, Tag::_SIZES));
}
}
}
}
/* ------------------------------------------------------------------------ */
void MasterElementInfoPerProc::synchronizeConnectivities() {
const auto & partition_num =
this->partition.getPartition(this->type, _not_ghost);
const auto & ghost_partition =
this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
std::vector<Array<UInt>> buffers(this->nb_proc);
- auto conn_it = this->mesh.getConnectivity(this->type, _not_ghost)
- .begin(this->nb_nodes_per_element);
- auto conn_end = this->mesh.getConnectivity(this->type, _not_ghost)
- .end(this->nb_nodes_per_element);
+ const auto & connectivities =
+ this->mesh.getConnectivity(this->type, _not_ghost);
/// copying the local connectivity
- auto part_it = partition_num.begin();
- for (; conn_it != conn_end; ++conn_it, ++part_it) {
- const auto & conn = *conn_it;
+ for (auto && part_conn :
+ zip(partition_num,
+ make_view(connectivities, this->nb_nodes_per_element))) {
+ auto && part = std::get<0>(part_conn);
+ auto && conn = std::get<1>(part_conn);
for (UInt i = 0; i < conn.size(); ++i) {
- buffers[*part_it].push_back(conn[i]);
+ buffers[part].push_back(conn[i]);
}
}
/// copying the connectivity of ghost element
- conn_it = this->mesh.getConnectivity(this->type, _not_ghost)
- .begin(this->nb_nodes_per_element);
- for (UInt el = 0; conn_it != conn_end; ++el, ++conn_it) {
+ for (auto && tuple :
+ enumerate(make_view(connectivities, this->nb_nodes_per_element))) {
+ auto && el = std::get<0>(tuple);
+ auto && conn = std::get<1>(tuple);
for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el);
++part) {
UInt proc = *part;
- const Vector<UInt> & conn = *conn_it;
for (UInt i = 0; i < conn.size(); ++i) {
buffers[proc].push_back(conn[i]);
}
}
}
#ifndef AKANTU_NDEBUG
- for (UInt p = 0; p < this->nb_proc; ++p) {
+ for (auto p : arange(this->nb_proc)) {
UInt size = this->nb_nodes_per_element *
(this->all_nb_local_element[p] + this->all_nb_ghost_element[p]);
AKANTU_DEBUG_ASSERT(
buffers[p].size() == size,
"The connectivity data packed in the buffer are not correct");
}
#endif
/// send all connectivity and ghost information to all processors
std::vector<CommunicationRequest> requests;
- for (UInt p = 0; p < this->nb_proc; ++p) {
+ for (auto p : arange(this->nb_proc)) {
if (p == this->root)
continue;
auto && tag =
Tag::genTag(this->rank, this->message_count, Tag::_CONNECTIVITY);
AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG(" << tag
<< ")");
requests.push_back(comm.asyncSend(buffers[p], p, tag));
}
Array<UInt> & old_nodes = this->getNodesGlobalIds();
/// create the renumbered connectivity
AKANTU_DEBUG_INFO("Renumbering local connectivities");
MeshUtils::renumberMeshNodes(mesh, buffers[root], all_nb_local_element[root],
all_nb_ghost_element[root], type, old_nodes);
-
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
}
/* ------------------------------------------------------------------------ */
void MasterElementInfoPerProc::synchronizePartitions() {
const auto & partition_num =
this->partition.getPartition(this->type, _not_ghost);
const auto & ghost_partition =
this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
std::vector<Array<UInt>> buffers(this->partition.getNbPartition());
/// splitting the partition information to send them to processors
Vector<UInt> count_by_proc(nb_proc, 0);
for (UInt el = 0; el < nb_element; ++el) {
UInt proc = partition_num(el);
buffers[proc].push_back(ghost_partition.getNbCols(el));
UInt i(0);
for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el);
++part, ++i) {
buffers[proc].push_back(*part);
}
}
for (UInt el = 0; el < nb_element; ++el) {
UInt i(0);
for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el);
++part, ++i) {
buffers[*part].push_back(partition_num(el));
}
}
#ifndef AKANTU_NDEBUG
for (UInt p = 0; p < this->nb_proc; ++p) {
AKANTU_DEBUG_ASSERT(buffers[p].size() == (this->all_nb_ghost_element[p] +
this->all_nb_element_to_send[p]),
"Data stored in the buffer are most probably wrong");
}
#endif
std::vector<CommunicationRequest> requests;
/// last data to compute the communication scheme
for (UInt p = 0; p < this->nb_proc; ++p) {
if (p == this->root)
continue;
auto && tag =
Tag::genTag(this->rank, this->message_count, Tag::_PARTITIONS);
AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG("
<< tag << ")");
requests.push_back(comm.asyncSend(buffers[p], p, tag));
}
if (Mesh::getSpatialDimension(this->type) ==
this->mesh.getSpatialDimension()) {
AKANTU_DEBUG_INFO("Creating communications scheme");
this->fillCommunicationScheme(buffers[this->rank]);
}
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
}
/* -------------------------------------------------------------------------- */
void MasterElementInfoPerProc::synchronizeTags() {
AKANTU_DEBUG_IN();
if (this->nb_tags == 0) {
AKANTU_DEBUG_OUT();
return;
}
UInt mesh_data_sizes_buffer_length;
auto & mesh_data = this->getMeshData();
/// tag info
std::vector<std::string> tag_names;
mesh_data.getTagNames(tag_names, type);
// Make sure the tags are sorted (or at least not in random order),
// because they come from a map !!
std::sort(tag_names.begin(), tag_names.end());
// Sending information about the tags in mesh_data: name, data type and
// number of components of the underlying array associated to the current
// type
DynamicCommunicationBuffer mesh_data_sizes_buffer;
for (auto && tag_name : tag_names) {
mesh_data_sizes_buffer << tag_name;
mesh_data_sizes_buffer << mesh_data.getTypeCode(tag_name);
mesh_data_sizes_buffer << mesh_data.getNbComponent(tag_name, type);
}
mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.size();
AKANTU_DEBUG_INFO(
"Broadcasting the size of the information about the mesh data tags: ("
<< mesh_data_sizes_buffer_length << ").");
comm.broadcast(mesh_data_sizes_buffer_length, root);
AKANTU_DEBUG_INFO(
"Broadcasting the information about the mesh data tags, addr "
<< (void *)mesh_data_sizes_buffer.storage());
if (mesh_data_sizes_buffer_length != 0)
comm.broadcast(mesh_data_sizes_buffer, root);
if (mesh_data_sizes_buffer_length != 0) {
// Sending the actual data to each processor
std::vector<DynamicCommunicationBuffer> buffers(nb_proc);
// Loop over each tag for the current type
for (auto && tag_name : tag_names) {
// Type code of the current tag (i.e. the tag named *names_it)
this->fillTagBuffer(buffers, tag_name);
}
std::vector<CommunicationRequest> requests;
for (UInt p = 0; p < nb_proc; ++p) {
if (p == root)
continue;
auto && tag =
Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA);
AKANTU_DEBUG_INFO("Sending " << buffers[p].size()
<< " bytes of mesh data to proc " << p
<< " TAG(" << tag << ")");
requests.push_back(comm.asyncSend(buffers[p], p, tag));
}
// Loop over each tag for the current type
for (auto && tag_name : tag_names) {
// Reinitializing the mesh data on the master
this->fillMeshData(buffers[root], tag_name,
mesh_data.getTypeCode(tag_name),
mesh_data.getNbComponent(tag_name, type));
}
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
requests.clear();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void MasterElementInfoPerProc::fillTagBufferTemplated(
std::vector<DynamicCommunicationBuffer> & buffers,
const std::string & tag_name) {
MeshData & mesh_data = this->getMeshData();
const auto & data = mesh_data.getElementalDataArray<T>(tag_name, type);
const auto & partition_num =
this->partition.getPartition(this->type, _not_ghost);
const auto & ghost_partition =
this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
// Not possible to use the iterator because it potentially triggers the
// creation of complex
// type templates (such as akantu::Vector< std::vector<Element> > which don't
// implement the right interface
// (e.g. operator<< in that case).
// typename Array<T>::template const_iterator< Vector<T> > data_it =
// data.begin(data.getNbComponent());
// typename Array<T>::template const_iterator< Vector<T> > data_end =
// data.end(data.getNbComponent());
const T * data_it = data.storage();
const T * data_end = data.storage() + data.size() * data.getNbComponent();
const UInt * part = partition_num.storage();
/// copying the data, element by element
for (; data_it != data_end; ++part) {
for (UInt j(0); j < data.getNbComponent(); ++j, ++data_it) {
buffers[*part] << *data_it;
}
}
data_it = data.storage();
/// copying the data for the ghost element
for (UInt el(0); data_it != data_end;
data_it += data.getNbComponent(), ++el) {
auto it = ghost_partition.begin(el);
auto end = ghost_partition.end(el);
for (; it != end; ++it) {
UInt proc = *it;
for (UInt j(0); j < data.getNbComponent(); ++j) {
buffers[proc] << data_it[j];
}
}
}
}
/* -------------------------------------------------------------------------- */
void MasterElementInfoPerProc::fillTagBuffer(
std::vector<DynamicCommunicationBuffer> & buffers,
const std::string & tag_name) {
MeshData & mesh_data = this->getMeshData();
#define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \
case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \
this->fillTagBufferTemplated<BOOST_PP_TUPLE_ELEM(2, 1, elem)>(buffers, \
tag_name); \
break; \
}
MeshDataTypeCode data_type_code = mesh_data.getTypeCode(tag_name);
switch (data_type_code) {
BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, ,
AKANTU_MESH_DATA_TYPES)
default:
AKANTU_DEBUG_ERROR("Could not obtain the type of tag" << tag_name << "!");
break;
}
#undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA
}
/* -------------------------------------------------------------------------- */
void MasterElementInfoPerProc::synchronizeGroups() {
AKANTU_DEBUG_IN();
std::vector<DynamicCommunicationBuffer> buffers(nb_proc);
using ElementToGroup = std::vector<std::vector<std::string>>;
ElementToGroup element_to_group(nb_element);
for (auto & eg : ElementGroupsIterable(mesh)) {
const auto & name = eg.getName();
for (const auto & element : eg.getElements(type, _not_ghost)) {
element_to_group[element].push_back(name);
}
auto eit = eg.begin(type, _not_ghost);
if (eit != eg.end(type, _not_ghost))
const_cast<Array<UInt> &>(eg.getElements(type)).empty();
}
const auto & partition_num =
this->partition.getPartition(this->type, _not_ghost);
const auto & ghost_partition =
this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
/// copying the data, element by element
for (auto && pair : zip(partition_num, element_to_group)) {
buffers[std::get<0>(pair)] << std::get<1>(pair);
}
/// copying the data for the ghost element
for (auto && pair : enumerate(element_to_group)) {
auto && el = std::get<0>(pair);
auto it = ghost_partition.begin(el);
auto end = ghost_partition.end(el);
for (; it != end; ++it) {
UInt proc = *it;
buffers[proc] << std::get<1>(pair);
}
}
std::vector<CommunicationRequest> requests;
for (UInt p = 0; p < this->nb_proc; ++p) {
if (p == this->rank)
continue;
auto && tag = Tag::genTag(this->rank, p, Tag::_ELEMENT_GROUP);
AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << tag
<< ")");
requests.push_back(comm.asyncSend(buffers[p], p, tag));
}
this->fillElementGroupsFromBuffer(buffers[this->rank]);
comm.waitAll(requests);
comm.freeCommunicationRequest(requests);
requests.clear();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/synchronizer/synchronizer_impl_tmpl.hh b/src/synchronizer/synchronizer_impl_tmpl.hh
index 7f3037156..8cd7a2cdf 100644
--- a/src/synchronizer/synchronizer_impl_tmpl.hh
+++ b/src/synchronizer/synchronizer_impl_tmpl.hh
@@ -1,305 +1,305 @@
/**
* @file synchronizer_impl_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Aug 24 13:29:47 2016
*
* @brief Implementation of the SynchronizerImpl
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "synchronizer_impl.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__
#define __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__
namespace akantu {
/* -------------------------------------------------------------------------- */
template <class Entity>
SynchronizerImpl<Entity>::SynchronizerImpl(const Communicator & comm, const ID & id, MemoryID memory_id)
: Synchronizer(comm, id, memory_id), communications(comm) {}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::synchronizeOnceImpl(
DataAccessor<Entity> & data_accessor,
const SynchronizationTag & tag) const {
// no need to synchronize
if (this->nb_proc == 1)
return;
using CommunicationRequests = std::vector<CommunicationRequest>;
using CommunicationBuffers = std::map<UInt, CommunicationBuffer>;
CommunicationRequests send_requests, recv_requests;
CommunicationBuffers send_buffers, recv_buffers;
auto postComm = [&](const CommunicationSendRecv & sr,
CommunicationBuffers & buffers,
CommunicationRequests & requests) -> void {
for (auto && pair : communications.iterateSchemes(sr)) {
auto & proc = pair.first;
const auto & scheme = pair.second;
auto & buffer = buffers[proc];
UInt buffer_size = data_accessor.getNbData(scheme, tag);
buffer.resize(buffer_size);
if (sr == _recv) {
requests.push_back(communicator.asyncReceive(
buffer, proc,
Tag::genTag(this->rank, 0, Tag::_SYNCHRONIZE, this->hash_id)));
} else {
data_accessor.packData(buffer, scheme, tag);
send_requests.push_back(communicator.asyncSend(
buffer, proc,
Tag::genTag(proc, 0, Tag::_SYNCHRONIZE, this->hash_id)));
}
}
};
// post the receive requests
postComm(_recv, recv_buffers, recv_requests);
// post the send data requests
postComm(_send, send_buffers, send_requests);
// treat the receive requests
UInt request_ready;
while ((request_ready = communicator.waitAny(recv_requests)) != UInt(-1)) {
CommunicationRequest & req = recv_requests[request_ready];
UInt proc = req.getSource();
CommunicationBuffer & buffer = recv_buffers[proc];
const auto & scheme = this->communications.getScheme(proc, _recv);
data_accessor.unpackData(buffer, scheme, tag);
req.free();
recv_requests.erase(recv_requests.begin() + request_ready);
}
communicator.waitAll(send_requests);
communicator.freeCommunicationRequest(send_requests);
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::asynchronousSynchronizeImpl(
const DataAccessor<Entity> & data_accessor,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
if (not this->communications.hasCommunicationSize(tag))
this->computeBufferSize(data_accessor, tag);
this->communications.incrementCounter(tag);
// Posting the receive -------------------------------------------------------
if (this->communications.hasPendingRecv(tag)) {
AKANTU_CUSTOM_EXCEPTION_INFO(
debug::CommunicationException(),
"There must still be some pending receive communications."
<< " Tag is " << tag << " Cannot start new ones");
}
for (auto && comm_desc : this->communications.iterateRecv(tag)) {
comm_desc.postRecv(this->hash_id);
}
// Posting the sends -------------------------------------------------------
if (communications.hasPendingSend(tag)) {
AKANTU_CUSTOM_EXCEPTION_INFO(
debug::CommunicationException(),
"There must be some pending sending communications."
<< " Tag is " << tag);
}
for (auto && comm_desc : this->communications.iterateSend(tag)) {
comm_desc.resetBuffer();
#ifndef AKANTU_NDEBUG
this->packSanityCheckData(comm_desc);
#endif
comm_desc.packData(data_accessor);
comm_desc.postSend(this->hash_id);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::waitEndSynchronizeImpl(
DataAccessor<Entity> & data_accessor, const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
#ifndef AKANTU_NDEBUG
if (this->communications.begin(tag, _recv) !=
this->communications.end(tag, _recv) &&
!this->communications.hasPendingRecv(tag))
AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(),
"No pending communication with the tag \""
<< tag);
#endif
auto recv_end = this->communications.end(tag, _recv);
decltype(recv_end) recv_it;
while ((recv_it = this->communications.waitAnyRecv(tag)) != recv_end) {
auto && comm_desc = *recv_it;
#ifndef AKANTU_NDEBUG
this->unpackSanityCheckData(comm_desc);
#endif
comm_desc.unpackData(data_accessor);
comm_desc.resetBuffer();
comm_desc.freeRequest();
}
this->communications.waitAllSend(tag);
this->communications.freeSendRequests(tag);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::computeAllBufferSizes(
const DataAccessor<Entity> & data_accessor) {
for (auto && tag : this->communications.iterateTags()) {
this->computeBufferSize(data_accessor, tag);
}
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::computeBufferSizeImpl(
const DataAccessor<Entity> & data_accessor,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
if (not this->communications.hasCommunication(tag)) {
this->communications.initializeCommunications(tag);
AKANTU_DEBUG_ASSERT(communications.hasCommunication(tag) == true,
"Communications where not properly initialized");
}
for (auto sr : iterate_send_recv) {
for (auto && pair : this->communications.iterateSchemes(sr)) {
auto proc = pair.first;
- auto & scheme = pair.second;
+ const auto & scheme = pair.second;
UInt size = 0;
#ifndef AKANTU_NDEBUG
size += this->sanityCheckDataSize(scheme, tag);
#endif
size += data_accessor.getNbData(scheme, tag);
AKANTU_DEBUG_INFO("I have "
<< size << "(" << printMemorySize<char>(size) << " - "
<< scheme.size() << " element(s)) data to "
<< std::string(sr == _recv ? "receive from" : "send to")
<< proc << " for tag " << tag);
this->communications.setCommunicationSize(tag, proc, size, sr);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <typename Entity>
template <typename Pred>
void SynchronizerImpl<Entity>::split(SynchronizerImpl<Entity> & in_synchronizer,
Pred && pred) {
AKANTU_DEBUG_IN();
auto filter_list = [&](auto & list, auto & new_list) {
auto copy = list;
list.resize(0);
new_list.resize(0);
for (auto && entity : copy) {
if (std::forward<Pred>(pred)(entity)) {
new_list.push_back(entity);
} else {
list.push_back(entity);
}
}
};
for (auto sr : iterate_send_recv) {
for (auto & scheme_pair :
in_synchronizer.communications.iterateSchemes(sr)) {
auto proc = scheme_pair.first;
auto & scheme = scheme_pair.second;
auto & new_scheme = communications.createScheme(proc, sr);
filter_list(new_scheme, scheme);
}
}
in_synchronizer.communications.invalidateSizes();
communications.invalidateSizes();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <typename Entity>
template <typename Updater>
void SynchronizerImpl<Entity>::updateSchemes(Updater && scheme_updater) {
for (auto sr : iterate_send_recv) {
for (auto & scheme_pair : communications.iterateSchemes(sr)) {
auto proc = scheme_pair.first;
auto & scheme = scheme_pair.second;
std::forward<Updater>(scheme_updater)(scheme, proc, sr);
}
}
communications.invalidateSizes();
}
/* -------------------------------------------------------------------------- */
template <class Entity>
UInt SynchronizerImpl<Entity>::sanityCheckDataSize(
const Array<Entity> &, const SynchronizationTag &) const {
return 0;
}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::packSanityCheckData(
CommunicationDescriptor<Entity> &) const {}
/* -------------------------------------------------------------------------- */
template <class Entity>
void SynchronizerImpl<Entity>::unpackSanityCheckData(
CommunicationDescriptor<Entity> &) const {}
/* -------------------------------------------------------------------------- */
} // namespace akantu
#endif /* __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ */