diff --git a/src/common/aka_array.cc b/src/common/aka_array.cc
index f8fa5c317..31bf3f4fb 100644
--- a/src/common/aka_array.cc
+++ b/src/common/aka_array.cc
@@ -1,121 +1,123 @@
/**
* @file aka_array.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Tue Aug 18 2015
*
* @brief Implementation of akantu::Array
*
* @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 <memory>
+#include <utility>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_array.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Functions ArrayBase */
/* -------------------------------------------------------------------------- */
-ArrayBase::ArrayBase(const ID & id)
- : id(id), allocated_size(0), size(0), nb_component(1), size_of_type(0) {}
+ArrayBase::ArrayBase(ID id)
+ : id(std::move(id)), allocated_size(0), size(0), nb_component(1),
+ size_of_type(0) {}
/* -------------------------------------------------------------------------- */
-ArrayBase::~ArrayBase() {}
+ArrayBase::~ArrayBase() = default;
/* -------------------------------------------------------------------------- */
void ArrayBase::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "ArrayBase [" << std::endl;
stream << space << " + size : " << size << std::endl;
stream << space << " + nb component : " << nb_component << std::endl;
stream << space << " + allocated size : " << allocated_size << std::endl;
Real mem_size = (allocated_size * nb_component * size_of_type) / 1024.;
stream << space << " + size of type : " << size_of_type << "B"
<< std::endl;
stream << space << " + memory allocated : " << mem_size << "kB" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
template <> Int Array<Real>::find(const Real & elem) const {
AKANTU_DEBUG_IN();
UInt i = 0;
Real epsilon = std::numeric_limits<Real>::epsilon();
for (; (i < size) && (fabs(values[i] - elem) <= epsilon); ++i)
;
AKANTU_DEBUG_OUT();
return (i == size) ? -1 : (Int)i;
}
/* -------------------------------------------------------------------------- */
template <>
Array<ElementType> & Array<ElementType>::operator*=(__attribute__((unused))
const ElementType & alpha) {
AKANTU_DEBUG_TO_IMPLEMENT();
return *this;
}
template <>
Array<ElementType> & Array<ElementType>::
operator-=(__attribute__((unused)) const Array<ElementType> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
return *this;
}
template <>
Array<ElementType> & Array<ElementType>::
operator+=(__attribute__((unused)) const Array<ElementType> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
return *this;
}
template <>
Array<char> & Array<char>::operator*=(__attribute__((unused))
const char & alpha) {
AKANTU_DEBUG_TO_IMPLEMENT();
return *this;
}
template <>
Array<char> & Array<char>::operator-=(__attribute__((unused))
const Array<char> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
return *this;
}
template <>
Array<char> & Array<char>::operator+=(__attribute__((unused))
const Array<char> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
return *this;
}
__END_AKANTU__
diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh
index eca903e93..557051299 100644
--- a/src/common/aka_array.hh
+++ b/src/common/aka_array.hh
@@ -1,397 +1,392 @@
/**
* @file aka_array.hh
*
* @author Till Junge <till.junge@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Jan 22 2016
*
* @brief Array container for Akantu
* This container differs from the std::vector from the fact it as 2 dimensions
* a main dimension and the size stored per entries
*
* @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_VECTOR_HH__
#define __AKANTU_VECTOR_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
#include <typeinfo>
#include <vector>
/* -------------------------------------------------------------------------- */
-__BEGIN_AKANTU__
+namespace akantu {
/// class that afford to store vectors in static memory
class ArrayBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- ArrayBase(const ID & id = "");
+ ArrayBase(ID id = "");
virtual ~ArrayBase();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// get the amount of space allocated in bytes
inline UInt getMemorySize() const;
/// set the size to zero without freeing the allocated space
inline void empty();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// Get the real size allocated in memory
AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt);
/// Get the Size of the Array
AKANTU_GET_MACRO(Size, size, UInt);
/// Get the number of components
AKANTU_GET_MACRO(NbComponent, nb_component, UInt);
/// Get the name of th array
AKANTU_GET_MACRO(ID, id, const ID &);
/// Set the name of th array
AKANTU_SET_MACRO(ID, id, const ID &);
// AKANTU_GET_MACRO(Tag, tag, const std::string &);
// AKANTU_SET_MACRO(Tag, tag, const std::string &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// id of the vector
ID id;
/// the size allocated
- UInt allocated_size;
+ UInt allocated_size{0};
/// the size used
- UInt size;
+ UInt size{0};
/// number of components
- UInt nb_component;
+ UInt nb_component{1};
/// size of the stored type
- UInt size_of_type;
-
- // /// User defined tag
- // std::string tag;
+ UInt size_of_type{0};
};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <typename T, bool is_scal> class Array : public ArrayBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- typedef T value_type;
- typedef value_type & reference;
- typedef value_type * pointer_type;
- typedef const value_type & const_reference;
+ using value_type = T;
+ using reference = value_type &;
+ using pointer_type = value_type *;
+ using const_reference = const value_type &;
/// Allocation of a new vector
inline Array(UInt size = 0, UInt nb_component = 1, const ID & id = "");
/// Allocation of a new vector with a default value
Array(UInt size, UInt nb_component, const value_type def_values[],
const ID & id = "");
/// Allocation of a new vector with a default value
Array(UInt size, UInt nb_component, const_reference value,
const ID & id = "");
/// Copy constructor (deep copy if deep=true)
Array(const Array<value_type, is_scal> & vect, bool deep = true,
const ID & id = "");
#ifndef SWIG
/// Copy constructor (deep copy)
Array(const std::vector<value_type> & vect);
#endif
- virtual inline ~Array();
+ inline ~Array() override;
Array & operator=(const Array & a) {
/// this is to let STL allocate and copy arrays in the case of
/// std::vector::resize
AKANTU_DEBUG_ASSERT(this->size == 0, "Cannot copy akantu::Array");
return const_cast<Array &>(a);
}
/* ------------------------------------------------------------------------ */
/* Iterator */
/* ------------------------------------------------------------------------ */
/// \todo protected: does not compile with intel check why
public:
template <class R, class IR = R, bool issame = is_same<IR, T>::value>
class iterator_internal;
public:
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
template <typename R = T> class const_iterator;
template <typename R = T> class iterator;
/* ------------------------------------------------------------------------ */
/// iterator for Array of nb_component = 1
- typedef iterator<T> scalar_iterator;
+ using scalar_iterator = iterator<T>;
/// const_iterator for Array of nb_component = 1
- typedef const_iterator<T> const_scalar_iterator;
+ using const_scalar_iterator = const_iterator<T>;
/// iterator rerturning Vectors of size n on entries of Array with
/// nb_component = n
- typedef iterator<Vector<T> > vector_iterator;
+ using vector_iterator = iterator<Vector<T>>;
/// const_iterator rerturning Vectors of n size on entries of Array with
/// nb_component = n
- typedef const_iterator<Vector<T> > const_vector_iterator;
+ using const_vector_iterator = const_iterator<Vector<T>>;
/// iterator rerturning Matrices of size (m, n) on entries of Array with
/// nb_component = m*n
- typedef iterator<Matrix<T> > matrix_iterator;
+ using matrix_iterator = iterator<Matrix<T>>;
/// const iterator rerturning Matrices of size (m, n) on entries of Array with
/// nb_component = m*n
- typedef const_iterator<Matrix<T> > const_matrix_iterator;
+ using const_matrix_iterator = const_iterator<Matrix<T>>;
/* ------------------------------------------------------------------------ */
-
/// Get an iterator that behaves like a pointer T * to the first entry
inline scalar_iterator begin();
/// Get an iterator that behaves like a pointer T * to the end of the Array
inline scalar_iterator end();
/// Get a const_iterator to the beginging of an Array of scalar
inline const_scalar_iterator begin() const;
/// Get a const_iterator to the end of an Array of scalar
inline const_scalar_iterator end() const;
- /// Get a scalar_iterator on the beginning of the Array considered of shape (new_size)
+ /// Get a scalar_iterator on the beginning of the Array considered of shape
+ /// (new_size)
inline scalar_iterator begin_reinterpret(UInt new_size);
- /// Get a scalar_iterator on the end of the Array considered of shape (new_size)
+ /// Get a scalar_iterator on the end of the Array considered of shape
+ /// (new_size)
inline scalar_iterator end_reinterpret(UInt new_size);
- /// Get a const_scalar_iterator on the beginning of the Array considered of shape (new_size)
+ /// Get a const_scalar_iterator on the beginning of the Array considered of
+ /// shape (new_size)
inline const_scalar_iterator begin_reinterpret(UInt new_size) const;
- /// Get a const_scalar_iterator on the end of the Array considered of shape (new_size)
+ /// Get a const_scalar_iterator on the end of the Array considered of shape
+ /// (new_size)
inline const_scalar_iterator end_reinterpret(UInt new_size) const;
/* ------------------------------------------------------------------------ */
/// Get a vector_iterator on the beginning of the Array
inline vector_iterator begin(UInt n);
/// Get a vector_iterator on the end of the Array
inline vector_iterator end(UInt n);
/// Get a vector_iterator on the beginning of the Array
inline const_vector_iterator begin(UInt n) const;
/// Get a vector_iterator on the end of the Array
inline const_vector_iterator end(UInt n) const;
/// Get a vector_iterator on the begining of the Array considered of shape
/// (new_size, n)
inline vector_iterator begin_reinterpret(UInt n, UInt new_size);
/// Get a vector_iterator on the end of the Array considered of shape
/// (new_size, n)
inline vector_iterator end_reinterpret(UInt n, UInt new_size);
/// Get a const_vector_iterator on the begining of the Array considered of
/// shape (new_size, n)
inline const_vector_iterator begin_reinterpret(UInt n, UInt new_size) const;
/// Get a const_vector_iterator on the end of the Array considered of shape
/// (new_size, n)
inline const_vector_iterator end_reinterpret(UInt n, UInt new_size) const;
/* ------------------------------------------------------------------------ */
/// Get a matrix_iterator on the begining of the Array (Matrices of size (m,
/// n))
inline matrix_iterator begin(UInt m, UInt n);
/// Get a matrix_iterator on the end of the Array (Matrices of size (m, n))
inline matrix_iterator end(UInt m, UInt n);
/// Get a const_matrix_iterator on the begining of the Array (Matrices of size
/// (m, n))
inline const_matrix_iterator begin(UInt m, UInt n) const;
/// Get a const_matrix_iterator on the end of the Array (Matrices of size (m,
/// n))
inline const_matrix_iterator end(UInt m, UInt n) const;
/// Get a matrix_iterator on the begining of the Array considered of shape
/// (new_size, m*n)
inline matrix_iterator begin_reinterpret(UInt m, UInt n, UInt size);
/// Get a matrix_iterator on the end of the Array considered of shape
/// (new_size, m*n)
inline matrix_iterator end_reinterpret(UInt m, UInt n, UInt size);
/// Get a const_matrix_iterator on the begining of the Array considered of
/// shape (new_size, m*n)
inline const_matrix_iterator begin_reinterpret(UInt m, UInt n,
UInt size) const;
/// Get a const_matrix_iterator on the end of the Array considered of shape
/// (new_size, m*n)
inline const_matrix_iterator end_reinterpret(UInt m, UInt n, UInt size) const;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// append a tuple of size nb_component containing value
inline void push_back(const_reference value);
/// append a vector
- //inline void push_back(const value_type new_elem[]);
+ // inline void push_back(const value_type new_elem[]);
/// append a Vector or a Matrix
template <template <typename> class C>
inline void push_back(const C<T> & new_elem);
/// append the value of the iterator
template <typename Ret> inline void push_back(const iterator<Ret> & it);
/// erase the value at position i
inline void erase(UInt i);
/// ask Nico, clarify
template <typename R> inline iterator<R> erase(const iterator<R> & it);
/// change the size of the Array
virtual void resize(UInt size);
/// change the size of the Array and initialize the values
virtual void resize(UInt size, const T & val);
/// 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
void extendComponentsInterlaced(UInt multiplicator, UInt stride);
/// search elem in the vector, return the position of the first occurrence or
/// -1 if not found
Int find(const_reference elem)
const; /// @see Array::find(const_reference elem) const
Int find(T elem[]) const;
/// @see Array::find(const_reference elem) const
template <template <typename> class C> inline Int find(const C<T> & elem);
/// set all entries of the array to 0
inline void clear() { std::fill_n(values, size * nb_component, T()); }
/// set all entries of the array to the value t
/// @param t value to fill the array with
inline void set(T t) { std::fill_n(values, size * nb_component, t); }
/// set all tuples of the array to a given vector or matrix
/// @param vm Matrix or Vector to fill the array with
template <template <typename> class C> inline void set(const C<T> & vm);
/// Append the content of the other array to the current one
void append(const Array<T> & other);
/// copy another Array in the current Array, the no_sanity_check allows you to
/// force the copy in cases where you know what you do with two non matching
/// Arrays in terms of n
void copy(const Array<T, is_scal> & other, bool no_sanity_check = false);
/// give the address of the memory allocated for this vector
T * storage() const { return values; };
/// function to print the containt of the class
- virtual void printself(std::ostream & stream, int indent = 0) const;
+ void printself(std::ostream & stream, int indent = 0) const override;
protected:
/// perform the allocation for the constructors
void allocate(UInt size, UInt nb_component = 1);
/// resize initializing with uninitialized_fill if fill is set
void resizeUnitialized(UInt new_size, bool fill, const T & val = T());
/* ------------------------------------------------------------------------ */
/* Operators */
/* ------------------------------------------------------------------------ */
public:
/// substraction entry-wise
Array<T, is_scal> & operator-=(const Array<T, is_scal> & other);
/// addition entry-wise
Array<T, is_scal> & operator+=(const Array<T, is_scal> & other);
/// multiply evry entry by alpha
Array<T, is_scal> & operator*=(const T & alpha);
/// check if the array are identical entry-wise
bool operator==(const Array<T, is_scal> & other) const;
/// @see Array::operator==(const Array<T, is_scal> & other) const
bool operator!=(const Array<T, is_scal> & other) const;
/// return a reference to the j-th entry of the i-th tuple
inline reference operator()(UInt i, UInt j = 0);
/// return a const reference to the j-th entry of the i-th tuple
inline const_reference operator()(UInt i, UInt j = 0) const;
/// return a reference to the ith component of the 1D array
inline reference operator[](UInt i);
/// return a const reference to the ith component of the 1D array
inline const_reference operator[](UInt i) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// array of values
T * values; // /!\ very dangerous
};
-#include "aka_array_tmpl.hh"
-
-__END_AKANTU__
-
-#include "aka_types.hh"
-
-__BEGIN_AKANTU__
-
/* -------------------------------------------------------------------------- */
/* Inline Functions Array<T, is_scal> */
/* -------------------------------------------------------------------------- */
template <typename T, bool is_scal>
inline std::ostream & operator<<(std::ostream & stream,
const Array<T, is_scal> & _this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
/* Inline Functions ArrayBase */
/* -------------------------------------------------------------------------- */
inline std::ostream & operator<<(std::ostream & stream,
const ArrayBase & _this) {
_this.printself(stream);
return stream;
}
-__END_AKANTU__
+} // namespace akantu
+
+#include "aka_array_tmpl.hh"
+#include "aka_types.hh"
#endif /* __AKANTU_VECTOR_HH__ */
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index ddbe3cc30..364cdd368 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1535 +1,1542 @@
/**
* @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> */
/* -------------------------------------------------------------------------- */
-
-__END_AKANTU__
-
+#include "aka_array.hh"
+/* -------------------------------------------------------------------------- */
#include <memory>
+/* -------------------------------------------------------------------------- */
+
+#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__
+#define __AKANTU_AKA_ARRAY_TMPL_HH__
-__BEGIN_AKANTU__
+namespace akantu {
/* -------------------------------------------------------------------------- */
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);
// }
/* -------------------------------------------------------------------------- */
/**
* 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>
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);
}
/* -------------------------------------------------------------------------- */
/**
* 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);
}
/* -------------------------------------------------------------------------- */
/**
* 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>
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);
}
}
/* -------------------------------------------------------------------------- */
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.getSize(), false);
T * tmp = values + nb_component * old_size;
std::uninitialized_copy(
other.storage(), other.storage() + other.getSize() * 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(NULL) {
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(NULL) {
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 = NULL;
} else {
values = static_cast<T *>(malloc(nb_component * size * sizeof(T)));
AKANTU_DEBUG_ASSERT(values != NULL,
"Cannot allocate "
<< printMemorySize<T>(size * nb_component) << " ("
<< id << ")");
}
if (values == NULL) {
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();
}
/* -------------------------------------------------------------------------- */
/**
* 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 = NULL;
} else {
- T * tmp_ptr = static_cast<T *>(
+ auto * tmp_ptr = static_cast<T *>(
realloc(values, new_size * nb_component * sizeof(T)));
if (tmp_ptr == NULL) {
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;
- T * tmp_ptr = static_cast<T *>(
+ auto * tmp_ptr = static_cast<T *>(
realloc(values, size_to_alloc * nb_component * sizeof(T)));
AKANTU_DEBUG_ASSERT(
tmp_ptr != NULL,
"Cannot allocate " << printMemorySize<T>(size_to_alloc * nb_component));
if (tmp_ptr == NULL) {
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>
Int Array<T, is_scal>::find(const T & elem) const {
AKANTU_DEBUG_IN();
UInt i = 0;
for (; (i < size) && (values[i] != elem); ++i)
;
AKANTU_DEBUG_OUT();
return (i == size) ? -1 : (Int)i;
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> Int 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 -1;
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <template <typename> class C>
inline Int Array<T, is_scal>::find(const C<T> & elem) {
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.getSize(); ++i) {
stream << "{";
for (UInt j = 0; j < vect.getNbComponent(); ++j) {
stream << vect(i, j);
if (j != vect.getNbComponent() - 1)
stream << ", ";
}
stream << "}";
if (i != vect.getSize() - 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; }
/* -------------------------------------------------------------------------- */
/* Iterators */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <class R, class IR, bool is_r_scal>
class Array<T, is_scal>::iterator_internal {
public:
- typedef R value_type;
- typedef R * pointer;
- typedef R & reference;
- typedef typename R::proxy proxy;
- typedef const typename R::proxy const_proxy;
- typedef const R & const_reference;
- typedef IR internal_value_type;
- typedef IR * internal_pointer;
- typedef std::ptrdiff_t difference_type;
- typedef std::random_access_iterator_tag iterator_category;
+ 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() : _offset(0), initial(NULL), ret(NULL), ret_ptr(NULL){};
+ iterator_internal() : initial(NULL), ret(NULL), ret_ptr(NULL){};
iterator_internal(pointer_type data, UInt _offset)
: _offset(_offset), initial(data), ret(NULL), ret_ptr(data) {
AKANTU_DEBUG_ERROR(
"The constructor should never be called it is just an ugly trick...");
}
iterator_internal(pointer wrapped)
: _offset(wrapped->size()), initial(wrapped->storage()),
ret(const_cast<internal_pointer>(wrapped)),
ret_ptr(wrapped->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 = new internal_value_type(*it.ret, false);
}
}
virtual ~iterator_internal() { delete ret; };
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 = new 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;
};
inline iterator_internal & operator++() {
ret_ptr += _offset;
return *this;
};
inline iterator_internal & operator--() {
ret_ptr -= _offset;
return *this;
};
inline iterator_internal & operator+=(const UInt n) {
ret_ptr += _offset * n;
return *this;
}
inline iterator_internal & operator-=(const UInt n) {
ret_ptr -= _offset * n;
return *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 iterator_internal operator+(difference_type n) {
iterator_internal tmp(*this);
tmp += n;
return tmp;
}
inline iterator_internal operator-(difference_type n) {
iterator_internal tmp(*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;
+ UInt _offset{0};
pointer_type initial;
internal_pointer ret;
pointer_type ret_ptr;
};
/* -------------------------------------------------------------------------- */
/**
* Specialization for scalar types
*/
template <class T, bool is_scal>
template <class R, class IR>
class Array<T, is_scal>::iterator_internal<R, IR, true> {
public:
- typedef R value_type;
- typedef R * pointer;
- typedef R & reference;
- typedef const R & const_reference;
- typedef IR internal_value_type;
- typedef IR * internal_pointer;
- typedef std::ptrdiff_t difference_type;
- typedef std::random_access_iterator_tag iterator_category;
+ 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 = NULL,
__attribute__((unused)) UInt _offset = 1)
: _offset(_offset), ret(data), initial(data){};
iterator_internal(const iterator_internal & it) {
if (this != &it) {
this->ret = it.ret;
this->initial = it.initial;
}
}
- virtual ~iterator_internal(){};
+ virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) {
if (this != &it) {
this->ret = it.ret;
this->initial = it.initial;
}
return *this;
}
UInt getCurrentIndex() {
return (this->ret - this->initial) / this->_offset;
};
inline reference operator*() { return *ret; };
inline const_reference operator*() const { return *ret; };
inline pointer operator->() { return ret; };
inline iterator_internal & operator++() {
++ret;
return *this;
};
inline iterator_internal & operator--() {
--ret;
return *this;
};
inline iterator_internal & operator+=(const UInt n) {
ret += n;
return *this;
}
inline iterator_internal & operator-=(const UInt n) {
ret -= n;
return *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 iterator_internal operator-(difference_type n) {
return iterator_internal(ret - n);
}
inline iterator_internal operator+(difference_type n) {
return iterator_internal(ret + n);
}
inline difference_type operator-(const iterator_internal & b) {
return ret - b.ret;
}
inline pointer data() const { return ret; }
inline difference_type offset() const { return _offset; }
protected:
difference_type _offset;
pointer ret;
pointer initial;
};
/* -------------------------------------------------------------------------- */
/* Begin/End functions implementation */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Vector<T> * to the
* first tuple of the array.
* @param n vector size. Has to be equal to nb_component. This unfortunate
* redundancy is necessary to distinguish it from ::begin() which it
* overloads. If compiled in debug mode, an incorrect value of n will result
* in an exception being thrown. Optimized code will fail in an unpredicted
* manner.
* @return a vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::vector_iterator
Array<T, is_scal>::begin(UInt n) {
AKANTU_DEBUG_ASSERT(nb_component == n,
"The iterator is not compatible with the type Array("
<< n << ")");
return vector_iterator(new Vector<T>(values, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Vector<T> * pointing
* *past* the last tuple of the array.
* @param n vector size. see Array::begin(UInt n) for more
* @return a vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::vector_iterator
Array<T, is_scal>::end(UInt n) {
AKANTU_DEBUG_ASSERT(nb_component == n,
"The iterator is not compatible with the type Array("
<< n << ")");
return vector_iterator(new Vector<T>(values + nb_component * size, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Vector<T> * to the
* first tuple of the array.
* @param n vector size. see Array::begin(UInt n) for more
* @return a vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_vector_iterator
Array<T, is_scal>::begin(UInt n) const {
AKANTU_DEBUG_ASSERT(nb_component == n,
"The iterator is not compatible with the type Array("
<< n << ")");
return const_vector_iterator(new Vector<T>(values, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Vector<T> * pointing
* *past* the last tuple of the array.
* @param n vector size. see Array::begin(UInt n) for more
* @return a const_vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_vector_iterator
Array<T, is_scal>::end(UInt n) const {
AKANTU_DEBUG_ASSERT(nb_component == n,
"The iterator is not compatible with the type Array("
<< n << ")");
return const_vector_iterator(new Vector<T>(values + nb_component * size, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Vector<T> * to the
* first tuple of the array.
*
* The reinterpret iterators allow to iterate over an array in any way that
* preserves the number of entries of the array. This can for instance be use
* full if the shape of the data in an array is not initially known.
* @param n vector size.
* @param size number of tuples in array. n times size must match the number
* of entries of the array. If compiled in debug mode, an incorrect
* combination of n and size will result
* in an exception being thrown. Optimized code will fail in an unpredicted
* manner.
* @return a vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::vector_iterator
Array<T, is_scal>::begin_reinterpret(UInt n,
__attribute__((unused)) UInt size) {
AKANTU_DEBUG_ASSERT(n * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << n
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return vector_iterator(new Vector<T>(values, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Vector<T> * pointing
* *past* the last tuple of the array.
* @param n vector size.
* @param size number of tuples in array. See Array::begin_reinterpret(UInt n,
* UInt size)
* @return a vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::vector_iterator
Array<T, is_scal>::end_reinterpret(UInt n, UInt size) {
AKANTU_DEBUG_ASSERT(n * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << n
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return vector_iterator(new Vector<T>(values + n * size, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Vector<T> * to the
* first tuple of the array.
* @param n vector size.
* @param size number of tuples in array. See Array::begin_reinterpret(UInt n,
* UInt size)
* @return a const_vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_vector_iterator
Array<T, is_scal>::begin_reinterpret(UInt n,
__attribute__((unused)) UInt size) const {
AKANTU_DEBUG_ASSERT(n * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << n
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return const_vector_iterator(new Vector<T>(values, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Vector<T> * pointing
* *past* the last tuple of the array.
* @param n vector size.
* @param size number of tuples in array. See Array::begin_reinterpret(UInt n,
* UInt size)
* @return a const_vector_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_vector_iterator
Array<T, is_scal>::end_reinterpret(UInt n, UInt size) const {
AKANTU_DEBUG_ASSERT(n * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << n
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return const_vector_iterator(new Vector<T>(values + n * size, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Matrix<T> * to the
* first tuple of the array.
* @param m number of rows
* @param n number of columns. m times n has to equal nb_component.
* If compiled in debug mode, an incorrect combination of m and n will result
* in an exception being thrown. Optimized code will fail in an unpredicted
* manner.
* @return a matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::matrix_iterator
Array<T, is_scal>::begin(UInt m, UInt n) {
AKANTU_DEBUG_ASSERT(nb_component == n * m,
"The iterator is not compatible with the type Matrix("
<< m << "," << n << ")");
return matrix_iterator(new Matrix<T>(values, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Matrix<T> * pointing
* *past* the last tuple of the array.
* @param m number of rows
* @param n number of columns. See Array::begin(UInt m, UInt n)
* @return a matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::matrix_iterator
Array<T, is_scal>::end(UInt m, UInt n) {
AKANTU_DEBUG_ASSERT(nb_component == n * m,
"The iterator is not compatible with the type Matrix("
<< m << "," << n << ")");
return matrix_iterator(new Matrix<T>(values + nb_component * size, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Matrix<T> * to the
* first tuple of the array.
* @param m number of rows
* @param n number of columns. See Array::begin(UInt m, UInt n)
* @return a matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_matrix_iterator
Array<T, is_scal>::begin(UInt m, UInt n) const {
AKANTU_DEBUG_ASSERT(nb_component == n * m,
"The iterator is not compatible with the type Matrix("
<< m << "," << n << ")");
return const_matrix_iterator(new Matrix<T>(values, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Matrix<T> * pointing
* *past* the last tuple of the array.
* @param m number of rows
* @param n number of columns. See Array::begin(UInt m, UInt n)
* @return a const_matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_matrix_iterator
Array<T, is_scal>::end(UInt m, UInt n) const {
AKANTU_DEBUG_ASSERT(nb_component == n * m,
"The iterator is not compatible with the type Matrix("
<< m << "," << n << ")");
return const_matrix_iterator(
new Matrix<T>(values + nb_component * size, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Matrix<T> * to the
* first tuple of the array.
*
* The reinterpret iterators allow to iterate over an array in any way that
* preserves the number of entries of the array. This can for instance be use
* full if the shape of the data in an array is not initially known.
* @param m number of rows
* @param n number of columns
* @param size number of tuples in array. m times n times size must match the
*number
* of entries of the array. If compiled in debug mode, an incorrect
* combination of m, n and size will result
* in an exception being thrown. Optimized code will fail in an unpredicted
* manner.
* @return a matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::matrix_iterator
Array<T, is_scal>::begin_reinterpret(UInt m, UInt n,
__attribute__((unused)) UInt size) {
AKANTU_DEBUG_ASSERT(n * m * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << m << "," << n
<< " = " << n * m
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return matrix_iterator(new Matrix<T>(values, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get an iterator that behaves like a pointer akantu::Matrix<T> * pointing
* *past* the last tuple of the array.
* @param m number of rows
* @param n number of columns
* @param size number of tuples in array. See Array::begin_reinterpret(UInt m,
* UInt n, UInt size)
* @return a matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::matrix_iterator
Array<T, is_scal>::end_reinterpret(UInt m, UInt n, UInt size) {
AKANTU_DEBUG_ASSERT(n * m * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << m << "," << n
<< " = " << n * m
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return matrix_iterator(new Matrix<T>(values + n * m * size, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Matrix<T> * to the
* first tuple of the array.
* @param m number of rows
* @param n number of columns
* @param size number of tuples in array. See Array::begin_reinterpret(UInt m,
* UInt n, UInt size)
* @return a const_matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_matrix_iterator
Array<T, is_scal>::begin_reinterpret(UInt m, UInt n,
__attribute__((unused)) UInt size) const {
AKANTU_DEBUG_ASSERT(n * m * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << m << "," << n
<< " = " << n * m
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return const_matrix_iterator(new Matrix<T>(values, m, n));
}
/* -------------------------------------------------------------------------- */
/**
* Get a const iterator that behaves like a pointer akantu::Matrix<T> * pointing
* *past* the last tuple of the array.
* @param m number of rows
* @param n number of columns
* @param size number of tuples in array. See Array::begin_reinterpret(UInt m,
* UInt n, UInt size)
* @return a const_matrix_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_matrix_iterator
Array<T, is_scal>::end_reinterpret(UInt m, UInt n, UInt size) const {
AKANTU_DEBUG_ASSERT(n * m * size == this->nb_component * this->size,
"The new values for size ("
<< size << ") and nb_component (" << m << "," << n
<< " = " << n * m
<< ") are not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return const_matrix_iterator(new Matrix<T>(values + n * m * size, m, n));
}
/* -------------------------------------------------------------------------- */
/** Get an iterator that behaves like a pointer T * to the
* first entry in the member array values
* @return a scalar_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::scalar_iterator Array<T, is_scal>::begin() {
AKANTU_DEBUG_ASSERT(
nb_component == 1,
"this iterator cannot be used on a vector which has nb_component != 1");
return scalar_iterator(values);
}
/* -------------------------------------------------------------------------- */
/*! Get an iterator that behaves like a pointer T * that points *past* the
* last entry in the member array values
* @return a scalar_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::scalar_iterator Array<T, is_scal>::end() {
AKANTU_DEBUG_ASSERT(
nb_component == 1,
"this iterator cannot be used on a array which has nb_component != 1");
return scalar_iterator(values + size);
}
/* -------------------------------------------------------------------------- */
/*! Get a const iterator that behaves like a pointer T * to the
* first entry in the member array values
* @return a const_scalar_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_scalar_iterator
Array<T, is_scal>::begin() const {
AKANTU_DEBUG_ASSERT(
nb_component == 1,
"this iterator cannot be used on a array which has nb_component != 1");
return const_scalar_iterator(values);
}
/* -------------------------------------------------------------------------- */
/*! Get a const iterator that behaves like a pointer T * that points *past* the
* last entry in the member array values
* @return a const_scalar_iterator
*/
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_scalar_iterator
Array<T, is_scal>::end() const {
AKANTU_DEBUG_ASSERT(
nb_component == 1,
"this iterator cannot be used on a array which has nb_component != 1");
return const_scalar_iterator(values + size);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline typename Array<T, is_scal>::scalar_iterator
-Array<T, is_scal>::begin_reinterpret(UInt new_size) {
+Array<T, is_scal>::begin_reinterpret(__attribute__((unused)) UInt new_size) {
AKANTU_DEBUG_ASSERT(new_size == this->nb_component * this->size,
"The new values for size ("
<< new_size
<< ") is not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return scalar_iterator(values);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline typename Array<T, is_scal>::scalar_iterator
Array<T, is_scal>::end_reinterpret(UInt new_size) {
AKANTU_DEBUG_ASSERT(new_size == this->nb_component * this->size,
"The new values for size ("
<< new_size
<< ") is not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
- return scalar_iterator(values + size);
+ return scalar_iterator(values + new_size);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_scalar_iterator
-Array<T, is_scal>::begin_reinterpret(UInt new_size) const {
+Array<T, is_scal>::begin_reinterpret(__attribute__((unused))
+ UInt new_size) const {
AKANTU_DEBUG_ASSERT(new_size == this->nb_component * this->size,
"The new values for size ("
<< new_size
<< ") is not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
return const_scalar_iterator(values);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline typename Array<T, is_scal>::const_scalar_iterator
Array<T, is_scal>::end_reinterpret(UInt new_size) const {
AKANTU_DEBUG_ASSERT(new_size == this->nb_component * this->size,
"The new values for size ("
<< new_size
<< ") is not compatible with the one of this array("
<< this->size << "," << this->nb_component << ")");
- return const_scalar_iterator(values + size);
+ return const_scalar_iterator(values + new_size);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::const_iterator : public iterator_internal<const R, R> {
public:
typedef iterator_internal<const R, R> parent;
- typedef typename parent::value_type value_type;
- typedef typename parent::pointer pointer;
- typedef typename parent::reference reference;
- typedef typename parent::difference_type difference_type;
- typedef typename parent::iterator_category iterator_category;
+ 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<R> & it) : parent(it) {}
inline const_iterator operator+(difference_type n) {
return parent::operator+(n);
}
inline const_iterator operator-(difference_type n) {
return parent::operator-(n);
}
inline difference_type operator-(const const_iterator & b) {
return parent::operator-(b);
}
inline const_iterator & operator++() {
parent::operator++();
return *this;
};
inline const_iterator & operator--() {
parent::operator--();
return *this;
};
inline const_iterator & operator+=(const UInt n) {
parent::operator+=(n);
return *this;
}
};
// #endif
// #if defined(AKANTU_CORE_CXX11)
// template<class R> using iterator = iterator_internal<R>;
// #else
template <class T, class R, bool issame = is_same<T, R>::value>
struct ConstConverterIteratorHelper {
- typedef typename Array<T>::template const_iterator<R> const_iterator;
- typedef typename Array<T>::template iterator<R> iterator;
+ 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(new R(*it, false));
}
};
template <class T, class R> struct ConstConverterIteratorHelper<T, R, true> {
- typedef typename Array<T>::template const_iterator<R> const_iterator;
- typedef typename Array<T>::template iterator<R> iterator;
+ 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(), it.offset());
}
};
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::iterator : public iterator_internal<R> {
public:
- typedef iterator_internal<R> parent;
- typedef typename parent::value_type value_type;
- typedef typename parent::pointer pointer;
- typedef typename parent::reference reference;
- typedef typename parent::difference_type difference_type;
- typedef typename parent::iterator_category iterator_category;
+ using parent = iterator_internal<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:
iterator() : parent(){};
iterator(pointer_type data, UInt offset) : parent(data, offset){};
iterator(pointer warped) : parent(warped) {}
iterator(const parent & it) : parent(it) {}
// iterator(const iterator<R> & it) : parent(it) {}
operator const_iterator<R>() {
return ConstConverterIteratorHelper<T, R>::convert(*this);
}
inline iterator operator+(difference_type n) {
return parent::operator+(n);
;
}
inline iterator operator-(difference_type n) {
return parent::operator-(n);
;
}
inline difference_type operator-(const iterator & b) {
return parent::operator-(b);
}
inline iterator & operator++() {
parent::operator++();
return *this;
};
inline iterator & operator--() {
parent::operator--();
return *this;
};
inline iterator & operator+=(const UInt n) {
parent::operator+=(n);
return *this;
}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
-inline Array<T, is_scal>::iterator<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
+
+} // akantu
+
+#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */
diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh
index 1e2cf4d15..e4a818674 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,447 +1,447 @@
/**
* @file aka_common.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Thu Jan 21 2016
*
* @brief common type descriptions for 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/>.
*
* @section DESCRIPTION
*
* All common things to be included in the projects files
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMON_HH__
#define __AKANTU_COMMON_HH__
/* -------------------------------------------------------------------------- */
#include <list>
#include <limits>
#include <type_traits>
#if __cplusplus < 201402L
namespace std {
template< bool B, class T = void >
using enable_if_t = typename enable_if<B,T>::type;
}
#endif
/* -------------------------------------------------------------------------- */
#define __BEGIN_AKANTU__ namespace akantu {
#define __END_AKANTU__ }
/* -------------------------------------------------------------------------- */
#define __BEGIN_AKANTU_DUMPER__ namespace dumper {
#define __END_AKANTU_DUMPER__ }
/* -------------------------------------------------------------------------- */
#if defined(WIN32)
#define __attribute__(x)
#endif
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
#include "aka_error.hh"
#include "aka_safe_enum.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Common types */
/* -------------------------------------------------------------------------- */
-typedef std::string ID;
+using ID = std::string;
#ifdef AKANTU_NDEBUG
static const Real REAL_INIT_VALUE = Real(0.);
#else
static const Real REAL_INIT_VALUE = std::numeric_limits<Real>::quiet_NaN();
#endif
/* -------------------------------------------------------------------------- */
/* Memory types */
/* -------------------------------------------------------------------------- */
-typedef UInt MemoryID;
+using MemoryID = UInt;
-typedef std::string Surface;
+using Surface = std::string;
typedef std::pair<Surface, Surface> SurfacePair;
-typedef std::list<SurfacePair> SurfacePairList;
+using SurfacePairList = std::list<SurfacePair>;
/* -------------------------------------------------------------------------- */
extern const UInt _all_dimensions;
/* -------------------------------------------------------------------------- */
/* Mesh/FEM/Model types */
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "aka_element_classes_info.hh"
__BEGIN_AKANTU__
/// small help to use names for directions
enum SpacialDirection { _x = 0, _y = 1, _z = 2 };
/// enum MeshIOType type of mesh reader/writer
enum MeshIOType {
_miot_auto, ///< Auto guess of the reader to use based on the extension
_miot_gmsh, ///< Gmsh files
_miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has
/// structures elements
_miot_diana, ///< TNO Diana mesh format
_miot_abaqus ///< Abaqus mesh format
};
/// enum AnalysisMethod type of solving method used to solve the equation of
/// motion
enum AnalysisMethod {
_static = 0,
_implicit_dynamic = 1,
_explicit_lumped_mass = 2,
_explicit_lumped_capacity = 2,
_explicit_consistent_mass = 3
};
/// enum DOFSupportType defines which kind of dof that can exists
enum DOFSupportType { _dst_nodal, _dst_generic };
/// Type of non linear resolution available in akantu
enum NonLinearSolverType {
_nls_linear, ///< No non linear convergence loop
_nls_newton_raphson, ///< Regular Newton-Raphson
_nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent
_nls_lumped, ///< Case of lumped mass or equivalent matrix
_nls_auto ///< This will take a default value that make sense in case of
/// model::getNewSolver
};
/// Define the node/dof type
enum NodeType : Int {
_nt_pure_gost = -3,
_nt_master = -2,
_nt_normal = -1
};
/// Type of time stepping solver
enum TimeStepSolverType {
_tsst_static, ///< Static solution
_tsst_dynamic, ///< Dynamic solver
_tsst_dynamic_lumped, ///< Dynamic solver with lumped mass
};
/// Type of integration scheme
enum IntegrationSchemeType {
_ist_pseudo_time, ///< Pseudo Time
_ist_forward_euler, ///< GeneralizedTrapezoidal(0)
_ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2)
_ist_backward_euler, ///< GeneralizedTrapezoidal(1)
_ist_central_difference, ///< NewmarkBeta(0, 1/2)
_ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2)
_ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2)
_ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2)
_ist_newmark_beta, ///< generic NewmarkBeta with user defined
/// alpha and beta
_ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user
/// defined alpha
};
/// enum SolveConvergenceCriteria different convergence criteria
enum SolveConvergenceCriteria {
_scc_residual, ///< Use residual to test the convergence
_scc_solution, ///< Use solution to test the convergence
_scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb
};
/// enum CohesiveMethod type of insertion of cohesive elements
enum CohesiveMethod { _intrinsic, _extrinsic };
/// @enum SparseMatrixType type of sparse matrix used
enum MatrixType { _unsymmetric, _symmetric };
/* -------------------------------------------------------------------------- */
/* Ghosts handling */
/* -------------------------------------------------------------------------- */
-typedef ID SynchronizerID;
+using SynchronizerID = ID;
/// @enum CommunicatorType type of communication method to use
enum CommunicatorType { _communicator_mpi, _communicator_dummy };
/// @enum SynchronizationTag type of synchronizations
enum SynchronizationTag {
//--- Generic tags ---
_gst_whatever,
_gst_update,
_gst_size,
//--- SolidMechanicsModel tags ---
_gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass
_gst_smm_for_gradu, ///< synchronization of the
/// SolidMechanicsModel.displacement
_gst_smm_boundary, ///< synchronization of the boundary, forces, velocities
/// and displacement
_gst_smm_uv, ///< synchronization of the nodal velocities and displacement
_gst_smm_res, ///< synchronization of the nodal residual
_gst_smm_init_mat, ///< synchronization of the data to initialize materials
_gst_smm_stress, ///< synchronization of the stresses to compute the internal
/// forces
_gst_smmc_facets, ///< synchronization of facet data to setup facet synch
_gst_smmc_facets_conn, ///< synchronization of facet global connectivity
_gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet
/// synch
_gst_smmc_damage, ///< synchronization of damage
// --- GlobalIdsUpdater tags ---
_gst_giu_global_conn, ///< synchronization of global connectivities
// --- CohesiveElementInserter tags ---
_gst_ce_groups, ///< synchronization of cohesive element insertion depending
/// on facet groups
// --- GroupManager tags ---
_gst_gm_clusters, ///< synchronization of clusters
// --- HeatTransfer tags ---
_gst_htm_capacity, ///< synchronization of the nodal heat capacity
_gst_htm_temperature, ///< synchronization of the nodal temperature
_gst_htm_gradient_temperature, ///< synchronization of the element gradient
/// temperature
// --- LevelSet tags ---
_gst_htm_phi, ///< synchronization of the nodal level set value phi
_gst_htm_gradient_phi, ///< synchronization of the element gradient phi
//--- Material non local ---
_gst_mnl_for_average, ///< synchronization of data to average in non local
/// material
_gst_mnl_weight, ///< synchronization of data for the weight computations
// --- NeighborhoodSynchronization tags ---
_gst_nh_criterion,
// --- General tags ---
_gst_test, ///< Test tag
_gst_user_1, ///< tag for user simulations
_gst_user_2, ///< tag for user simulations
_gst_material_id, ///< synchronization of the material ids
_gst_for_dump, ///< everything that needs to be synch before dump
// --- Contact & Friction ---
_gst_cf_nodal, ///< synchronization of disp, velo, and current position
_gst_cf_incr, ///< synchronization of increment
// --- Solver tags ---
_gst_solver_solution ///< synchronization of the solution obained with the
/// PETSc solver
};
/// standard output stream operator for SynchronizationTag
inline std::ostream & operator<<(std::ostream & stream,
SynchronizationTag type);
/// @enum GhostType type of ghost
enum GhostType {
_not_ghost,
_ghost,
_casper // not used but a real cute ghost
};
/* -------------------------------------------------------------------------- */
struct GhostType_def {
- typedef GhostType type;
+ using type = GhostType;
static const type _begin_ = _not_ghost;
static const type _end_ = _casper;
};
-typedef safe_enum<GhostType_def> ghost_type_t;
+using ghost_type_t = safe_enum<GhostType_def>;
extern ghost_type_t ghost_types;
/// standard output stream operator for GhostType
inline std::ostream & operator<<(std::ostream & stream, GhostType type);
/// @enum SynchronizerOperation reduce operation that the synchronizer can
/// perform
enum SynchronizerOperation {
_so_sum,
_so_min,
_so_max,
_so_prod,
_so_land,
_so_band,
_so_lor,
_so_bor,
_so_lxor,
_so_bxor,
_so_min_loc,
_so_max_loc,
_so_null
};
/* -------------------------------------------------------------------------- */
/* Global defines */
/* -------------------------------------------------------------------------- */
#define AKANTU_MIN_ALLOCATION 2000
#define AKANTU_INDENT " "
#define AKANTU_INCLUDE_INLINE_IMPL
/* -------------------------------------------------------------------------- */
template <class T> struct is_scalar {
enum { value = false };
};
#define AKANTU_SPECIFY_IS_SCALAR(type) \
template <> struct is_scalar<type> { \
enum { value = true }; \
}
AKANTU_SPECIFY_IS_SCALAR(Real);
AKANTU_SPECIFY_IS_SCALAR(UInt);
AKANTU_SPECIFY_IS_SCALAR(Int);
AKANTU_SPECIFY_IS_SCALAR(bool);
template <typename T1, typename T2> struct is_same {
enum { value = false }; // is_same represents a bool.
};
template <typename T> struct is_same<T, T> {
enum { value = true };
};
/* -------------------------------------------------------------------------- */
#define AKANTU_SET_MACRO(name, variable, type) \
inline void set##name(type variable) { this->variable = variable; }
#define AKANTU_GET_MACRO(name, variable, type) \
inline type get##name() const { return variable; }
#define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \
inline type get##name() { return variable; }
#define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \
inline con Array<type> & get##name( \
const support & el_type, const GhostType & ghost_type = _not_ghost) \
con { \
return variable(el_type, ghost_type); \
}
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, )
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const)
#define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, )
#define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const)
/* -------------------------------------------------------------------------- */
/// initialize the static part of akantu
void initialize(int & argc, char **& argv);
/// initialize the static part of akantu and read the global input_file
void initialize(const std::string & input_file, int & argc, char **& argv);
/* -------------------------------------------------------------------------- */
/// finilize correctly akantu and clean the memory
void finalize();
/* -------------------------------------------------------------------------- */
/// Read an new input file
void readInputFile(const std::string & input_file);
/* -------------------------------------------------------------------------- */
/*
* For intel compiler annoying remark
*/
// #if defined(__INTEL_COMPILER)
// /// remark #981: operands are evaluated in unspecified order
// #pragma warning(disable : 981)
// /// remark #383: value copied to temporary, reference to temporary used
// #pragma warning(disable : 383)
// #endif // defined(__INTEL_COMPILER)
/* -------------------------------------------------------------------------- */
/* string manipulation */
/* -------------------------------------------------------------------------- */
inline std::string to_lower(const std::string & str);
/* -------------------------------------------------------------------------- */
inline std::string trim(const std::string & to_trim);
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/// give a string representation of the a human readable size in bit
template <typename T> std::string printMemorySize(UInt size);
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "aka_fwd.hh"
__BEGIN_AKANTU__
/// get access to the internal argument parser
cppargparse::ArgumentParser & getStaticArgumentParser();
/// get access to the internal input file parser
Parser & getStaticParser();
/// get access to the user part of the internal input file parser
const ParserSection & getUserParser();
__END_AKANTU__
#include "aka_common_inline_impl.cc"
/* -------------------------------------------------------------------------- */
#if AKANTU_INTEGER_SIZE == 4
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9
#elif AKANTU_INTEGER_SIZE == 8
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL
#endif
namespace std {
/**
* Hashing function for pairs based on hash_combine from boost The magic number
* is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f]
* @f[\frac{2^32}{\phi} = 0x9e3779b9@f]
* http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine
* http://burtleburtle.net/bob/hash/doobs.html
*/
template <typename a, typename b> struct hash<std::pair<a, b> > {
public:
hash() : ah(), bh() {}
size_t operator()(const std::pair<a, b> & p) const {
size_t seed = ah(p.first);
return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) +
(seed >> 2);
}
private:
const hash<a> ah;
const hash<b> bh;
};
} //std
#endif /* __AKANTU_COMMON_HH__ */
diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in
index 505ec52cc..f5e77c1fd 100644
--- a/src/common/aka_config.hh.in
+++ b/src/common/aka_config.hh.in
@@ -1,115 +1,104 @@
/**
* @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 \
+ AKANTU_VERSION_PATCH)
@AKANTU_TYPES_EXTRA_INCLUDES@
namespace akantu {
- typedef @AKANTU_FLOAT_TYPE@ Real;
- typedef @AKANTU_SIGNED_INTEGER_TYPE@ Int;
- typedef @AKANTU_UNSIGNED_INTEGER_TYPE@ UInt;
-
- // template<class Key, class Ty>
- // struct unordered_map {
- // typedef typename @AKANTU_UNORDERED_MAP_TYPE@<Key, Ty> type;
- // };
-
- // template<class T>
- // std::size_t hash(const T & t) {
- // typedef typename @AKANTU_HASH_TYPE@<T> hash_type;
- // return hash_type()(t);
- // }
-}
+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_NLOPT
#cmakedefine AKANTU_USE_CPPARRAY
#cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY
#cmakedefine AKANTU_EXTRA_MATERIALS
#cmakedefine AKANTU_STUDENTS_EXTRA_PACKAGE
#cmakedefine AKANTU_DAMAGE_NON_LOCAL
#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/common/aka_element_classes_info.hh.in b/src/common/aka_element_classes_info.hh.in
index 1b4ba794b..e71f120fc 100644
--- a/src/common/aka_element_classes_info.hh.in
+++ b/src/common/aka_element_classes_info.hh.in
@@ -1,199 +1,199 @@
/**
* @file aka_element_classes_info.hh.in
*
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Sun Jul 19 2015
* @date last modification: Fri Oct 16 2015
*
* @brief Declaration of the enums for the element classes
*
* @section LICENSE
*
* Copyright (©) 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 <boost/preprocessor.hpp>
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Element Types */
/* -------------------------------------------------------------------------- */
/// @enum ElementType type of elements
enum ElementType {
_not_defined,
@AKANTU_ELEMENT_TYPES_ENUM@
_max_element_type
};
@AKANTU_ELEMENT_TYPES_BOOST_SEQ@
@AKANTU_ALL_ELEMENT_BOOST_SEQ@
/* -------------------------------------------------------------------------- */
/* Element Kinds */
/* -------------------------------------------------------------------------- */
@AKANTU_ELEMENT_KINDS_BOOST_SEQ@
@AKANTU_ELEMENT_KIND_BOOST_SEQ@
#ifndef SWIG
enum ElementKind {
BOOST_PP_SEQ_ENUM(AKANTU_ELEMENT_KIND),
_ek_not_defined
};
/* -------------------------------------------------------------------------- */
struct ElementKind_def {
- typedef ElementKind type;
+ using type = ElementKind;
static const type _begin_ = BOOST_PP_SEQ_HEAD(AKANTU_ELEMENT_KIND);
static const type _end_ = _ek_not_defined;
};
-typedef safe_enum<ElementKind_def> element_kind_t;
+using element_kind_t = safe_enum<ElementKind_def> ;
#else
enum ElementKind;
#endif
/* -------------------------------------------------------------------------- */
/// @enum GeometricalType type of element potentially contained in a Mesh
enum GeometricalType {
@AKANTU_GEOMETRICAL_TYPES_ENUM@
_gt_not_defined
};
/* -------------------------------------------------------------------------- */
/* Interpolation Types */
/* -------------------------------------------------------------------------- */
@AKANTU_INTERPOLATION_TYPES_BOOST_SEQ@
#ifndef SWIG
/// @enum InterpolationType type of elements
enum InterpolationType {
BOOST_PP_SEQ_ENUM(AKANTU_INTERPOLATION_TYPES),
_itp_not_defined
};
#else
enum InterpolationType;
#endif
/* -------------------------------------------------------------------------- */
/* Some sub types less probable to change */
/* -------------------------------------------------------------------------- */
/// @enum GeometricalShapeType types of shapes to define the contains
/// function in the element classes
enum GeometricalShapeType {
@AKANTU_GEOMETRICAL_SHAPES_ENUM@
_gst_not_defined
};
/* -------------------------------------------------------------------------- */
/// @enum GaussIntegrationType classes of types using common
/// description of the gauss point position and weights
enum GaussIntegrationType {
@AKANTU_GAUSS_INTEGRATION_TYPES_ENUM@
_git_not_defined
};
/* -------------------------------------------------------------------------- */
/// @enum InterpolationKind the family of interpolation types
enum InterpolationKind {
@AKANTU_INTERPOLATION_KIND_ENUM@
_itk_not_defined
};
/* -------------------------------------------------------------------------- */
// BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING
#define AKANTU_BOOST_CASE_MACRO(r,macro,_type) \
case _type : { macro(_type); break; }
#define AKANTU_BOOST_LIST_SWITCH(macro1, list1, var) \
do { \
switch(var) { \
BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \
default: { \
AKANTU_DEBUG_ERROR("Type (" << var << ") not handled by this function"); \
} \
} \
} while(0)
#define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1) \
AKANTU_BOOST_LIST_SWITCH(macro1, list1, type)
#define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro) \
AKANTU_BOOST_ELEMENT_SWITCH(macro, \
AKANTU_ALL_ELEMENT_TYPE)
#define AKANTU_BOOST_LIST_MACRO(r, macro, type) \
macro(type)
#define AKANTU_BOOST_APPLY_ON_LIST(macro, list) \
BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list)
#define AKANTU_BOOST_ALL_ELEMENT_LIST(macro) \
AKANTU_BOOST_APPLY_ON_LIST(macro, \
AKANTU_ALL_ELEMENT_TYPE)
#define AKANTU_GET_ELEMENT_LIST(kind) \
AKANTU##kind##_ELEMENT_TYPE
#define AKANTU_BOOST_KIND_ELEMENT_SWITCH(macro, kind) \
AKANTU_BOOST_ELEMENT_SWITCH(macro, \
AKANTU_GET_ELEMENT_LIST(kind))
// BOOST_PP_SEQ_TO_LIST does not exists in Boost < 1.49
#define AKANTU_GENERATE_KIND_LIST(seq) \
BOOST_PP_TUPLE_TO_LIST(BOOST_PP_SEQ_SIZE(seq), \
BOOST_PP_SEQ_TO_TUPLE(seq))
#define AKANTU_ELEMENT_KIND_BOOST_LIST AKANTU_GENERATE_KIND_LIST(AKANTU_ELEMENT_KIND)
#define AKANTU_BOOST_ALL_KIND_LIST(macro, list) \
BOOST_PP_LIST_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list)
#define AKANTU_BOOST_ALL_KIND(macro) \
AKANTU_BOOST_ALL_KIND_LIST(macro, AKANTU_ELEMENT_KIND_BOOST_LIST)
#define AKANTU_BOOST_ALL_KIND_SWITCH(macro) \
AKANTU_BOOST_LIST_SWITCH(macro, \
AKANTU_ELEMENT_KIND, \
kind)
@AKANTU_ELEMENT_KINDS_BOOST_MACROS@
// /// define kept for compatibility reasons (they are most probably not needed
// /// anymore) \todo check if they can be removed
// #define AKANTU_REGULAR_ELEMENT_TYPE AKANTU_ek_regular_ELEMENT_TYPE
// #define AKANTU_COHESIVE_ELEMENT_TYPE AKANTU_ek_cohesive_ELEMENT_TYPE
// #define AKANTU_STRUCTURAL_ELEMENT_TYPE AKANTU_ek_structural_ELEMENT_TYPE
// #define AKANTU_IGFEM_ELEMENT_TYPE AKANTU_ek_igfem_ELEMENT_TYPE
/* -------------------------------------------------------------------------- */
/* Lists of interests for FEEngineTemplate functions */
/* -------------------------------------------------------------------------- */
@AKANTU_FE_ENGINE_LISTS@
__END_AKANTU__
#include "aka_element_classes_info_inline_impl.cc"
diff --git a/src/common/aka_error.hh b/src/common/aka_error.hh
index eb667dc7c..d10e1a59e 100644
--- a/src/common/aka_error.hh
+++ b/src/common/aka_error.hh
@@ -1,339 +1,340 @@
/**
* @file aka_error.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Mon Jul 13 2015
*
* @brief error management and internal exceptions
*
* @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_ERROR_HH__
#define __AKANTU_ERROR_HH__
/* -------------------------------------------------------------------------- */
#include <cstdlib>
#include <ostream>
#include <sstream>
#include <typeinfo>
+#include <utility>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
enum DebugLevel {
dbl0 = 0,
dblError = 0,
dblAssert = 0,
dbl1 = 1,
dblException = 1,
dblCritical = 1,
dbl2 = 2,
dblMajor = 2,
dbl3 = 3,
dblCall = 3,
dblSecondary = 3,
dblHead = 3,
dbl4 = 4,
dblWarning = 4,
dbl5 = 5,
dblInfo = 5,
dbl6 = 6,
dblIn = 6,
dblOut = 6,
dbl7 = 7,
dbl8 = 8,
dblTrace = 8,
dbl9 = 9,
dblAccessory = 9,
dbl10 = 10,
dblDebug = 42,
dbl100 = 100,
dblDump = 100,
dblTest = 1337
};
/* -------------------------------------------------------------------------- */
#define AKANTU_LOCATION \
"(" << __func__ << "(): " << __FILE__ << ":" << __LINE__ << ")"
/* -------------------------------------------------------------------------- */
namespace debug {
void setDebugLevel(const DebugLevel & level);
const DebugLevel & getDebugLevel();
void initSignalHandler();
std::string demangle(const char * symbol);
std::string exec(std::string cmd);
void printBacktrace(int sig);
void exit(int status) __attribute__((noreturn));
/* ------------------------------------------------------------------------ */
/// exception class that can be thrown by akantu
class Exception : public std::exception {
/* ---------------------------------------------------------------------- */
/* Constructors/Destructors */
/* ---------------------------------------------------------------------- */
protected:
- Exception(std::string info = "") : _info(info), _file(""), _line(0) {}
+ Exception(std::string info = "") : _info(std::move(info)), _file("") {}
public:
//! full constructor
Exception(std::string info, std::string file, unsigned int line)
- : _info(info), _file(file), _line(line) {}
+ : _info(std::move(info)), _file(std::move(file)), _line(line) {}
//! destructor
- virtual ~Exception() throw(){};
+ ~Exception() throw() override = default;
/* ---------------------------------------------------------------------- */
/* Methods */
/* ---------------------------------------------------------------------- */
public:
- virtual const char * what() const throw() { return _info.c_str(); }
+ const char * what() const throw() override { return _info.c_str(); }
virtual const char * info() const throw() {
std::stringstream stream;
stream << debug::demangle(typeid(*this).name()) << " : " << _info << " ["
<< _file << ":" << _line << "]";
return stream.str().c_str();
}
public:
void setInfo(std::string info) { _info = info; }
void setFile(std::string file) { _file = file; }
void setLine(unsigned int line) { _line = line; }
/* ---------------------------------------------------------------------- */
/* Class Members */
/* ---------------------------------------------------------------------- */
private:
/// exception description and additionals
std::string _info;
/// file it is thrown from
std::string _file;
/// ligne it is thrown from
- unsigned int _line;
+ unsigned int _line{0};
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Exception & _this) {
stream << _this.what();
return stream;
}
/* --------------------------------------------------------------------------
*/
class Debugger {
public:
Debugger();
virtual ~Debugger();
void exit(int status) __attribute__((noreturn));
void throwException(const std::string & info, const std::string & file,
unsigned int line, __attribute__((unused)) bool silent,
__attribute__((unused))
const std::string & location) const
throw(akantu::debug::Exception) __attribute__((noreturn));
/*----------------------------------------------------------------------- */
template <class Except>
void throwCustomException(const Except & ex, const std::string & info,
const std::string & file, unsigned int line) const
throw(Except) __attribute__((noreturn));
/*----------------------------------------------------------------------- */
template <class Except>
void throwCustomException(const Except & ex, const std::string & file,
unsigned int line) const throw(Except)
__attribute__((noreturn));
void printMessage(const std::string & prefix, const DebugLevel & level,
const std::string & info) const;
void setOutStream(std::ostream & out) { cout = &out; }
std::ostream & getOutStream() { return *cout; }
public:
void setParallelContext(int rank, int size);
void setDebugLevel(const DebugLevel & level);
const DebugLevel & getDebugLevel() const;
void setLogFile(const std::string & filename);
std::ostream & getOutputStream();
inline bool testLevel(const DebugLevel & level) const {
return (this->level >= (level));
}
void printBacktrace(bool on_off) { this->print_backtrace = on_off; }
bool printBacktrace() { return this->print_backtrace; }
private:
std::string parallel_context;
std::ostream * cout;
bool file_open;
DebugLevel level;
bool print_backtrace;
};
extern Debugger debugger;
}
/* -------------------------------------------------------------------------- */
#define AKANTU_STRINGSTREAM_IN(_str, _sstr) \
; \
do { \
std::stringstream _dbg_s_info; \
_dbg_s_info << _sstr; \
_str = _dbg_s_info.str(); \
} while (false)
/* -------------------------------------------------------------------------- */
#define AKANTU_EXCEPTION(info) AKANTU_EXCEPTION_(info, false)
#define AKANTU_SILENT_EXCEPTION(info) AKANTU_EXCEPTION_(info, true)
#define AKANTU_EXCEPTION_(info, silent) \
do { \
std::stringstream _dbg_str; \
_dbg_str << info; \
std::stringstream _dbg_loc; \
_dbg_loc << AKANTU_LOCATION; \
::akantu::debug::debugger.throwException( \
_dbg_str.str(), __FILE__, __LINE__, silent, _dbg_loc.str()); \
} while (false)
#define AKANTU_CUSTOM_EXCEPTION_INFO(ex, info) \
do { \
std::stringstream _dbg_str; \
_dbg_str << info; \
::akantu::debug::debugger.throwCustomException(ex, _dbg_str.str(), \
__FILE__, __LINE__); \
} while (false)
#define AKANTU_CUSTOM_EXCEPTION(ex) \
do { \
::akantu::debug::debugger.throwCustomException(ex, __FILE__, __LINE__); \
} while (false)
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_NDEBUG
#define AKANTU_DEBUG_TEST(level) (false)
#define AKANTU_DEBUG_LEVEL_IS_TEST() \
(::akantu::debug::debugger.testLevel(dblTest))
#define AKANTU_DEBUG(level, info)
#define AKANTU_DEBUG_(pref, level, info)
#define AKANTU_DEBUG_IN()
#define AKANTU_DEBUG_OUT()
#define AKANTU_DEBUG_INFO(info)
#define AKANTU_DEBUG_WARNING(info)
#define AKANTU_DEBUG_TRACE(info)
#define AKANTU_DEBUG_ASSERT(test, info)
#define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info)
#define AKANTU_DEBUG_TO_IMPLEMENT() ::akantu::debug::debugger.exit(EXIT_FAILURE)
/* -------------------------------------------------------------------------- */
#else
#define AKANTU_DEBUG(level, info) AKANTU_DEBUG_(" ", level, info)
#define AKANTU_DEBUG_(pref, level, info) \
do { \
std::string _dbg_str; \
AKANTU_STRINGSTREAM_IN(_dbg_str, info << " " << AKANTU_LOCATION); \
::akantu::debug::debugger.printMessage(pref, level, _dbg_str); \
} while (false)
#define AKANTU_DEBUG_TEST(level) (::akantu::debug::debugger.testLevel(level))
#define AKANTU_DEBUG_LEVEL_IS_TEST() \
(::akantu::debug::debugger.testLevel(dblTest))
#define AKANTU_DEBUG_IN() \
AKANTU_DEBUG_("==>", ::akantu::dblIn, __func__ << "()")
#define AKANTU_DEBUG_OUT() \
AKANTU_DEBUG_("<==", ::akantu::dblOut, __func__ << "()")
#define AKANTU_DEBUG_INFO(info) AKANTU_DEBUG_("---", ::akantu::dblInfo, info)
#define AKANTU_DEBUG_WARNING(info) \
AKANTU_DEBUG_("/!\\", ::akantu::dblWarning, info)
#define AKANTU_DEBUG_TRACE(info) AKANTU_DEBUG_(">>>", ::akantu::dblTrace, info)
#define AKANTU_DEBUG_ASSERT(test, info) \
do { \
if (!(test)) { \
AKANTU_DEBUG_("!!! ", ::akantu::dblAssert, "assert [" << #test << "] " \
<< info); \
::akantu::debug::debugger.exit(EXIT_FAILURE); \
} \
} while (false)
#define AKANTU_DEBUG_ERROR(info) \
do { \
AKANTU_DEBUG_("!!! ", ::akantu::dblError, info); \
::akantu::debug::debugger.exit(EXIT_FAILURE); \
} while (false)
#define AKANTU_DEBUG_TO_IMPLEMENT() \
AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !")
#endif // AKANTU_NDEBUG
/* -------------------------------------------------------------------------- */
namespace debug {
/* ---------------------------------------------------------------------- */
template <class Except>
void Debugger::throwCustomException(const Except & ex,
const std::string & info,
const std::string & file,
unsigned int line) const throw(Except) {
- Except & nc_ex = const_cast<Except &>(ex);
+ auto & nc_ex = const_cast<Except &>(ex);
nc_ex.setInfo(info);
nc_ex.setFile(file);
nc_ex.setLine(line);
throw ex;
}
/* ----------------------------------------------------------------------- */
template <class Except>
void Debugger::throwCustomException(const Except & ex,
const std::string & file,
unsigned int line) const throw(Except) {
- Except & nc_ex = const_cast<Except &>(ex);
+ auto & nc_ex = const_cast<Except &>(ex);
nc_ex.setFile(file);
nc_ex.setLine(line);
throw ex;
}
}
}
#endif /* __AKANTU_ERROR_HH__ */
diff --git a/src/common/aka_math.hh b/src/common/aka_math.hh
index c1b5840c2..44389f462 100644
--- a/src/common/aka_math.hh
+++ b/src/common/aka_math.hh
@@ -1,289 +1,291 @@
/**
* @file aka_math.hh
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Marion Estelle Chambart <marion.chambart@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 Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Aug 04 2010
* @date last modification: Fri May 15 2015
*
* @brief mathematical operations
*
* @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_MATH_H__
#define __AKANTU_AKA_MATH_H__
/* -------------------------------------------------------------------------- */
+#include <utility>
+
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template <typename T, bool is_scal> class Array;
class Math {
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Matrix algebra */
/* ------------------------------------------------------------------------ */
/// @f$ y = A*x @f$
static void matrix_vector(UInt m, UInt n, const Array<Real, true> & A,
const Array<Real, true> & x, Array<Real, true> & y,
Real alpha = 1.);
/// @f$ y = A*x @f$
static inline void matrix_vector(UInt m, UInt n, Real * A, Real * x, Real * y,
Real alpha = 1.);
/// @f$ y = A^t*x @f$
static inline void matrixt_vector(UInt m, UInt n, Real * A, Real * x,
Real * y, Real alpha = 1.);
/// @f$ C = A*B @f$
static void matrix_matrix(UInt m, UInt n, UInt k, const Array<Real, true> & A,
const Array<Real, true> & B, Array<Real, true> & C,
Real alpha = 1.);
/// @f$ C = A*B^t @f$
static void matrix_matrixt(UInt m, UInt n, UInt k,
const Array<Real, true> & A,
const Array<Real, true> & B, Array<Real, true> & C,
Real alpha = 1.);
/// @f$ C = A*B @f$
static inline void matrix_matrix(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
/// @f$ C = A^t*B @f$
static inline void matrixt_matrix(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
/// @f$ C = A*B^t @f$
static inline void matrix_matrixt(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
/// @f$ C = A^t*B^t @f$
static inline void matrixt_matrixt(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
template <bool tr_A, bool tr_B>
static inline void matMul(UInt m, UInt n, UInt k, Real alpha, Real * A,
Real * B, Real beta, Real * C);
template <bool tr_A>
static inline void matVectMul(UInt m, UInt n, Real alpha, Real * A, Real * x,
Real beta, Real * y);
static inline void aXplusY(UInt n, Real alpha, Real * x, Real * y);
static inline void matrix33_eigenvalues(Real * A, Real * Adiag);
static inline void matrix22_eigenvalues(Real * A, Real * Adiag);
template <UInt dim> static inline void eigenvalues(Real * A, Real * d);
/// solve @f$ A x = \Lambda x @f$ and return d and V such as @f$ A V[i:] =
/// d[i] V[i:]@f$
template <typename T>
- static void matrixEig(UInt n, T * A, T * d, T * V = NULL);
+ static void matrixEig(UInt n, T * A, T * d, T * V = nullptr);
/// determinent of a 2x2 matrix
static inline Real det2(const Real * mat);
/// determinent of a 3x3 matrix
static inline Real det3(const Real * mat);
/// determinent of a nxn matrix
template <UInt n> static inline Real det(const Real * mat);
/// determinent of a nxn matrix
template <typename T> static inline T det(UInt n, const T * mat);
/// inverse a nxn matrix
template <UInt n> static inline void inv(const Real * mat, Real * inv);
/// inverse a nxn matrix
template <typename T> static inline void inv(UInt n, const T * mat, T * inv);
/// inverse a 3x3 matrix
static inline void inv3(const Real * mat, Real * inv);
/// inverse a 2x2 matrix
static inline void inv2(const Real * mat, Real * inv);
/// solve A x = b using a LU factorization
template <typename T>
static inline void solve(UInt n, const T * A, T * x, const T * b);
/// return the double dot product between 2 tensors in 2d
static inline Real matrixDoubleDot22(Real * A, Real * B);
/// return the double dot product between 2 tensors in 3d
static inline Real matrixDoubleDot33(Real * A, Real * B);
/// extension of the double dot product to two 2nd order tensor in dimension n
static inline Real matrixDoubleDot(UInt n, Real * A, Real * B);
/* ------------------------------------------------------------------------ */
/* Array algebra */
/* ------------------------------------------------------------------------ */
/// vector cross product
static inline void vectorProduct3(const Real * v1, const Real * v2,
Real * res);
/// normalize a vector
static inline void normalize2(Real * v);
/// normalize a vector
static inline void normalize3(Real * v);
/// return norm of a 2-vector
static inline Real norm2(const Real * v);
/// return norm of a 3-vector
static inline Real norm3(const Real * v);
/// return norm of a vector
static inline Real norm(UInt n, const Real * v);
/// return the dot product between 2 vectors in 2d
static inline Real vectorDot2(const Real * v1, const Real * v2);
/// return the dot product between 2 vectors in 3d
static inline Real vectorDot3(const Real * v1, const Real * v2);
/// return the dot product between 2 vectors
static inline Real vectorDot(Real * v1, Real * v2, UInt n);
/* ------------------------------------------------------------------------ */
/* Geometry */
/* ------------------------------------------------------------------------ */
/// compute normal a normal to a vector
static inline void normal2(const Real * v1, Real * res);
/// compute normal a normal to a vector
static inline void normal3(const Real * v1, const Real * v2, Real * res);
/// compute the tangents to an array of normal vectors
static void compute_tangents(const Array<Real> & normals,
Array<Real> & tangents);
/// distance in 2D between x and y
static inline Real distance_2d(const Real * x, const Real * y);
/// distance in 3D between x and y
static inline Real distance_3d(const Real * x, const Real * y);
/// radius of the in-circle of a triangle
static inline Real triangle_inradius(const Real * coord1, const Real * coord2,
const Real * coord3);
/// radius of the in-circle of a tetrahedron
static inline Real tetrahedron_inradius(const Real * coord1,
const Real * coord2,
const Real * coord3,
const Real * coord4);
/// volume of a tetrahedron
static inline Real tetrahedron_volume(const Real * coord1,
const Real * coord2,
const Real * coord3,
const Real * coord4);
/// compute the barycenter of n points
static inline void barycenter(const Real * coord, UInt nb_points,
UInt spatial_dimension, Real * barycenter);
/// vector between x and y
static inline void vector_2d(const Real * x, const Real * y, Real * vec);
/// vector pointing from x to y in 3 spatial dimension
static inline void vector_3d(const Real * x, const Real * y, Real * vec);
/// test if two scalar are equal within a given tolerance
static inline bool are_float_equal(Real x, Real y);
/// test if two vectors are equal within a given tolerance
static inline bool are_vector_equal(UInt n, Real * x, Real * y);
#ifdef isnan
#error \
"You probably included <math.h> which is incompatible with aka_math please use\
<cmath> or add a \"#undef isnan\" before akantu includes"
#endif
/// test if a real is a NaN
static inline bool isnan(Real x);
/// test if the line x and y intersects each other
static inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max);
/// test if a is in the range [x_min, x_max]
static inline bool is_in_range(Real a, Real x_min, Real x_max);
static inline Real getTolerance() { return tolerance; };
static inline void setTolerance(Real tol) { tolerance = tol; };
template <UInt p, typename T> static inline T pow(T x);
/// reduce all the values of an array, the summation is done in place and the
/// array is modified
static Real reduce(Array<Real> & array);
class NewtonRaphson {
public:
NewtonRaphson(Real tolerance, Real max_iteration)
: tolerance(tolerance), max_iteration(max_iteration) {}
template <class Functor> Real solve(const Functor & funct, Real x_0);
private:
Real tolerance;
Real max_iteration;
};
struct NewtonRaphsonFunctor {
- NewtonRaphsonFunctor(const std::string & name) : name(name) {}
+ NewtonRaphsonFunctor(std::string name) : name(std::move(name)) {}
virtual Real f(Real x) const = 0;
virtual Real f_prime(Real x) const = 0;
std::string name;
};
private:
/// tolerance for functions that need one
static Real tolerance;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "aka_math_tmpl.hh"
__END_AKANTU__
#endif /* __AKANTU_AKA_MATH_H__ */
diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh
index 5813430e3..18f64791e 100644
--- a/src/common/aka_math_tmpl.hh
+++ b/src/common/aka_math_tmpl.hh
@@ -1,786 +1,786 @@
/**
* @file aka_math_tmpl.hh
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Mathilde Radiguet <mathilde.radiguet@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Aug 04 2010
* @date last modification: Wed Oct 21 2015
*
* @brief Implementation of the inline functions of the math toolkit
*
* @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/>.
*
*/
__END_AKANTU__
#include <cmath>
#include <cstring>
#include <typeinfo>
#include "aka_blas_lapack.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
inline void Math::matrix_vector(UInt im, UInt in, Real * A, Real * x, Real * y,
Real alpha) {
#ifdef AKANTU_USE_BLAS
/// y = alpha*op(A)*x + beta*y
char tran_A = 'N';
int incx = 1;
int incy = 1;
double beta = 0.;
int m = im;
int n = in;
aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy);
#else
memset(y, 0, im * sizeof(Real));
for (UInt i = 0; i < im; ++i) {
for (UInt j = 0; j < in; ++j) {
y[i] += A[i + j * im] * x[j];
}
y[i] *= alpha;
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrixt_vector(UInt im, UInt in, Real * A, Real * x, Real * y,
Real alpha) {
#ifdef AKANTU_USE_BLAS
/// y = alpha*op(A)*x + beta*y
char tran_A = 'T';
int incx = 1;
int incy = 1;
double beta = 0.;
int m = im;
int n = in;
aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy);
#else
memset(y, 0, in * sizeof(Real));
for (UInt i = 0; i < im; ++i) {
for (UInt j = 0; j < in; ++j) {
y[j] += A[j * im + i] * x[i];
}
y[i] *= alpha;
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix_matrix(UInt im, UInt in, UInt ik, Real * A, Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'N';
char trans_b = 'N';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &k, &beta, C, &m);
#else
memset(C, 0, im * in * sizeof(Real));
for (UInt j = 0; j < in; ++j) {
UInt _jb = j * ik;
UInt _jc = j * im;
for (UInt i = 0; i < im; ++i) {
for (UInt l = 0; l < ik; ++l) {
UInt _la = l * im;
C[i + _jc] += A[i + _la] * B[l + _jb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrixt_matrix(UInt im, UInt in, UInt ik, Real * A, Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'T';
char trans_b = 'N';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &k, &beta, C, &m);
#else
memset(C, 0, im * in * sizeof(Real));
for (UInt j = 0; j < in; ++j) {
UInt _jc = j * im;
UInt _jb = j * ik;
for (UInt i = 0; i < im; ++i) {
UInt _ia = i * ik;
for (UInt l = 0; l < ik; ++l) {
C[i + _jc] += A[l + _ia] * B[l + _jb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix_matrixt(UInt im, UInt in, UInt ik, Real * A, Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'N';
char trans_b = 'T';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &n, &beta, C, &m);
#else
memset(C, 0, im * in * sizeof(Real));
for (UInt j = 0; j < in; ++j) {
UInt _jc = j * im;
for (UInt i = 0; i < im; ++i) {
for (UInt l = 0; l < ik; ++l) {
UInt _la = l * im;
UInt _lb = l * in;
C[i + _jc] += A[i + _la] * B[j + _lb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrixt_matrixt(UInt im, UInt in, UInt ik, Real * A, Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'T';
char trans_b = 'T';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &n, &beta, C, &m);
#else
memset(C, 0, im * in * sizeof(Real));
for (UInt j = 0; j < in; ++j) {
UInt _jc = j * im;
for (UInt i = 0; i < im; ++i) {
UInt _ia = i * ik;
for (UInt l = 0; l < ik; ++l) {
UInt _lb = l * in;
C[i + _jc] += A[l + _ia] * B[j + _lb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::aXplusY(UInt n, Real alpha, Real * x, Real * y) {
#ifdef AKANTU_USE_BLAS
/// y := alpha x + y
int incx = 1, incy = 1;
aka_axpy(&n, &alpha, x, &incx, y, &incy);
#else
for (UInt i = 0; i < n; ++i)
*(y++) += alpha * *(x++);
#endif
}
/* -------------------------------------------------------------------------- */
inline Real Math::vectorDot(Real * v1, Real * v2, UInt in) {
#ifdef AKANTU_USE_BLAS
/// d := v1 . v2
int incx = 1, incy = 1, n = in;
Real d = aka_dot(&n, v1, &incx, v2, &incy);
#else
Real d = 0;
for (UInt i = 0; i < in; ++i) {
d += v1[i] * v2[i];
}
#endif
return d;
}
/* -------------------------------------------------------------------------- */
template <bool tr_A, bool tr_B>
inline void Math::matMul(UInt m, UInt n, UInt k, Real alpha, Real * A, Real * B,
__attribute__((unused)) Real beta, Real * C) {
if (tr_A) {
if (tr_B)
matrixt_matrixt(m, n, k, A, B, C, alpha);
else
matrixt_matrix(m, n, k, A, B, C, alpha);
} else {
if (tr_B)
matrix_matrixt(m, n, k, A, B, C, alpha);
else
matrix_matrix(m, n, k, A, B, C, alpha);
}
}
/* -------------------------------------------------------------------------- */
template <bool tr_A>
inline void Math::matVectMul(UInt m, UInt n, Real alpha, Real * A, Real * x,
__attribute__((unused)) Real beta, Real * y) {
if (tr_A) {
matrixt_vector(m, n, A, x, y, alpha);
} else {
matrix_vector(m, n, A, x, y, alpha);
}
}
/* -------------------------------------------------------------------------- */
template <typename T> inline void Math::matrixEig(UInt n, T * A, T * d, T * V) {
// Matrix A is row major, so the lapack function in fortran will process
// A^t. Asking for the left eigenvectors of A^t will give the transposed right
// eigenvectors of A so in the C++ code the right eigenvectors.
char jobvl;
- if (V != NULL)
+ if (V != nullptr)
jobvl = 'V'; // compute left eigenvectors
else
jobvl = 'N'; // compute left eigenvectors
char jobvr('N'); // compute right eigenvectors
- T * di = new T[n]; // imaginary part of the eigenvalues
+ auto * di = new T[n]; // imaginary part of the eigenvalues
int info;
int N = n;
T wkopt;
int lwork = -1;
// query and allocate the optimal workspace
- aka_geev<T>(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, &wkopt, &lwork,
- &info);
+ aka_geev<T>(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, nullptr, &N, &wkopt,
+ &lwork, &info);
lwork = int(wkopt);
- T * work = new T[lwork];
+ auto * work = new T[lwork];
// solve the eigenproblem
- aka_geev<T>(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, work, &lwork,
- &info);
+ aka_geev<T>(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, nullptr, &N, work,
+ &lwork, &info);
AKANTU_DEBUG_ASSERT(
info == 0,
"Problem computing eigenvalues/vectors. DGEEV exited with the value "
<< info);
delete[] work;
delete[] di; // I hope for you that there was no complex eigenvalues !!!
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix22_eigenvalues(Real * A, Real * Adiag) {
/// d = determinant of Matrix A
Real d = det2(A);
/// b = trace of Matrix A
Real b = A[0] + A[3];
Real c = sqrt(b * b - 4 * d);
Adiag[0] = .5 * (b + c);
Adiag[1] = .5 * (b - c);
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix33_eigenvalues(Real * A, Real * Adiag) {
matrixEig(3, A, Adiag);
}
/* -------------------------------------------------------------------------- */
template <UInt dim> inline void Math::eigenvalues(Real * A, Real * d) {
if (dim == 1) {
d[0] = A[0];
} else if (dim == 2) {
matrix22_eigenvalues(A, d);
}
// else if(dim == 3) { matrix33_eigenvalues(A, d); }
else
matrixEig(dim, A, d);
}
/* -------------------------------------------------------------------------- */
inline Real Math::det2(const Real * mat) {
return mat[0] * mat[3] - mat[1] * mat[2];
}
/* -------------------------------------------------------------------------- */
inline Real Math::det3(const Real * mat) {
return mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) -
mat[3] * (mat[1] * mat[8] - mat[7] * mat[2]) +
mat[6] * (mat[1] * mat[5] - mat[4] * mat[2]);
}
/* -------------------------------------------------------------------------- */
template <UInt n> inline Real Math::det(const Real * mat) {
if (n == 1)
return *mat;
else if (n == 2)
return det2(mat);
else if (n == 3)
return det3(mat);
else
return det(n, mat);
}
/* -------------------------------------------------------------------------- */
template <typename T> inline T Math::det(UInt n, const T * A) {
int N = n;
int info;
- int * ipiv = new int[N + 1];
+ auto * ipiv = new int[N + 1];
- T * LU = new T[N * N];
+ auto * LU = new T[N * N];
std::copy(A, A + N * N, LU);
// LU factorization of A
aka_getrf(&N, &N, LU, &N, ipiv, &info);
if (info > 0) {
AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info
<< " )");
}
// det(A) = det(L) * det(U) = 1 * det(U) = product_i U_{ii}
T det = 1.;
for (int i = 0; i < N; ++i)
det *= (2 * (ipiv[i] == i) - 1) * LU[i * n + i];
delete[] ipiv;
delete[] LU;
return det;
}
/* -------------------------------------------------------------------------- */
inline void Math::normal2(const Real * vec, Real * normal) {
normal[0] = vec[1];
normal[1] = -vec[0];
Math::normalize2(normal);
}
/* -------------------------------------------------------------------------- */
inline void Math::normal3(const Real * vec1, const Real * vec2, Real * normal) {
Math::vectorProduct3(vec1, vec2, normal);
Math::normalize3(normal);
}
/* -------------------------------------------------------------------------- */
inline void Math::normalize2(Real * vec) {
Real norm = Math::norm2(vec);
vec[0] /= norm;
vec[1] /= norm;
}
/* -------------------------------------------------------------------------- */
inline void Math::normalize3(Real * vec) {
Real norm = Math::norm3(vec);
vec[0] /= norm;
vec[1] /= norm;
vec[2] /= norm;
}
/* -------------------------------------------------------------------------- */
inline Real Math::norm2(const Real * vec) {
return sqrt(vec[0] * vec[0] + vec[1] * vec[1]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::norm3(const Real * vec) {
return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::norm(UInt n, const Real * vec) {
Real norm = 0.;
for (UInt i = 0; i < n; ++i) {
norm += vec[i] * vec[i];
}
return sqrt(norm);
}
/* -------------------------------------------------------------------------- */
inline void Math::inv2(const Real * mat, Real * inv) {
Real det_mat = det2(mat);
inv[0] = mat[3] / det_mat;
inv[1] = -mat[1] / det_mat;
inv[2] = -mat[2] / det_mat;
inv[3] = mat[0] / det_mat;
}
/* -------------------------------------------------------------------------- */
inline void Math::inv3(const Real * mat, Real * inv) {
Real det_mat = det3(mat);
inv[0] = (mat[4] * mat[8] - mat[7] * mat[5]) / det_mat;
inv[1] = (mat[2] * mat[7] - mat[8] * mat[1]) / det_mat;
inv[2] = (mat[1] * mat[5] - mat[4] * mat[2]) / det_mat;
inv[3] = (mat[5] * mat[6] - mat[8] * mat[3]) / det_mat;
inv[4] = (mat[0] * mat[8] - mat[6] * mat[2]) / det_mat;
inv[5] = (mat[2] * mat[3] - mat[5] * mat[0]) / det_mat;
inv[6] = (mat[3] * mat[7] - mat[6] * mat[4]) / det_mat;
inv[7] = (mat[1] * mat[6] - mat[7] * mat[0]) / det_mat;
inv[8] = (mat[0] * mat[4] - mat[3] * mat[1]) / det_mat;
}
/* -------------------------------------------------------------------------- */
template <UInt n> inline void Math::inv(const Real * A, Real * Ainv) {
if (n == 1)
*Ainv = 1. / *A;
else if (n == 2)
inv2(A, Ainv);
else if (n == 3)
inv3(A, Ainv);
else
inv(n, A, Ainv);
}
/* -------------------------------------------------------------------------- */
template <typename T> inline void Math::inv(UInt n, const T * A, T * invA) {
int N = n;
int info;
- int * ipiv = new int[N + 1];
+ auto * ipiv = new int[N + 1];
int lwork = N * N;
- T * work = new T[lwork];
+ auto * work = new T[lwork];
std::copy(A, A + n * n, invA);
aka_getrf(&N, &N, invA, &N, ipiv, &info);
if (info > 0) {
AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info
<< " )");
}
aka_getri(&N, invA, &N, ipiv, work, &lwork, &info);
if (info != 0) {
AKANTU_DEBUG_ERROR("Cannot invert the matrix (info: " << info << " )");
}
delete[] ipiv;
delete[] work;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Math::solve(UInt n, const T * A, T * x, const T * b) {
int N = n;
int info;
- int * ipiv = new int[N];
- T * lu_A = new T[N * N];
+ auto * ipiv = new int[N];
+ auto * lu_A = new T[N * N];
std::copy(A, A + N * N, lu_A);
aka_getrf(&N, &N, lu_A, &N, ipiv, &info);
if (info > 0) {
AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info
<< " )");
exit(EXIT_FAILURE);
}
char trans = 'N';
int nrhs = 1;
std::copy(b, b + N, x);
aka_getrs(&trans, &N, &nrhs, lu_A, &N, ipiv, x, &N, &info);
if (info != 0) {
AKANTU_DEBUG_ERROR("Cannot solve the system (info: " << info << " )");
exit(EXIT_FAILURE);
}
delete[] ipiv;
delete[] lu_A;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
inline Real Math::matrixDoubleDot22(Real * A, Real * B) {
Real d;
d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3];
return d;
}
/* -------------------------------------------------------------------------- */
inline Real Math::matrixDoubleDot33(Real * A, Real * B) {
Real d;
d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3] + A[4] * B[4] +
A[5] * B[5] + A[6] * B[6] + A[7] * B[7] + A[8] * B[8];
return d;
}
/* -------------------------------------------------------------------------- */
inline Real Math::matrixDoubleDot(UInt n, Real * A, Real * B) {
Real d = 0.;
for (UInt i = 0; i < n; ++i) {
for (UInt j = 0; j < n; ++j) {
d += A[i * n + j] * B[i * n + j];
}
}
return d;
}
/* -------------------------------------------------------------------------- */
inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) {
res[0] = v1[1] * v2[2] - v1[2] * v2[1];
res[1] = v1[2] * v2[0] - v1[0] * v2[2];
res[2] = v1[0] * v2[1] - v1[1] * v2[0];
}
/* -------------------------------------------------------------------------- */
inline Real Math::vectorDot2(const Real * v1, const Real * v2) {
return (v1[0] * v2[0] + v1[1] * v2[1]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::vectorDot3(const Real * v1, const Real * v2) {
return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::distance_2d(const Real * x, const Real * y) {
return sqrt((y[0] - x[0]) * (y[0] - x[0]) + (y[1] - x[1]) * (y[1] - x[1]));
}
/* -------------------------------------------------------------------------- */
inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2,
const Real * coord3) {
/**
* @f{eqnarray*}{
* r &=& A / s \\
* A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\
* s &=& \frac{a + b + c}{2}
* @f}
*/
Real a, b, c;
a = distance_2d(coord1, coord2);
b = distance_2d(coord2, coord3);
c = distance_2d(coord1, coord3);
Real s;
s = (a + b + c) * 0.5;
return sqrt((s - a) * (s - b) * (s - c) / s);
}
/* -------------------------------------------------------------------------- */
inline Real Math::distance_3d(const Real * x, const Real * y) {
return sqrt((y[0] - x[0]) * (y[0] - x[0]) + (y[1] - x[1]) * (y[1] - x[1]) +
(y[2] - x[2]) * (y[2] - x[2]));
}
/* -------------------------------------------------------------------------- */
inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2,
const Real * coord3, const Real * coord4) {
Real xx[9], vol;
xx[0] = coord2[0];
xx[1] = coord2[1];
xx[2] = coord2[2];
xx[3] = coord3[0];
xx[4] = coord3[1];
xx[5] = coord3[2];
xx[6] = coord4[0];
xx[7] = coord4[1];
xx[8] = coord4[2];
vol = det3(xx);
xx[0] = coord1[0];
xx[1] = coord1[1];
xx[2] = coord1[2];
xx[3] = coord3[0];
xx[4] = coord3[1];
xx[5] = coord3[2];
xx[6] = coord4[0];
xx[7] = coord4[1];
xx[8] = coord4[2];
vol -= det3(xx);
xx[0] = coord1[0];
xx[1] = coord1[1];
xx[2] = coord1[2];
xx[3] = coord2[0];
xx[4] = coord2[1];
xx[5] = coord2[2];
xx[6] = coord4[0];
xx[7] = coord4[1];
xx[8] = coord4[2];
vol += det3(xx);
xx[0] = coord1[0];
xx[1] = coord1[1];
xx[2] = coord1[2];
xx[3] = coord2[0];
xx[4] = coord2[1];
xx[5] = coord2[2];
xx[6] = coord3[0];
xx[7] = coord3[1];
xx[8] = coord3[2];
vol -= det3(xx);
vol /= 6;
return vol;
}
/* -------------------------------------------------------------------------- */
inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2,
const Real * coord3,
const Real * coord4) {
Real l12, l13, l14, l23, l24, l34;
l12 = distance_3d(coord1, coord2);
l13 = distance_3d(coord1, coord3);
l14 = distance_3d(coord1, coord4);
l23 = distance_3d(coord2, coord3);
l24 = distance_3d(coord2, coord4);
l34 = distance_3d(coord3, coord4);
Real s1, s2, s3, s4;
s1 = (l12 + l23 + l13) * 0.5;
s1 = sqrt(s1 * (s1 - l12) * (s1 - l23) * (s1 - l13));
s2 = (l12 + l24 + l14) * 0.5;
s2 = sqrt(s2 * (s2 - l12) * (s2 - l24) * (s2 - l14));
s3 = (l23 + l34 + l24) * 0.5;
s3 = sqrt(s3 * (s3 - l23) * (s3 - l34) * (s3 - l24));
s4 = (l13 + l34 + l14) * 0.5;
s4 = sqrt(s4 * (s4 - l13) * (s4 - l34) * (s4 - l14));
Real volume = Math::tetrahedron_volume(coord1, coord2, coord3, coord4);
return 3 * volume / (s1 + s2 + s3 + s4);
}
/* -------------------------------------------------------------------------- */
inline void Math::barycenter(const Real * coord, UInt nb_points,
UInt spatial_dimension, Real * barycenter) {
memset(barycenter, 0, spatial_dimension * sizeof(Real));
for (UInt n = 0; n < nb_points; ++n) {
UInt offset = n * spatial_dimension;
for (UInt i = 0; i < spatial_dimension; ++i) {
barycenter[i] += coord[offset + i] / (Real)nb_points;
}
}
}
/* -------------------------------------------------------------------------- */
inline void Math::vector_2d(const Real * x, const Real * y, Real * res) {
res[0] = y[0] - x[0];
res[1] = y[1] - x[1];
}
/* -------------------------------------------------------------------------- */
inline void Math::vector_3d(const Real * x, const Real * y, Real * res) {
res[0] = y[0] - x[0];
res[1] = y[1] - x[1];
res[2] = y[2] - x[2];
}
/* -------------------------------------------------------------------------- */
/// Combined absolute and relative tolerance test proposed in
/// Real-time collision detection by C. Ericson (2004)
inline bool Math::are_float_equal(const Real x, const Real y) {
Real abs_max = std::max(std::abs(x), std::abs(y));
abs_max = std::max(abs_max, Real(1.));
return std::abs(x - y) <= (tolerance * abs_max);
}
/* -------------------------------------------------------------------------- */
inline bool Math::isnan(Real x) {
#if defined(__INTEL_COMPILER)
#pragma warning(push)
#pragma warning(disable : 1572)
#endif // defined(__INTEL_COMPILER)
// x = x return false means x = quiet_NaN
return !(x == x);
#if defined(__INTEL_COMPILER)
#pragma warning(pop)
#endif // defined(__INTEL_COMPILER)
}
/* -------------------------------------------------------------------------- */
inline bool Math::are_vector_equal(UInt n, Real * x, Real * y) {
bool test = true;
for (UInt i = 0; i < n; ++i) {
test &= are_float_equal(x[i], y[i]);
}
return test;
}
/* -------------------------------------------------------------------------- */
inline bool Math::intersects(Real x_min, Real x_max, Real y_min, Real y_max) {
return !((x_max <= y_min) || (x_min >= y_max));
}
/* -------------------------------------------------------------------------- */
inline bool Math::is_in_range(Real a, Real x_min, Real x_max) {
return ((a >= x_min) && (a <= x_max));
}
/* -------------------------------------------------------------------------- */
template <UInt p, typename T> inline T Math::pow(T x) {
return (pow<p - 1, T>(x) * x);
}
template <> inline UInt Math::pow<0, UInt>(__attribute__((unused)) UInt x) {
return (1);
}
template <> inline Real Math::pow<0, Real>(__attribute__((unused)) Real x) {
return (1.);
}
/* -------------------------------------------------------------------------- */
template <class Functor>
Real Math::NewtonRaphson::solve(const Functor & funct, Real x_0) {
Real x = x_0;
Real f_x = funct.f(x);
UInt iter = 0;
while (std::abs(f_x) > this->tolerance && iter < this->max_iteration) {
x -= f_x / funct.f_prime(x);
f_x = funct.f(x);
iter++;
}
AKANTU_DEBUG_ASSERT(iter < this->max_iteration,
"Newton Raphson ("
<< funct.name << ") solve did not converge in "
<< this->max_iteration << " iterations (tolerance: "
<< this->tolerance << ")");
return x;
}
diff --git a/src/common/aka_safe_enum.hh b/src/common/aka_safe_enum.hh
index 9c7dd64b6..fcbf9dee5 100644
--- a/src/common/aka_safe_enum.hh
+++ b/src/common/aka_safe_enum.hh
@@ -1,82 +1,83 @@
/**
* @file aka_safe_enum.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 21 2013
* @date last modification: Sun Oct 19 2014
*
* @brief Safe enums type (see More C++ Idioms/Type Safe Enum on Wikibooks
* http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Type_Safe_Enum)
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_SAFE_ENUM_HH__
#define __AKANTU_AKA_SAFE_ENUM_HH__
namespace akantu {
/// Safe enumerated type
template<typename def, typename inner = typename def::type>
class safe_enum : public def {
- typedef typename def::type type;
+ using type = typename def::type;
+
public:
explicit safe_enum(type v) : val(v) {}
inner underlying() const { return val; }
bool operator == (const safe_enum & s) const { return this->val == s.val; }
bool operator != (const safe_enum & s) const { return this->val != s.val; }
bool operator < (const safe_enum & s) const { return this->val < s.val; }
bool operator <= (const safe_enum & s) const { return this->val <= s.val; }
bool operator > (const safe_enum & s) const { return this->val > s.val; }
bool operator >= (const safe_enum & s) const { return this->val >= s.val; }
operator inner() { return val; };
public:
// Works only if enumerations are contiguous.
class iterator {
public:
explicit iterator(type v) : it(v) { }
void operator++() { ++it; }
safe_enum operator*() { return safe_enum(static_cast<type>(it)); }
bool operator!=(iterator const & it) { return it.it != this->it; }
private:
int it;
};
static iterator begin() {
return iterator(def::_begin_);
}
static iterator end() {
return iterator(def::_end_);
}
protected:
inner val;
};
} // akantu
#endif /* __AKANTU_AKA_SAFE_ENUM_HH__ */
diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index 5d7990535..3d5f70746 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,1185 +1,1185 @@
/**
* @file aka_types.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 17 2011
* @date last modification: Fri Jan 22 2016
*
* @brief description of the "simple" types
*
* @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 "aka_error.hh"
#include "aka_fwd.hh"
#include "aka_math.hh"
/* -------------------------------------------------------------------------- */
#include <iomanip>
#include <initializer_list>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_TYPES_HH__
#define __AKANTU_AKA_TYPES_HH__
__BEGIN_AKANTU__
enum NormType { L_1 = 1, L_2 = 2, L_inf = UInt(-1) };
/**
* DimHelper is a class to generalize the setup of a dim array from 3
* values. This gives a common interface in the TensorStorage class
* independently of its derived inheritance (Vector, Matrix, Tensor3)
* @tparam dim
*/
template <UInt dim> struct DimHelper {
static inline void setDims(UInt m, UInt n, UInt p, UInt dims[dim]);
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<1> {
static inline void setDims(UInt m, __attribute__((unused)) UInt n,
__attribute__((unused)) UInt p, UInt dims[1]) {
dims[0] = m;
}
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<2> {
static inline void setDims(UInt m, UInt n, __attribute__((unused)) UInt p,
UInt dims[2]) {
dims[0] = m;
dims[1] = n;
}
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<3> {
static inline void setDims(UInt m, UInt n, UInt p, UInt dims[3]) {
dims[0] = m;
dims[1] = n;
dims[2] = p;
}
};
/* -------------------------------------------------------------------------- */
template <typename T, UInt ndim, class RetType> class TensorStorage;
/* -------------------------------------------------------------------------- */
/* Proxy classes */
/* -------------------------------------------------------------------------- */
/**
* @class TensorProxy aka_types.hh
* @desc The TensorProxy class is a proxy class to the TensorStorage it handles
* the
* wrapped case. That is to say if an accessor should give access to a Tensor
* wrapped on some data, like the Array<T>::iterator they can return a
* TensorProxy that will be automatically transformed as a TensorStorage wrapped
* on the same data
* @tparam T stored type
* @tparam ndim order of the tensor
* @tparam RetType real derived type
*/
template <typename T, UInt ndim, class RetType> class TensorProxy {
protected:
TensorProxy(T * data, UInt m, UInt n, UInt p) {
DimHelper<ndim>::setDims(m, n, p, this->n);
this->values = data;
}
TensorProxy(const TensorProxy & other) {
this->values = other.storage();
for (UInt i = 0; i < ndim; ++i)
this->n[i] = other.n[i];
}
inline TensorProxy(const TensorStorage<T, ndim, RetType> & other);
- typedef typename RetType::proxy RetTypeProxy;
+ using RetTypeProxy = typename RetType::proxy;
public:
operator RetType() { return RetType(static_cast<RetTypeProxy &>(*this)); }
UInt size(UInt i) const {
AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
<< " dimensions, not "
<< (i + 1));
return n[i];
}
inline UInt size() const {
UInt _size = 1;
for (UInt d = 0; d < ndim; ++d)
_size *= this->n[d];
return _size;
}
T * storage() const { return values; }
inline RetTypeProxy & operator=(const RetType & src) {
AKANTU_DEBUG_ASSERT(
src.size() == this->size(),
"You are trying to copy two tensors with different sizes");
memcpy(this->values, src.storage(), this->size() * sizeof(T));
return *this;
}
inline RetTypeProxy & operator=(const RetTypeProxy & src) {
AKANTU_DEBUG_ASSERT(
src.size() == this->size(),
"You are trying to copy two tensors with different sizes");
memcpy(this->values, src.storage(), this->size() * sizeof(T));
return *this;
}
template <typename O> inline RetTypeProxy & operator*=(const O & o) {
RetType(*this) *= o;
return static_cast<RetTypeProxy &>(*this);
}
template <typename O> inline RetTypeProxy & operator/=(const O & o) {
RetType(*this) /= o;
return static_cast<RetTypeProxy &>(*this);
}
protected:
T * values;
UInt n[ndim];
};
/* -------------------------------------------------------------------------- */
template <typename T> class VectorProxy : public TensorProxy<T, 1, Vector<T> > {
typedef TensorProxy<T, 1, Vector<T> > parent;
- typedef Vector<T> type;
+ using type = Vector<T>;
public:
VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {}
VectorProxy(const VectorProxy & src) : parent(src) {}
VectorProxy(const Vector<T> & src) : parent(src) {}
T & operator()(UInt index) { return this->values[index]; };
const T & operator()(UInt index) const { return this->values[index]; };
};
template <typename T> class MatrixProxy : public TensorProxy<T, 2, Matrix<T> > {
typedef TensorProxy<T, 2, Matrix<T> > parent;
- typedef Matrix<T> type;
+ using type = Matrix<T>;
public:
MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
MatrixProxy(const MatrixProxy & src) : parent(src) {}
MatrixProxy(const type & src) : parent(src) {}
};
template <typename T>
class Tensor3Proxy : public TensorProxy<T, 3, Tensor3<T> > {
typedef TensorProxy<T, 3, Tensor3<T> > parent;
- typedef Tensor3<T> type;
+ using type = Tensor3<T>;
public:
Tensor3Proxy(T * data, UInt m, UInt n, UInt k) : parent(data, m, n, k) {}
Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {}
Tensor3Proxy(const Tensor3<T> & src) : parent(src) {}
};
/* -------------------------------------------------------------------------- */
/* Tensor base class */
/* -------------------------------------------------------------------------- */
template <typename T, UInt ndim, class RetType> class TensorStorage {
public:
- typedef T value_type;
+ using value_type = T;
protected:
template <class TensorType> void copySize(const TensorType & src) {
for (UInt d = 0; d < ndim; ++d)
this->n[d] = src.size(d);
this->_size = src.size();
}
- TensorStorage() : values(NULL), wrapped(false) {
+ TensorStorage() : values(NULL) {
for (UInt d = 0; d < ndim; ++d)
this->n[d] = 0;
_size = 0;
}
TensorStorage(const TensorProxy<T, ndim, RetType> & proxy) {
this->copySize(proxy);
this->values = proxy.storage();
this->wrapped = true;
}
protected:
TensorStorage(const TensorStorage & src) = delete;
public:
TensorStorage(const TensorStorage & src, bool deep_copy)
: values(NULL), wrapped(false) {
if (deep_copy)
this->deepCopy(src);
else
this->shallowCopy(src);
}
protected:
TensorStorage(UInt m, UInt n, UInt p, const T & def) {
DimHelper<ndim>::setDims(m, n, p, this->n);
this->computeSize();
this->values = new T[this->_size];
this->set(def);
this->wrapped = false;
}
TensorStorage(T * data, UInt m, UInt n, UInt p) {
DimHelper<ndim>::setDims(m, n, p, this->n);
this->computeSize();
this->values = data;
this->wrapped = true;
}
public:
/* ------------------------------------------------------------------------ */
template <class TensorType> inline void shallowCopy(const TensorType & src) {
this->copySize(src);
if (!this->wrapped)
delete[] this->values;
this->values = src.storage();
this->wrapped = true;
}
/* ------------------------------------------------------------------------ */
template <class TensorType> inline void deepCopy(const TensorType & src) {
this->copySize(src);
if (!this->wrapped)
delete[] this->values;
this->values = new T[this->_size];
memcpy((void *)this->values, (void *)src.storage(),
this->_size * sizeof(T));
this->wrapped = false;
}
virtual ~TensorStorage() {
if (!this->wrapped)
delete[] this->values;
}
/* ------------------------------------------------------------------------ */
inline TensorStorage & operator=(const RetType & src) {
if (this != &src) {
if (this->wrapped) {
// this test is not sufficient for Tensor of order higher than 1
AKANTU_DEBUG_ASSERT(this->_size == src.size(),
"Tensors of different size");
memcpy((void *)this->values, (void *)src.storage(),
this->_size * sizeof(T));
} else {
deepCopy(src);
}
}
return *this;
}
/* ------------------------------------------------------------------------ */
template <class R>
inline RetType & operator+=(const TensorStorage<T, ndim, R> & other) {
T * a = this->storage();
T * b = other.storage();
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be subtracted");
for (UInt i = 0; i < _size; ++i)
*(a++) += *(b++);
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
template <class R>
inline RetType & operator-=(const TensorStorage<T, ndim, R> & other) {
T * a = this->storage();
T * b = other.storage();
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be subtracted");
for (UInt i = 0; i < _size; ++i)
*(a++) -= *(b++);
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator+=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i)
*(a++) += x;
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator-=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i)
*(a++) -= x;
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator*=(const T & x) {
T * a = this->storage();
for (UInt i = 0; i < _size; ++i)
*(a++) *= x;
return *(static_cast<RetType *>(this));
}
/* ---------------------------------------------------------------------- */
inline RetType & operator/=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i)
*(a++) /= x;
return *(static_cast<RetType *>(this));
}
/// Y = \alpha X + Y
inline RetType & aXplusY(const TensorStorage & other, const T & alpha = 1.) {
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be subtracted");
Math::aXplusY(this->_size, alpha, other.storage(), this->storage());
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
T * storage() const { return values; }
UInt size() const { return _size; }
UInt size(UInt i) const {
AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
<< " dimensions, not "
<< (i + 1));
return n[i];
};
/* ------------------------------------------------------------------------ */
inline void clear() { memset(values, 0, _size * sizeof(T)); };
inline void set(const T & t) { std::fill_n(values, _size, t); };
template <class TensorType> inline void copy(const TensorType & other) {
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be copied");
memcpy(values, other.storage(), _size * sizeof(T));
}
bool isWrapped() const { return this->wrapped; }
protected:
friend class Array<T>;
inline void computeSize() {
_size = 1;
for (UInt d = 0; d < ndim; ++d)
_size *= this->n[d];
}
protected:
template <typename R, NormType norm_type> struct NormHelper {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it)
_norm += std::pow(std::abs(*it), norm_type);
return std::pow(_norm, 1. / norm_type);
}
};
template <typename R> struct NormHelper<R, L_1> {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it)
_norm += std::abs(*it);
return _norm;
}
};
template <typename R> struct NormHelper<R, L_2> {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it)
_norm += *it * *it;
return sqrt(_norm);
}
};
template <typename R> struct NormHelper<R, L_inf> {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it)
_norm = std::max(std::abs(*it), _norm);
return _norm;
}
};
public:
/*----------------------------------------------------------------------- */
/// "Entrywise" norm norm<L_p> @f[ \|\boldsymbol{T}\|_p = \left(
/// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}}
/// @f]
template <NormType norm_type> inline T norm() const {
return NormHelper<T, norm_type>::norm(*this);
}
protected:
UInt n[ndim];
UInt _size;
T * values;
- bool wrapped;
+ bool wrapped{false};
};
template <typename T, UInt ndim, class RetType>
inline TensorProxy<T, ndim, RetType>::TensorProxy(
const TensorStorage<T, ndim, RetType> & other) {
this->values = other.storage();
for (UInt i = 0; i < ndim; ++i)
this->n[i] = other.size(i);
}
/* -------------------------------------------------------------------------- */
/* Vector */
/* -------------------------------------------------------------------------- */
template <typename T> class Vector : public TensorStorage<T, 1, Vector<T> > {
typedef TensorStorage<T, 1, Vector<T> > parent;
public:
- typedef typename parent::value_type value_type;
- typedef VectorProxy<T> proxy;
+ using value_type = typename parent::value_type;
+ using proxy = VectorProxy<T>;
public:
Vector() : parent() {}
explicit Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) {}
Vector(T * data, UInt n) : parent(data, n, 0, 0) {}
Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) {}
Vector(const VectorProxy<T> & src) : parent(src) {}
Vector(std::initializer_list<T> list) : parent(list.size(), 0, 0, T()) {
UInt i = 0;
for(auto val : list) {
operator()(i++) = val;
}
}
public:
- virtual ~Vector(){};
+ ~Vector() override = default;
/* ------------------------------------------------------------------------ */
inline Vector & operator=(const Vector & src) {
parent::operator=(src);
return *this;
}
/* ------------------------------------------------------------------------ */
inline T & operator()(UInt i) {
AKANTU_DEBUG_ASSERT((i < this->n[0]),
"Access out of the vector! "
<< "Index (" << i
<< ") is out of the vector of size (" << this->n[0]
<< ")");
return *(this->values + i);
}
inline const T & operator()(UInt i) const {
AKANTU_DEBUG_ASSERT((i < this->n[0]),
"Access out of the vector! "
<< "Index (" << i
<< ") is out of the vector of size (" << this->n[0]
<< ")");
return *(this->values + i);
}
inline T & operator[](UInt i) { return this->operator()(i); }
inline const T & operator[](UInt i) const { return this->operator()(i); }
/* ------------------------------------------------------------------------ */
inline Vector<T> & operator*=(Real x) { return parent::operator*=(x); }
inline Vector<T> & operator/=(Real x) { return parent::operator/=(x); }
/* ------------------------------------------------------------------------ */
inline Vector<T> & operator*=(const Vector<T> & vect) {
T * a = this->storage();
T * b = vect.storage();
for (UInt i = 0; i < this->_size; ++i)
*(a++) *= *(b++);
return *this;
}
/* ------------------------------------------------------------------------ */
inline Real dot(const Vector<T> & vect) const {
return Math::vectorDot(this->values, vect.storage(), this->_size);
}
/* ------------------------------------------------------------------------ */
inline Real mean() const {
Real mean = 0;
T * a = this->storage();
for (UInt i = 0; i < this->_size; ++i)
mean += *(a++);
return mean / this->_size;
}
/* ------------------------------------------------------------------------ */
inline Vector & crossProduct(const Vector<T> & v1, const Vector<T> & v2) {
AKANTU_DEBUG_ASSERT(this->size() == 3,
"crossProduct is only defined in 3D (n=" << this->size()
<< ")");
AKANTU_DEBUG_ASSERT(
this->size() == v1.size() && this->size() == v2.size(),
"crossProduct is not a valid operation non matching size vectors");
Math::vectorProduct3(v1.storage(), v2.storage(), this->values);
return *this;
}
/* ------------------------------------------------------------------------ */
inline void solve(Matrix<T> & A, const Vector<T> & b) {
AKANTU_DEBUG_ASSERT(
this->size() == A.rows() && this->_size = A.cols(),
"The size of the solution vector mismatches the size of the matrix");
AKANTU_DEBUG_ASSERT(
this->_size == b._size,
"The rhs vector has a mismatch in size with the matrix");
Math::solve(this->_size, A.storage(), this->values, b.storage());
}
/* ------------------------------------------------------------------------ */
template <bool tr_A>
inline void mul(const Matrix<T> & A, const Vector<T> & x, Real alpha = 1.0);
/* ------------------------------------------------------------------------ */
inline Real norm() const { return parent::template norm<L_2>(); }
template <NormType nt> inline Real norm() const {
return parent::template norm<nt>();
}
/* ------------------------------------------------------------------------ */
inline void normalize() {
Real n = norm();
operator/=(n);
}
/* ------------------------------------------------------------------------ */
/// norm of (*this - x)
inline Real distance(const Vector<T> & y) const {
Real * vx = this->values;
Real * vy = y.storage();
Real sum_2 = 0;
for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy)
sum_2 += (*vx - *vy) * (*vx - *vy);
return sqrt(sum_2);
}
/* ------------------------------------------------------------------------ */
inline bool equal(const Vector<T> & v,
Real tolerance = Math::getTolerance()) const {
T * a = this->storage();
T * b = v.storage();
UInt i = 0;
while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance))
++i;
return i == this->_size;
}
/* ------------------------------------------------------------------------ */
inline short compare(const Vector<T> & v,
Real tolerance = Math::getTolerance()) const {
T * a = this->storage();
T * b = v.storage();
for (UInt i(0); i < this->_size; ++i, ++a, ++b) {
if (std::abs(*a - *b) > tolerance)
return (((*a - *b) > tolerance) ? 1 : -1);
}
return 0;
}
/* ------------------------------------------------------------------------ */
inline bool operator==(const Vector<T> & v) const { return equal(v); }
inline bool operator!=(const Vector<T> & v) const { return ! operator==(v); }
inline bool operator<(const Vector<T> & v) const { return compare(v) == -1; }
inline bool operator>(const Vector<T> & v) const { return compare(v) == 1; }
/* ------------------------------------------------------------------------ */
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << "[";
for (UInt i = 0; i < this->_size; ++i) {
if (i != 0)
stream << ", ";
stream << this->values[i];
}
stream << "]";
}
// friend class ::akantu::Array<T>;
};
-typedef Vector<Real> RVector;
+using RVector = Vector<Real>;
/* ------------------------------------------------------------------------ */
template <>
inline bool Vector<UInt>::equal(const Vector<UInt> & v,
__attribute__((unused)) Real tolerance) const {
UInt * a = this->storage();
UInt * b = v.storage();
UInt i = 0;
while (i < this->_size && (*(a++) == *(b++)))
++i;
return i == this->_size;
}
/* ------------------------------------------------------------------------ */
/* Matrix */
/* ------------------------------------------------------------------------ */
template <typename T> class Matrix : public TensorStorage<T, 2, Matrix<T> > {
typedef TensorStorage<T, 2, Matrix<T> > parent;
public:
- typedef typename parent::value_type value_type;
- typedef MatrixProxy<T> proxy;
+ using value_type = typename parent::value_type;
+ using proxy = MatrixProxy<T>;
public:
Matrix() : parent() {}
Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) {}
Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) {}
Matrix(const MatrixProxy<T> & src) : parent(src) {}
Matrix(std::initializer_list<std::initializer_list<T>> list) {
std::size_t n = 0;
std::size_t m = list.size();
for(auto row : list) {
n = std::max(n, row.size());
}
DimHelper<2>::setDims(m, n, 0, this->n);
this->computeSize();
this->values = new T[this->_size];
this->set(0);
UInt i = 0, j = 0;
for(auto row : list) {
for(auto val : row) {
at(i, j++) = val;
}
++i;
j = 0;
}
}
- virtual ~Matrix() {}
+ virtual ~Matrix() = default;
/* ------------------------------------------------------------------------ */
inline Matrix & operator=(const Matrix & src) {
parent::operator=(src);
return *this;
}
public:
/* ---------------------------------------------------------------------- */
UInt rows() const { return this->n[0]; }
UInt cols() const { return this->n[1]; }
/* ---------------------------------------------------------------------- */
inline T & at(UInt i, UInt j) {
AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])),
"Access out of the matrix! "
<< "Index (" << i << ", " << j
<< ") is out of the matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return *(this->values + i + j * this->n[0]);
}
inline const T & at(UInt i, UInt j) const {
AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])),
"Access out of the matrix! "
<< "Index (" << i << ", " << j
<< ") is out of the matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return *(this->values + i + j * this->n[0]);
}
/* ---------------------------------------------------------------------- */
inline T & operator()(UInt i, UInt j) { return this->at(i, j); }
inline const T & operator()(UInt i, UInt j) const { return this->at(i, j); }
/// give a line vector wrapped on the column i
inline VectorProxy<T> operator()(UInt j) {
AKANTU_DEBUG_ASSERT(j < this->n[1],
"Access out of the matrix! "
<< "You are trying to access the column vector "
<< j << " in a matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return VectorProxy<T>(this->values + j * this->n[0], this->n[0]);
}
inline const VectorProxy<T> operator()(UInt j) const {
AKANTU_DEBUG_ASSERT(j < this->n[1],
"Access out of the matrix! "
<< "You are trying to access the column vector "
<< j << " in a matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return VectorProxy<T>(this->values + j * this->n[0], this->n[0]);
}
inline T & operator[](UInt idx) { return *(this->values + idx); };
inline const T & operator[](UInt idx) const { return *(this->values + idx); };
/* ---------------------------------------------------------------------- */
inline Matrix operator*(const Matrix & B) {
Matrix C(this->rows(), B.cols());
C.mul<false, false>(*this, B);
return C;
}
/* ----------------------------------------------------------------------- */
inline Matrix & operator*=(const T & x) { return parent::operator*=(x); }
inline Matrix & operator*=(const Matrix & B) {
Matrix C(*this);
this->mul<false, false>(C, B);
return *this;
}
/* ---------------------------------------------------------------------- */
template <bool tr_A, bool tr_B>
inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) {
UInt k = A.cols();
if (tr_A)
k = A.rows();
#ifndef AKANTU_NDEBUG
if (tr_B) {
AKANTU_DEBUG_ASSERT(k == B.cols(),
"matrices to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->cols() == B.rows(),
"matrices to multiply have no fit dimensions");
} else {
AKANTU_DEBUG_ASSERT(k == B.rows(),
"matrices to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->cols() == B.cols(),
"matrices to multiply have no fit dimensions");
}
if (tr_A) {
AKANTU_DEBUG_ASSERT(this->rows() == A.cols(),
"matrices to multiply have no fit dimensions");
} else {
AKANTU_DEBUG_ASSERT(this->rows() == A.rows(),
"matrices to multiply have no fit dimensions");
}
#endif // AKANTU_NDEBUG
Math::matMul<tr_A, tr_B>(this->rows(), this->cols(), k, alpha, A.storage(),
B.storage(), 0., this->storage());
}
/* ---------------------------------------------------------------------- */
inline void outerProduct(const Vector<T> & A, const Vector<T> & B) {
AKANTU_DEBUG_ASSERT(
A.size() == this->rows() && B.size() == this->cols(),
"A and B are not compatible with the size of the matrix");
for (UInt i = 0; i < this->rows(); ++i) {
for (UInt j = 0; j < this->cols(); ++j) {
this->values[i + j * this->rows()] += A[i] * B[j];
}
}
}
private:
class EigenSorter {
public:
EigenSorter(const Vector<T> & eigs) : eigs(eigs) {}
bool operator()(const UInt & a, const UInt & b) const {
return (eigs(a) > eigs(b));
}
private:
const Vector<T> & eigs;
};
public:
/* ---------------------------------------------------------------------- */
inline void eig(Vector<T> & eigenvalues, Matrix<T> & eigenvectors) const {
AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
"eig is not a valid operation on a rectangular matrix");
AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(),
"eigenvalues should be of size " << this->cols()
<< ".");
#ifndef AKANTU_NDEBUG
if (eigenvectors.storage() != NULL)
AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) &&
(eigenvectors.rows() == this->cols()),
"Eigenvectors needs to be a square matrix of size "
<< this->cols() << " x " << this->cols() << ".");
#endif
Matrix<T> tmp = *this;
Vector<T> tmp_eigs(eigenvalues.size());
Matrix<T> tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols());
if (tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0)
Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage());
else
Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(),
tmp_eig_vects.storage());
Vector<UInt> perm(eigenvalues.size());
for (UInt i = 0; i < perm.size(); ++i)
perm(i) = i;
std::sort(perm.storage(), perm.storage() + perm.size(),
EigenSorter(tmp_eigs));
for (UInt i = 0; i < perm.size(); ++i)
eigenvalues(i) = tmp_eigs(perm(i));
if (tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0)
for (UInt i = 0; i < perm.size(); ++i) {
for (UInt j = 0; j < eigenvectors.rows(); ++j) {
eigenvectors(j, i) = tmp_eig_vects(j, perm(i));
}
}
}
/* ---------------------------------------------------------------------- */
inline void eig(Vector<T> & eigenvalues) const {
Matrix<T> empty;
eig(eigenvalues, empty);
}
/* ---------------------------------------------------------------------- */
inline void eye(T alpha = 1.) {
AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
"eye is not a valid operation on a rectangular matrix");
this->clear();
for (UInt i = 0; i < this->cols(); ++i) {
this->values[i + i * this->rows()] = alpha;
}
}
/* ---------------------------------------------------------------------- */
static inline Matrix<T> eye(UInt m, T alpha = 1.) {
Matrix<T> tmp(m, m);
tmp.eye(alpha);
return tmp;
}
/* ---------------------------------------------------------------------- */
inline T trace() const {
AKANTU_DEBUG_ASSERT(
this->cols() == this->rows(),
"trace is not a valid operation on a rectangular matrix");
T trace = 0.;
for (UInt i = 0; i < this->rows(); ++i) {
trace += this->values[i + i * this->rows()];
}
return trace;
}
/* ---------------------------------------------------------------------- */
inline Matrix transpose() const {
Matrix tmp(this->cols(), this->rows());
for (UInt i = 0; i < this->rows(); ++i) {
for (UInt j = 0; j < this->cols(); ++j) {
tmp(j, i) = operator()(i, j);
}
}
return tmp;
}
/* ---------------------------------------------------------------------- */
inline void inverse(const Matrix & A) {
AKANTU_DEBUG_ASSERT(A.cols() == A.rows(),
"inv is not a valid operation on a rectangular matrix");
AKANTU_DEBUG_ASSERT(this->cols() == A.cols(),
"the matrix should have the same size as its inverse");
if (this->cols() == 1)
*this->values = 1. / *A.storage();
else if (this->cols() == 2)
Math::inv2(A.storage(), this->values);
else if (this->cols() == 3)
Math::inv3(A.storage(), this->values);
else
Math::inv(this->cols(), A.storage(), this->values);
}
/* --------------------------------------------------------------------- */
inline T det() const {
AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
"inv is not a valid operation on a rectangular matrix");
if (this->cols() == 1)
return *(this->values);
else if (this->cols() == 2)
return Math::det2(this->values);
else if (this->cols() == 3)
return Math::det3(this->values);
else
return Math::det(this->cols(), this->values);
}
/* --------------------------------------------------------------------- */
inline T doubleDot(const Matrix<T> & other) const {
AKANTU_DEBUG_ASSERT(
this->cols() == this->rows(),
"doubleDot is not a valid operation on a rectangular matrix");
if (this->cols() == 1)
return *(this->values) * *(other.storage());
else if (this->cols() == 2)
return Math::matrixDoubleDot22(this->values, other.storage());
else if (this->cols() == 3)
return Math::matrixDoubleDot33(this->values, other.storage());
else
AKANTU_DEBUG_ERROR("doubleDot is not defined for other spatial dimensions"
<< " than 1, 2 or 3.");
return T();
}
/* ---------------------------------------------------------------------- */
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << "[";
for (UInt i = 0; i < this->n[0]; ++i) {
if (i != 0)
stream << ", ";
stream << "[";
for (UInt j = 0; j < this->n[1]; ++j) {
if (j != 0)
stream << ", ";
stream << operator()(i, j);
}
stream << "]";
}
stream << "]";
};
};
/* ------------------------------------------------------------------------ */
template <typename T>
template <bool tr_A>
inline void Vector<T>::mul(const Matrix<T> & A, const Vector<T> & x,
Real alpha) {
#ifndef AKANTU_NDEBUG
UInt n = x.size();
if (tr_A) {
AKANTU_DEBUG_ASSERT(n == A.rows(),
"matrix and vector to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->size() == A.cols(),
"matrix and vector to multiply have no fit dimensions");
} else {
AKANTU_DEBUG_ASSERT(n == A.cols(),
"matrix and vector to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->size() == A.rows(),
"matrix and vector to multiply have no fit dimensions");
}
#endif
Math::matVectMul<tr_A>(A.rows(), A.cols(), alpha, A.storage(), x.storage(),
0., this->storage());
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::ostream & operator<<(std::ostream & stream,
const Matrix<T> & _this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::ostream & operator<<(std::ostream & stream,
const Vector<T> & _this) {
_this.printself(stream);
return stream;
}
/* ------------------------------------------------------------------------ */
/* Tensor3 */
/* ------------------------------------------------------------------------ */
template <typename T> class Tensor3 : public TensorStorage<T, 3, Tensor3<T> > {
typedef TensorStorage<T, 3, Tensor3<T> > parent;
public:
- typedef typename parent::value_type value_type;
- typedef Tensor3Proxy<T> proxy;
+ using value_type = typename parent::value_type;
+ using proxy = Tensor3Proxy<T>;
public:
Tensor3() : parent(){};
Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) {}
Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) {}
Tensor3(const Tensor3 & src, bool deep_copy = true)
: parent(src, deep_copy) {}
public:
/* ------------------------------------------------------------------------ */
inline Tensor3 & operator=(const Tensor3 & src) {
parent::operator=(src);
return *this;
}
/* ---------------------------------------------------------------------- */
inline T & operator()(UInt i, UInt j, UInt k) {
AKANTU_DEBUG_ASSERT(
(i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the element "
<< "(" << i << ", " << j << ", " << k << ") in a tensor of size ("
<< this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")");
return *(this->values + (k * this->n[0] + i) * this->n[1] + j);
}
inline const T & operator()(UInt i, UInt j, UInt k) const {
AKANTU_DEBUG_ASSERT(
(i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the element "
<< "(" << i << ", " << j << ", " << k << ") in a tensor of size ("
<< this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")");
return *(this->values + (k * this->n[0] + i) * this->n[1] + j);
}
inline MatrixProxy<T> operator()(UInt k) {
AKANTU_DEBUG_ASSERT((k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the slice " << k
<< " in a tensor3 of size (" << this->n[0] << ", "
<< this->n[1] << ", " << this->n[2] << ")");
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline const MatrixProxy<T> operator()(UInt k) const {
AKANTU_DEBUG_ASSERT((k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the slice " << k
<< " in a tensor3 of size (" << this->n[0] << ", "
<< this->n[1] << ", " << this->n[2] << ")");
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline MatrixProxy<T> operator[](UInt k) {
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline const MatrixProxy<T> operator[](UInt k) const {
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
};
/* -------------------------------------------------------------------------- */
// support operations for the creation of other vectors
/* -------------------------------------------------------------------------- */
template <typename T>
Vector<T> operator*(const T & scalar, const Vector<T> & a) {
Vector<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Vector<T> operator*(const Vector<T> & a, const T & scalar) {
Vector<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Vector<T> operator/(const Vector<T> & a, const T & scalar) {
Vector<T> r(a);
r /= scalar;
return r;
}
template <typename T>
Vector<T> operator+(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r += b;
return r;
}
template <typename T>
Vector<T> operator-(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r -= b;
return r;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Matrix<T> operator*(const T & scalar, const Matrix<T> & a) {
Matrix<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Matrix<T> operator*(const Matrix<T> & a, const T & scalar) {
Matrix<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Matrix<T> operator/(const Matrix<T> & a, const T & scalar) {
Matrix<T> r(a);
r /= scalar;
return r;
}
template <typename T>
Matrix<T> operator+(const Matrix<T> & a, const Matrix<T> & b) {
Matrix<T> r(a);
r += b;
return r;
}
template <typename T>
Matrix<T> operator-(const Matrix<T> & a, const Matrix<T> & b) {
Matrix<T> r(a);
r -= b;
return r;
}
__END_AKANTU__
#endif /* __AKANTU_AKA_TYPES_HH__ */