diff --git a/cmake/AkantuExtraCompilationProfiles.cmake b/cmake/AkantuExtraCompilationProfiles.cmake
index b2e5bc731..b96e88ec6 100644
--- a/cmake/AkantuExtraCompilationProfiles.cmake
+++ b/cmake/AkantuExtraCompilationProfiles.cmake
@@ -1,34 +1,34 @@
#Profiling
set(CMAKE_CXX_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2"
CACHE STRING "Flags used by the compiler during profiling builds")
set(CMAKE_C_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2"
CACHE STRING "Flags used by the compiler during profiling builds")
set(CMAKE_Fortran_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2"
CACHE STRING "Flags used by the compiler during profiling builds")
set(CMAKE_EXE_LINKER_FLAGS_PROFILING "-pg"
CACHE STRING "Flags used by the linker during profiling builds")
set(CMAKE_SHARED_LINKER_FLAGS_PROFILING "-pg"
CACHE STRING "Flags used by the linker during profiling builds")
# Sanitize the code
if ((CMAKE_CXX_COMPILER_ID STREQUAL "GNU" AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER "5.2") OR
CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
- set(CMAKE_CXX_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined"
+ set(CMAKE_CXX_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer"
CACHE STRING "Flags used by the compiler during sanitining builds")
- set(CMAKE_C_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined"
+ set(CMAKE_C_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer"
CACHE STRING "Flags used by the compiler during sanitining builds")
- set(CMAKE_Fortran_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined"
+ set(CMAKE_Fortran_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer"
CACHE STRING "Flags used by the compiler during sanitining builds")
- set(CMAKE_EXE_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined"
+ set(CMAKE_EXE_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer"
CACHE STRING "Flags used by the linker during sanitining builds")
- set(CMAKE_SHARED_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined"
+ set(CMAKE_SHARED_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer"
CACHE STRING "Flags used by the linker during sanitining builds")
mark_as_advanced(CMAKE_CXX_FLAGS_PROFILING CMAKE_C_FLAGS_PROFILING
CMAKE_Fortran_FLAGS_PROFILING CMAKE_EXE_LINKER_FLAGS_PROFILING
CMAKE_SHARED_LINKER_FLAGS_PROFILING CMAKE_SHARED_LINKER_FLAGS_SANITIZE
CMAKE_CXX_FLAGS_SANITIZE CMAKE_C_FLAGS_SANITIZE
CMAKE_Fortran_FLAGS_SANITIZE CMAKE_EXE_LINKER_FLAGS_SANITIZE
)
endif()
diff --git a/src/common/aka_array.cc b/src/common/aka_array.cc
index b31c3248d..d1601dbdc 100644
--- a/src/common/aka_array.cc
+++ b/src/common/aka_array.cc
@@ -1,117 +1,101 @@
/**
* @file aka_array.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Tue Feb 20 2018
*
* @brief Implementation of akantu::Array
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <memory>
#include <utility>
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "aka_common.hh"
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Functions ArrayBase */
/* -------------------------------------------------------------------------- */
-/* -------------------------------------------------------------------------- */
-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 <> UInt Array<Real>::find(const Real & elem) const {
AKANTU_DEBUG_IN();
Real epsilon = std::numeric_limits<Real>::epsilon();
auto it = std::find_if(begin(), end(), [&elem, &epsilon](auto && a) {
return std::abs(a - elem) <= epsilon;
});
AKANTU_DEBUG_OUT();
return (it != end()) ? end() - it : UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <>
Array<ElementType> & Array<ElementType>::operator*=(__attribute__((unused))
const ElementType & alpha) {
AKANTU_TO_IMPLEMENT();
return *this;
}
template <>
Array<ElementType> & Array<ElementType>::
operator-=(__attribute__((unused)) const Array<ElementType> & vect) {
AKANTU_TO_IMPLEMENT();
return *this;
}
template <>
Array<ElementType> & Array<ElementType>::
operator+=(__attribute__((unused)) const Array<ElementType> & vect) {
AKANTU_TO_IMPLEMENT();
return *this;
}
template <>
Array<char> & Array<char>::operator*=(__attribute__((unused))
const char & alpha) {
AKANTU_TO_IMPLEMENT();
return *this;
}
template <>
Array<char> & Array<char>::operator-=(__attribute__((unused))
const Array<char> & vect) {
AKANTU_TO_IMPLEMENT();
return *this;
}
template <>
Array<char> & Array<char>::operator+=(__attribute__((unused))
const Array<char> & vect) {
AKANTU_TO_IMPLEMENT();
return *this;
}
} // namespace akantu
diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh
index 2a36fd1e7..c29ee5fa5 100644
--- a/src/common/aka_array.hh
+++ b/src/common/aka_array.hh
@@ -1,381 +1,439 @@
/**
* @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: Tue Jan 16 2018
*
* @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-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_VECTOR_HH__
#define __AKANTU_VECTOR_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
#include <typeinfo>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace akantu {
/// class that afford to store vectors in static memory
class ArrayBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
explicit ArrayBase(ID id = "") : id(std::move(id)) {}
+ ArrayBase(const ArrayBase & other, const ID & id = "") {
+ this->id = (id == "") ? other.id : id;
+ }
ArrayBase(const ArrayBase & other) = default;
ArrayBase(ArrayBase && other) = default;
-
ArrayBase & operator=(const ArrayBase & other) = default;
- //ArrayBase & operator=(ArrayBase && other) = default;
+ // ArrayBase & operator=(ArrayBase && other) = default;
virtual ~ArrayBase() = default;
-
-
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// get the amount of space allocated in bytes
- inline UInt getMemorySize() const;
+ virtual UInt getMemorySize() const = 0;
/// 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;
+ virtual void printself(std::ostream & stream, int indent = 0) const = 0;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
- /// Get the real size allocated in memory
- AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt);
/// Get the Size of the Array
- UInt getSize() const __attribute__((deprecated)) { return size_; }
UInt size() const { return size_; }
/// 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 &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// id of the vector
- ID id;
-
- /// the size allocated
- UInt allocated_size{0};
+ ID id{""};
/// the size used
UInt size_{0};
/// number of components
UInt nb_component{1};
-
- /// size of the stored type
- UInt size_of_type{0};
};
/* -------------------------------------------------------------------------- */
namespace {
template <std::size_t dim, typename T> struct IteratorHelper {};
template <typename T> struct IteratorHelper<0, T> { using type = T; };
template <typename T> struct IteratorHelper<1, T> { using type = Vector<T>; };
template <typename T> struct IteratorHelper<2, T> { using type = Matrix<T>; };
template <typename T> struct IteratorHelper<3, T> {
using type = Tensor3<T>;
};
template <std::size_t dim, typename T>
using IteratorHelper_t = typename IteratorHelper<dim, T>::type;
} // namespace
/* -------------------------------------------------------------------------- */
-template <typename T, bool is_scal> class Array : public ArrayBase {
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
+/* Memory handling layer */
+/* -------------------------------------------------------------------------- */
+enum class ArrayAllocationType {
+ _default,
+ _pod,
+};
+
+template <typename T>
+struct ArrayAllocationTrait
+ : public std::conditional_t<
+ std::is_scalar<T>::value,
+ std::integral_constant<ArrayAllocationType,
+ ArrayAllocationType::_pod>,
+ std::integral_constant<ArrayAllocationType,
+ ArrayAllocationType::_default>> {};
+
+/* -------------------------------------------------------------------------- */
+template <typename T,
+ ArrayAllocationType allocation_trait = ArrayAllocationTrait<T>::value>
+class ArrayDataLayer : public ArrayBase {
public:
using value_type = T;
using reference = value_type &;
using pointer_type = value_type *;
using const_reference = const value_type &;
- inline ~Array() override;
+public:
+ virtual ~ArrayDataLayer() = default;
/// Allocation of a new vector
- explicit inline Array(UInt size = 0, UInt nb_component = 1,
- const ID & id = "");
+ explicit ArrayDataLayer(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 = "");
+ ArrayDataLayer(UInt size, UInt nb_component, const_reference value,
+ const ID & id = "");
+
+ /// Copy constructor (deep copy)
+ ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "");
+
+#ifndef SWIG
+ /// Copy constructor (deep copy)
+ explicit ArrayDataLayer(const std::vector<value_type> & vect);
+#endif
+
+ // copy operator
+ ArrayDataLayer & operator=(const ArrayDataLayer & other);
+
+ // move constructor
+ ArrayDataLayer(ArrayDataLayer && other);
+
+ // move assign
+ ArrayDataLayer & operator=(ArrayDataLayer && other);
+
+protected:
+ // allocate the memory
+ void allocate(UInt size, UInt nb_component);
+
+ // allocate and initialize the memory
+ void allocate(UInt size, UInt nb_component, const T & value);
+
+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[]);
+
+#ifndef SWIG
+ /// append a Vector or a Matrix
+ template <template <typename> class C,
+ typename = std::enable_if_t<is_tensor<C<T>>::value>>
+ inline void push_back(const C<T> & new_elem);
+#endif
+
+ /// changes the allocated size but not the size
+ void reserve(UInt size);
+
+ /// change the size of the Array
+ void resize(UInt size);
+
+ /// change the size of the Array and initialize the values
+ void resize(UInt size, const T & val);
+
+ /// get the amount of space allocated in bytes
+ inline UInt getMemorySize() const override final;
+
+ /// Get the real size allocated in memory
+ inline UInt getAllocatedSize() const;
+
+ /// give the address of the memory allocated for this vector
+ T * storage() const { return values; };
+
+protected:
+ /// allocation type agnostic data access
+ T * values{nullptr};
+
+ /// data storage
+ std::vector<T> data_storage;
+};
+
+/* -------------------------------------------------------------------------- */
+/* Actual Array */
+/* -------------------------------------------------------------------------- */
+template <typename T, bool is_scal> class Array : public ArrayDataLayer<T> {
+private:
+ using parent = ArrayDataLayer<T>;
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+public:
+ using value_type = typename parent::value_type;
+ using reference = typename parent::reference;
+ using pointer_type = typename parent::pointer_type;
+ using const_reference = typename parent::const_reference;
+
+ ~Array() override;
+
+ /// Allocation of a new vector
+ 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_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 = "");
+ Array(const Array & vect, const ID & id = "");
#ifndef SWIG
/// Copy constructor (deep copy)
- explicit Array(const std::vector<value_type> & vect);
+ explicit Array(const std::vector<T> & vect);
#endif
// copy operator
Array & operator=(const Array & other);
// move constructor
- Array(Array && other);
+ Array(Array && other) = default;
// move assign
- Array & operator=(Array && other);
+ Array & operator=(Array && other) = default;
#ifndef SWIG
/* ------------------------------------------------------------------------ */
/* Iterator */
/* ------------------------------------------------------------------------ */
/// \todo protected: does not compile with intel check why
public:
template <class R, class it, class IR = R,
bool is_tensor_ = is_tensor<R>::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
using scalar_iterator = iterator<T>;
/// const_iterator for Array of nb_component = 1
using const_scalar_iterator = const_iterator<T>;
/// iterator returning Vectors of size n on entries of Array with
/// nb_component = n
using vector_iterator = iterator<Vector<T>>;
/// const_iterator returning Vectors of n size on entries of Array with
/// nb_component = n
using const_vector_iterator = const_iterator<Vector<T>>;
/// iterator returning Matrices of size (m, n) on entries of Array with
/// nb_component = m*n
using matrix_iterator = iterator<Matrix<T>>;
/// const iterator returning Matrices of size (m, n) on entries of Array with
/// nb_component = m*n
using const_matrix_iterator = const_iterator<Matrix<T>>;
/// iterator returning Tensor3 of size (m, n, k) on entries of Array with
/// nb_component = m*n*k
using tensor3_iterator = iterator<Tensor3<T>>;
/// const iterator returning Tensor3 of size (m, n, k) on entries of Array
/// with nb_component = m*n*k
using const_tensor3_iterator = const_iterator<Tensor3<T>>;
/* ------------------------------------------------------------------------ */
template <typename... Ns> inline decltype(auto) begin(Ns &&... n);
template <typename... Ns> inline decltype(auto) end(Ns &&... n);
template <typename... Ns> inline decltype(auto) begin(Ns &&... n) const;
template <typename... Ns> inline decltype(auto) end(Ns &&... n) const;
template <typename... Ns> inline decltype(auto) begin_reinterpret(Ns &&... n);
template <typename... Ns> inline decltype(auto) end_reinterpret(Ns &&... n);
template <typename... Ns>
inline decltype(auto) begin_reinterpret(Ns &&... n) const;
template <typename... Ns>
inline decltype(auto) end_reinterpret(Ns &&... n) const;
#endif // SWIG
/* ------------------------------------------------------------------------ */
/* 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[]);
+ /// search elem in the vector, return the position of the first occurrence or
+ /// -1 if not found
+ UInt find(const_reference elem) const;
+
+ /// @see Array::find(const_reference elem) const
+ UInt find(T elem[]) const;
+
+ inline void push_back(const_reference value) { parent::push_back(value); }
#ifndef SWIG
/// append a Vector or a Matrix
template <template <typename> class C,
typename = std::enable_if_t<is_tensor<C<T>>::value>>
- 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);
+ inline void push_back(const C<T> & new_elem) {
+ parent::push_back(new_elem);
+ }
+
+ template <typename Ret> inline void push_back(const iterator<Ret> & it) {
+ push_back(*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);
-#endif
-
- /// changes the allocated size but not the size
- virtual void reserve(UInt size);
-
- /// 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
- UInt find(const_reference elem) const;
- /// @see Array::find(const_reference elem) const
- UInt find(T elem[]) const;
-
-#ifndef SWIG
/// @see Array::find(const_reference elem) const
template <template <typename> class C,
typename = std::enable_if_t<is_tensor<C<T>>::value>>
inline UInt find(const C<T> & elem);
#endif
- /// 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); }
+ inline void set(T t) {
+ std::fill_n(this->values, this->size_ * this->nb_component, t);
+ }
+
+ /// set all entries of the array to 0
+ inline void clear() { set(T()); }
#ifndef SWIG
/// set all tuples of the array to a given vector or matrix
/// @param vm Matrix or Vector to fill the array with
template <template <typename> class C,
typename = std::enable_if_t<is_tensor<C<T>>::value>>
inline void set(const C<T> & vm);
#endif
/// 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
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
};
/* -------------------------------------------------------------------------- */
/* 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;
}
} // 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 44dc7e3e3..d57fdd2f7 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1307 +1,1359 @@
/**
* @file aka_array_tmpl.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Jul 15 2010
* @date last modification: Tue Feb 20 2018
*
* @brief Inline functions of the classes Array<T> and ArrayBase
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
-/* Inline Functions Array<T> */
+/* Inline Functions Array<T> */
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <memory>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__
#define __AKANTU_AKA_ARRAY_TMPL_HH__
namespace akantu {
namespace debug {
struct ArrayException : public Exception {};
} // namespace debug
/* -------------------------------------------------------------------------- */
-template <class T, bool is_scal>
-inline T & Array<T, is_scal>::operator()(UInt i, UInt j) {
- AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
- AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
- "The value at position ["
- << i << "," << j << "] is out of range in array \""
- << id << "\"");
- return values[i * nb_component + j];
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(UInt size,
+ UInt nb_component,
+ const ID & id)
+ : ArrayBase(id) {
+ allocate(size, nb_component);
}
/* -------------------------------------------------------------------------- */
-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 <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(UInt size,
+ UInt nb_component,
+ const_reference value,
+ const ID & id)
+ : ArrayBase(id) {
+ allocate(size, nb_component, value);
}
-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 <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(const ArrayDataLayer & vect,
+ const ID & id)
+ : ArrayBase(vect, id) {
+ this->data_storage = vect.data_storage;
+ this->size_ = vect.size_;
+ this->nb_component = vect.nb_component;
+ this->values = this->data_storage.data();
}
+#ifndef SWIG
/* -------------------------------------------------------------------------- */
-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];
+template <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(
+ const std::vector<value_type> & vect) {
+ this->data_storage = vect;
+ this->size_ = vect.size();
+ this->nb_component = 1;
+ this->values = this->data_storage.data();
+}
+#endif
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait> & ArrayDataLayer<T, allocation_trait>::
+operator=(const ArrayDataLayer & other) {
+ if (this != &other) {
+ this->data_storage = other.data_storage;
+ this->nb_component = other.nb_component;
+ this->size_ = other.size_;
+ this->values = this->data_storage.data();
+ }
+ return *this;
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(ArrayDataLayer && other) =
+ default;
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+ArrayDataLayer<T, allocation_trait> & ArrayDataLayer<T, allocation_trait>::
+operator=(ArrayDataLayer && other) = default;
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+void ArrayDataLayer<T, allocation_trait>::allocate(UInt new_size,
+ UInt nb_component) {
+ this->nb_component = nb_component;
+ this->resize(new_size);
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+void ArrayDataLayer<T, allocation_trait>::allocate(UInt new_size,
+ UInt nb_component,
+ const T & val) {
+ this->nb_component = nb_component;
+ this->resize(new_size, val);
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+void ArrayDataLayer<T, allocation_trait>::resize(UInt new_size) {
+ this->data_storage.resize(new_size * this->nb_component);
+ this->values = this->data_storage.data();
+ this->size_ = new_size;
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+void ArrayDataLayer<T, allocation_trait>::resize(UInt new_size,
+ const T & value) {
+ this->data_storage.resize(new_size * this->nb_component, value);
+ this->values = this->data_storage.data();
+ this->size_ = new_size;
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+void ArrayDataLayer<T, allocation_trait>::reserve(UInt size) {
+ this->data_storage.reserve(size * this->nb_component);
+ this->values = this->data_storage.data();
}
/* -------------------------------------------------------------------------- */
/**
* 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);
+template <typename T, ArrayAllocationType allocation_trait>
+inline void ArrayDataLayer<T, allocation_trait>::push_back(const T & value) {
+ this->data_storage.push_back(value);
+ this->values = this->data_storage.data();
+ this->size_ += 1;
}
-/* -------------------------------------------------------------------------- */
-/**
- * append a tuple to the array
- * @param new_elem a C-array containing the values to be copied to the end of
- * the array */
-// template <class T, bool is_scal>
-// inline void Array<T, is_scal>::push_back(const T new_elem[]) {
-// UInt pos = size_;
-
-// resizeUnitialized(size_ + 1, false);
-
-// T * tmp = values + nb_component * pos;
-// std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
-// }
-
/* -------------------------------------------------------------------------- */
#ifndef SWIG
/**
* append a matrix or a vector to the array
* @param new_elem a reference to a Matrix<T> or Vector<T> */
-template <class T, bool is_scal>
+template <typename T, ArrayAllocationType allocation_trait>
template <template <typename> class C, typename>
-inline void Array<T, is_scal>::push_back(const C<T> & new_elem) {
+inline void
+ArrayDataLayer<T, allocation_trait>::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);
+ for (UInt i = 0; i < new_elem.size(); ++i) {
+ this->data_storage.push_back(new_elem[i]);
+ }
+ this->values = this->data_storage.data();
+ this->size_ += 1;
+}
+#endif
- T * tmp = values + nb_component * pos;
- std::uninitialized_copy(new_elem.storage(), new_elem.storage() + nb_component,
- tmp);
+/* -------------------------------------------------------------------------- */
+template <typename T, ArrayAllocationType allocation_trait>
+inline UInt ArrayDataLayer<T, allocation_trait>::getAllocatedSize() const {
+ return this->data_storage.capacity() / this->nb_component;
}
/* -------------------------------------------------------------------------- */
-/**
- * append a tuple to the array
- * @param it an iterator to the tuple to be copied to the end of the array */
+template <typename T, ArrayAllocationType allocation_trait>
+inline UInt ArrayDataLayer<T, allocation_trait>::getMemorySize() const {
+ return this->data_storage.capacity() * sizeof(T);
+}
+
+/* -------------------------------------------------------------------------- */
+
+/* -------------------------------------------------------------------------- */
+template <typename T>
+class ArrayDataLayer<T, ArrayAllocationType::_pod> : public ArrayBase {
+public:
+ using value_type = T;
+ using reference = value_type &;
+ using pointer_type = value_type *;
+ using const_reference = const value_type &;
+
+public:
+ virtual ~ArrayDataLayer() { free(values); }
+
+ /// Allocation of a new vector
+ ArrayDataLayer(UInt size = 0, UInt nb_component = 1, const ID & id = "")
+ : ArrayBase(id) {
+ allocate(size, nb_component);
+ }
+
+ /// Allocation of a new vector with a default value
+ ArrayDataLayer(UInt size, UInt nb_component, const_reference value,
+ const ID & id = "")
+ : ArrayBase(id) {
+ allocate(size, nb_component, value);
+ }
+
+ /// Copy constructor (deep copy)
+ ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "")
+ : ArrayBase(vect, id) {
+ allocate(vect.size(), vect.getNbComponent());
+ std::copy_n(vect.storage(), this->size_ * this->nb_component, values);
+ }
+
+#ifndef SWIG
+ /// Copy constructor (deep copy)
+ explicit ArrayDataLayer(const std::vector<value_type> & vect) {
+ allocate(vect.size(), 1);
+ std::copy_n(vect.data(), this->size_ * this->nb_component, values);
+ }
+#endif
+
+ // copy operator
+ inline ArrayDataLayer & operator=(const ArrayDataLayer & other) {
+ if (this != &other) {
+ allocate(other.size(), other.getNbComponent());
+ std::copy_n(other.storage(), this->size_ * this->nb_component, values);
+ }
+ return *this;
+ }
+
+ // move constructor
+ inline ArrayDataLayer(ArrayDataLayer && other) = default;
+
+ // move assign
+ inline ArrayDataLayer & operator=(ArrayDataLayer && other) = default;
+
+protected:
+ // allocate the memory
+ inline void allocate(UInt size, UInt nb_component) {
+ this->values =
+ reinterpret_cast<T *>(malloc(nb_component * size * sizeof(T)));
+ if (this->values == nullptr and size != 0) {
+ throw std::bad_alloc();
+ }
+ this->nb_component = nb_component;
+ this->allocated_size = this->size_ = size;
+ }
+
+ // allocate and initialize the memory
+ inline void allocate(UInt size, UInt nb_component, const T & value) {
+ allocate(size, nb_component);
+ std::fill_n(values, size * nb_component, value);
+ }
+
+public:
+ /// append a tuple of size nb_component containing value
+ inline void push_back(const_reference value) {
+ resize(this->size_ + 1, value);
+ }
+
+#ifndef SWIG
+ /// append a Vector or a Matrix
+ template <template <typename> class C,
+ typename = std::enable_if_t<is_tensor<C<T>>::value>>
+ inline void 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 << ").");
+ this->resize(this->size_ + 1);
+ std::copy_n(new_elem.storage(), new_elem.size(),
+ values + this->nb_component * (this->size_ - 1));
+ }
+#endif
+
+ /// changes the allocated size but not the size
+ void reserve(UInt size) {
+ UInt tmp_size = this->size_;
+ this->resize(size);
+ this->size_ = tmp_size;
+ }
+
+ /// change the size of the Array
+ void resize(UInt size) {
+ if (size == 0) {
+ free(values);
+ values = nullptr;
+ this->allocated_size = 0;
+ } else {
+ Int diff = size - allocated_size;
+ UInt size_to_allocate = (std::abs(diff) > AKANTU_MIN_ALLOCATION)
+ ? size
+ : (diff > 0)
+ ? allocated_size + AKANTU_MIN_ALLOCATION
+ : allocated_size;
+
+ auto * tmp_ptr = reinterpret_cast<T *>(realloc(
+ this->values, size_to_allocate * this->nb_component * sizeof(T)));
+ if (tmp_ptr == nullptr) {
+ throw std::bad_alloc();
+ }
+
+ this->values = tmp_ptr;
+ this->allocated_size = size_to_allocate;
+ }
+
+ this->size_ = size;
+ }
+
+ /// change the size of the Array and initialize the values
+ void resize(UInt size, const T & val) {
+ UInt tmp_size = this->size_;
+ this->resize(size);
+ if (size > tmp_size) {
+ std::fill_n(values + this->nb_component * tmp_size, (size - tmp_size),
+ val);
+ }
+ }
+
+ /// get the amount of space allocated in bytes
+ inline UInt getMemorySize() const override final {
+ return this->allocated_size * this->nb_component * sizeof(T);
+ }
+
+ /// Get the real size allocated in memory
+ inline UInt getAllocatedSize() const { return this->allocated_size; }
+
+ /// give the address of the memory allocated for this vector
+ T * storage() const { return values; };
+
+protected:
+ /// allocation type agnostic data access
+ T * values{nullptr};
+
+ UInt allocated_size{0};
+};
+/* -------------------------------------------------------------------------- */
+
+/* -------------------------------------------------------------------------- */
+/* template <class T> class AllocatorMalloc : public Allocator<T> {
+public:
+ T * allocate(UInt size, UInt nb_component) override final {
+ auto * ptr = reinterpret_cast<T *>(malloc(nb_component * size * sizeof(T)));
+
+ if (ptr == nullptr and size != 0) {
+ throw std::bad_alloc();
+ }
+ return ptr;
+ }
+
+ void deallocate(T * ptr, UInt size, UInt ,
+ UInt nb_component) override final {
+ if (ptr) {
+ if (not is_scalar<T>::value) {
+ for (UInt i = 0; i < size * nb_component; ++i) {
+ (ptr + i)->~T();
+ }
+ }
+ free(ptr);
+ }
+ }
+
+ std::tuple<T *, UInt> resize(UInt new_size, UInt size, UInt allocated_size,
+ UInt nb_component, T * ptr) override final {
+ UInt size_to_alloc = 0;
+
+ if (not is_scalar<T>::value and (new_size < size)) {
+ for (UInt i = new_size * nb_component; i < size * nb_component; ++i) {
+ (ptr + i)->~T();
+ }
+ }
+
+ // free some memory
+ if (new_size == 0) {
+ free(ptr);
+ return std::make_tuple(nullptr, 0);
+ }
+
+ if (new_size <= allocated_size) {
+ if (allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
+ size_to_alloc = new_size;
+ } else {
+ return std::make_tuple(ptr, allocated_size);
+ }
+ } else {
+ // allocate more memory
+ size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION)
+ ? allocated_size + AKANTU_MIN_ALLOCATION
+ : new_size;
+ }
+
+ auto * tmp_ptr = reinterpret_cast<T *>(
+ realloc(ptr, size_to_alloc * nb_component * sizeof(T)));
+ if (tmp_ptr == nullptr) {
+ throw std::bad_alloc();
+ }
+
+ return std::make_tuple(tmp_ptr, size_to_alloc);
+ }
+};
+*/
+
+/* -------------------------------------------------------------------------- */
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_;
+inline auto Array<T, is_scal>::operator()(UInt i, UInt j) -> reference {
+ AKANTU_DEBUG_ASSERT(this->size_ > 0,
+ "The array \"" << this->id << "\" is empty");
+ AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component),
+ "The value at position ["
+ << i << "," << j << "] is out of range in array \""
+ << this->id << "\"");
+ return this->values[i * this->nb_component + j];
+}
- resizeUnitialized(size_ + 1, false);
+/* -------------------------------------------------------------------------- */
+template <class T, bool is_scal>
+inline auto Array<T, is_scal>::operator()(UInt i, UInt j) const
+ -> const_reference {
+ AKANTU_DEBUG_ASSERT(this->size_ > 0,
+ "The array \"" << this->id << "\" is empty");
+ AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component),
+ "The value at position ["
+ << i << "," << j << "] is out of range in array \""
+ << this->id << "\"");
+ return this->values[i * this->nb_component + j];
+}
- T * tmp = values + nb_component * pos;
- T * new_elem = it.data();
- std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
+template <class T, bool is_scal>
+inline auto Array<T, is_scal>::operator[](UInt i) -> reference {
+ AKANTU_DEBUG_ASSERT(this->size_ > 0,
+ "The array \"" << this->id << "\" is empty");
+ AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component),
+ "The value at position ["
+ << i << "] is out of range in array \"" << this->id
+ << "\"");
+ return this->values[i];
}
-#endif
+
+/* -------------------------------------------------------------------------- */
+template <class T, bool is_scal>
+inline auto Array<T, is_scal>::operator[](UInt i) const -> const_reference {
+ AKANTU_DEBUG_ASSERT(this->size_ > 0,
+ "The array \"" << this->id << "\" is empty");
+ AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component),
+ "The value at position ["
+ << i << "] is out of range in array \"" << this->id
+ << "\"");
+ return this->values[i];
+}
+
/* -------------------------------------------------------------------------- */
/**
* 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];
+ AKANTU_DEBUG_ASSERT((this->size_ > 0), "The array is empty");
+ AKANTU_DEBUG_ASSERT((i < this->size_), "The element at position ["
+ << i << "] is out of range (" << i
+ << ">=" << this->size_ << ")");
+
+ if (i != (this->size_ - 1)) {
+ for (UInt j = 0; j < this->nb_component; ++j) {
+ this->values[i * this->nb_component + j] =
+ this->values[(this->size_ - 1) * this->nb_component + j];
}
}
- resize(size_ - 1);
+ this->resize(this->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),
+ AKANTU_DEBUG_ASSERT((this->size_ == vect.size_) &&
+ (this->nb_component == vect.nb_component),
"The too array don't have the same sizes");
- T * a = values;
+ T * a = this->values;
T * b = vect.storage();
- for (UInt i = 0; i < size_ * nb_component; ++i) {
+ for (UInt i = 0; i < this->size_ * this->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),
+ AKANTU_DEBUG_ASSERT((this->size_ == vect.size()) &&
+ (this->nb_component == vect.nb_component),
"The too array don't have the same sizes");
- T * a = values;
+ T * a = this->values;
T * b = vect.storage();
- for (UInt i = 0; i < size_ * nb_component; ++i) {
+ for (UInt i = 0; i < this->size_ * this->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
*/
#ifndef SWIG
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) {
+ T * a = this->values;
+ for (UInt i = 0; i < this->size_ * this->nb_component; ++i) {
*a++ *= alpha;
}
return *this;
}
#endif
/* -------------------------------------------------------------------------- */
/**
* 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;
+ bool equal = this->nb_component == array.nb_component &&
+ this->size_ == array.size_ && this->id == array.id;
if (!equal)
return false;
- if (values == array.storage())
+ if (this->values == array.storage())
return true;
else
- return std::equal(values, values + size_ * nb_component, array.storage());
+ return std::equal(this->values,
+ this->values + this->size_ * this->nb_component,
+ array.storage());
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
bool Array<T, is_scal>::operator!=(const Array<T, is_scal> & array) const {
return !operator==(array);
}
/* -------------------------------------------------------------------------- */
#ifndef SWIG
/**
* set all tuples of the array to a given vector or matrix
* @param vm Matrix or Vector to fill the array with
*/
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline void Array<T, is_scal>::set(const C<T> & vm) {
AKANTU_DEBUG_ASSERT(
- nb_component == vm.size(),
+ this->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);
+ for (T * it = this->values;
+ it < this->values + this->nb_component * this->size_;
+ it += this->nb_component) {
+ std::copy_n(vm.storage(), this->nb_component, it);
}
}
#endif
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::append(const Array<T> & other) {
AKANTU_DEBUG_ASSERT(
- nb_component == other.nb_component,
+ this->nb_component == other.nb_component,
"Cannot append an array with a different number of component");
UInt old_size = this->size_;
- this->resizeUnitialized(this->size_ + other.size(), false);
+ this->resize(this->size_ + other.size());
- T * tmp = values + nb_component * old_size;
- std::uninitialized_copy(other.storage(),
- other.storage() + other.size() * nb_component, tmp);
+ T * tmp = this->values + this->nb_component * old_size;
+ std::copy_n(other.storage(), other.size() * this->nb_component, tmp);
}
/* -------------------------------------------------------------------------- */
/* Functions Array<T, is_scal> */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const ID & id)
- : ArrayBase(id), values(nullptr) {
- AKANTU_DEBUG_IN();
- allocate(size, nb_component);
-
- if (!is_scal) {
- T val = T();
- std::uninitialized_fill(values, values + size * nb_component, val);
- }
+ : parent(size, nb_component, id) {}
- AKANTU_DEBUG_OUT();
-}
+template <>
+inline Array<std::string, false>::Array(UInt size, UInt nb_component,
+ const ID & id)
+ : parent(size, nb_component, "", id) {}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
-Array<T, is_scal>::Array(UInt size, UInt nb_component, const T def_values[],
+Array<T, is_scal>::Array(UInt size, UInt nb_component, const_reference value,
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();
-}
+ : parent(size, nb_component, value, id) {}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
-Array<T, is_scal>::Array(UInt size, UInt nb_component, const T & value,
- const ID & id)
- : ArrayBase(id), values(nullptr) {
- AKANTU_DEBUG_IN();
- allocate(size, nb_component);
-
- std::uninitialized_fill_n(values, size * nb_component, value);
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-template <class T, bool is_scal>
-Array<T, is_scal>::Array(const Array<T, is_scal> & vect, const ID & id)
- : ArrayBase(vect) {
- AKANTU_DEBUG_IN();
- this->id = (id == "") ? vect.id : id;
-
- // if (deep) {
- allocate(vect.size_, vect.nb_component);
- std::uninitialized_copy(vect.storage(), vect.storage() + size_ * nb_component,
- values);
- // } 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();
-}
+Array<T, is_scal>::Array(const Array & vect, const ID & id)
+ : parent(vect, id) {}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator=(const Array<T, is_scal> & other) {
AKANTU_DEBUG_WARNING("You are copying the array "
- << id << " are you sure it is on purpose");
+ << this->id << " are you sure it is on purpose");
if (&other == this)
return *this;
- allocate(other.size_, other.nb_component);
- std::uninitialized_copy(other.storage(),
- other.storage() + size_ * nb_component, values);
-
- return *this;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class T, bool is_scal>
-Array<T, is_scal>::Array(Array<T, is_scal> && other)
- : ArrayBase(other), values(std::exchange(other.values, nullptr)) {}
-
-/* -------------------------------------------------------------------------- */
-template <class T, bool is_scal>
-Array<T, is_scal> & Array<T, is_scal>::operator=(Array && other) {
- if (&other == this)
- return *this;
-
- ArrayBase::operator=(other);
- values = std::exchange(other.values, nullptr);
+ parent::operator=(other);
return *this;
}
/* -------------------------------------------------------------------------- */
#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();
-}
+Array<T, is_scal>::Array(const std::vector<T> & vect) : parent(vect) {}
#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 (not 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();
-}
-
-/* -------------------------------------------------------------------------- */
-template <class T, bool is_scal>
-void Array<T, is_scal>::reserve(UInt new_size) {
- UInt tmp_size = this->size_;
- resizeUnitialized(new_size, false);
- this->size_ = tmp_size;
-}
-
-/* -------------------------------------------------------------------------- */
-/**
- * change the size of the array and allocate or free memory if needed. If the
- * size increases, the new tuples are filled with zeros
- * @param new_size new number of tuples contained in the array */
-template <class T, bool is_scal> void Array<T, is_scal>::resize(UInt new_size) {
- resizeUnitialized(new_size, !is_scal);
-}
-
-/* -------------------------------------------------------------------------- */
-/**
- * change the size of the array and allocate or free memory if needed. If the
- * size increases, the new tuples are filled with zeros
- * @param new_size new number of tuples contained in the array */
-template <class T, bool is_scal>
-void Array<T, is_scal>::resize(UInt new_size, const T & val) {
- this->resizeUnitialized(new_size, true, val);
-}
-
-/* -------------------------------------------------------------------------- */
-/**
- * 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();
-
- values = static_cast<T *>(malloc(nb_component * size * sizeof(T)));
-
- if (values == nullptr and size != 0) {
- this->size_ = this->allocated_size = 0;
- AKANTU_EXCEPTION("Cannot allocate "
- << printMemorySize<T>(size * nb_component) << " (" << id
- << ")");
- } 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.
- * @param new_size new number of tuples contained in the array */
-template <class T, bool is_scal>
-void Array<T, is_scal>::resizeUnitialized(UInt new_size, bool fill,
- const T & val) {
- // AKANTU_DEBUG_IN();
- // free some memory
- if (new_size <= allocated_size) {
- if (!is_scal) {
- T * old_values = values;
- if (new_size < size_) {
- for (UInt i = new_size * nb_component; i < size_ * nb_component; ++i) {
- T * obj = old_values + i;
- obj->~T();
- }
- }
- }
-
- if (allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
- AKANTU_DEBUG(dblAccessory,
- "Freeing " << printMemorySize<T>((allocated_size - size_) *
- nb_component)
- << " (" << id << ")");
-
- // Normally there are no allocation problem when reducing an array
- if (new_size == 0) {
- free(values);
- values = nullptr;
- } else {
- auto * tmp_ptr = static_cast<T *>(
- realloc(values, new_size * nb_component * sizeof(T)));
-
- if (tmp_ptr == nullptr) {
- AKANTU_EXCEPTION("Cannot free data ("
- << id << ")"
- << " [current allocated size : " << allocated_size
- << " | "
- << "requested size : " << new_size << "]");
- }
- values = tmp_ptr;
- }
- allocated_size = new_size;
- }
- } else {
- // allocate more memory
- UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION)
- ? allocated_size + AKANTU_MIN_ALLOCATION
- : new_size;
-
- auto * tmp_ptr = static_cast<T *>(
- realloc(values, size_to_alloc * nb_component * sizeof(T)));
- AKANTU_DEBUG_ASSERT(
- tmp_ptr != nullptr,
- "Cannot allocate " << printMemorySize<T>(size_to_alloc * nb_component));
- if (tmp_ptr == nullptr) {
- AKANTU_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();
-}
+template <class T, bool is_scal> Array<T, is_scal>::~Array() = default;
/* -------------------------------------------------------------------------- */
/**
* search elem in the array, return the position of the first occurrence or
* -1 if not found
* @param elem the element to look for
* @return index of the first occurrence of elem or -1 if elem is not present
*/
template <class T, bool is_scal>
-UInt Array<T, is_scal>::find(const T & elem) const {
+UInt Array<T, is_scal>::find(const_reference elem) const {
AKANTU_DEBUG_IN();
auto begin = this->begin();
auto end = this->end();
auto it = std::find(begin, end, elem);
AKANTU_DEBUG_OUT();
return (it != end) ? it - begin : UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> UInt Array<T, is_scal>::find(T elem[]) const {
AKANTU_DEBUG_IN();
- T * it = values;
+ T * it = this->values;
UInt i = 0;
- for (; i < size_; ++i) {
+ for (; i < this->size_; ++i) {
if (*it == elem[0]) {
T * cit = it;
UInt c = 0;
- for (; (c < nb_component) && (*cit == elem[c]); ++c, ++cit)
+ for (; (c < this->nb_component) && (*cit == elem[c]); ++c, ++cit)
;
- if (c == nb_component) {
+ if (c == this->nb_component) {
AKANTU_DEBUG_OUT();
return i;
}
}
- it += nb_component;
+ it += this->nb_component;
}
return UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline UInt Array<T, is_scal>::find(const C<T> & elem) {
- AKANTU_DEBUG_ASSERT(elem.size() == nb_component,
+ AKANTU_DEBUG_ASSERT(elem.size() == this->nb_component,
"Cannot find an element with a wrong size ("
- << elem.size() << ") != " << nb_component);
+ << elem.size() << ") != " << this->nb_component);
return this->find(elem.storage());
}
/* -------------------------------------------------------------------------- */
/**
* copy the content of another array. This overwrites the current content.
* @param other Array to copy into this array. It has to have the same
* nb_component as this. If compiled in debug mode, an incorrect other will
* result in an exception being thrown. Optimised code may result in
* unpredicted behaviour.
*/
template <class T, bool is_scal>
void Array<T, is_scal>::copy(const Array<T, is_scal> & vect,
bool no_sanity_check) {
AKANTU_DEBUG_IN();
if (!no_sanity_check)
- if (vect.nb_component != nb_component)
+ if (vect.nb_component != this->nb_component)
AKANTU_ERROR("The two arrays do not have the same number of components");
- resize((vect.size_ * vect.nb_component) / nb_component);
+ this->resize((vect.size_ * vect.nb_component) / this->nb_component);
- T * tmp = values;
- std::uninitialized_copy(vect.storage(), vect.storage() + size_ * nb_component,
- tmp);
+ std::copy_n(vect.storage(), this->size_ * this->nb_component, this->values);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <bool is_scal> class ArrayPrintHelper {
public:
template <typename T>
static void print_content(const Array<T> & vect, std::ostream & stream,
int indent) {
if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << " + values : {";
for (UInt i = 0; i < vect.size(); ++i) {
stream << "{";
for (UInt j = 0; j < vect.getNbComponent(); ++j) {
stream << vect(i, j);
if (j != vect.getNbComponent() - 1)
stream << ", ";
}
stream << "}";
if (i != vect.size() - 1)
stream << ", ";
}
stream << "}" << std::endl;
}
}
};
template <> class ArrayPrintHelper<false> {
public:
template <typename T>
static void print_content(__attribute__((unused)) const Array<T> & vect,
__attribute__((unused)) std::ostream & stream,
__attribute__((unused)) int indent) {}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
std::streamsize prec = stream.precision();
std::ios_base::fmtflags ff = stream.flags();
stream.setf(std::ios_base::showbase);
stream.precision(2);
stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> ["
<< std::endl;
stream << space << " + id : " << this->id << std::endl;
stream << space << " + size : " << this->size_ << std::endl;
stream << space << " + nb_component : " << this->nb_component << std::endl;
- stream << space << " + allocated size : " << this->allocated_size
+ stream << space << " + allocated size : " << this->getAllocatedSize()
+ << std::endl;
+ stream << space
+ << " + memory size : " << printMemorySize<T>(this->getMemorySize())
<< 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; }
+inline void ArrayBase::empty() { this->size_ = 0; }
#ifndef SWIG
/* -------------------------------------------------------------------------- */
/* Iterators */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <class R, class daughter, class IR, bool is_tensor>
class Array<T, is_scal>::iterator_internal {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using proxy = typename R::proxy;
using const_proxy = const typename R::proxy;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
iterator_internal() = default;
iterator_internal(pointer_type data, UInt _offset)
: _offset(_offset), initial(data), ret(nullptr), ret_ptr(data) {
AKANTU_ERROR(
"The constructor should never be called it is just an ugly trick...");
}
iterator_internal(std::unique_ptr<internal_value_type> && wrapped)
: _offset(wrapped->size()), initial(wrapped->storage()),
ret(std::move(wrapped)), ret_ptr(ret->storage()) {}
iterator_internal(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
}
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
if (this->ret)
this->ret->shallowCopy(*it.ret);
else
this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
return *this;
}
UInt getCurrentIndex() {
return (this->ret_ptr - this->initial) / this->_offset;
};
inline reference operator*() {
ret->values = ret_ptr;
return *ret;
};
inline const_reference operator*() const {
ret->values = ret_ptr;
return *ret;
};
inline pointer operator->() {
ret->values = ret_ptr;
return ret.get();
};
inline daughter & operator++() {
ret_ptr += _offset;
return static_cast<daughter &>(*this);
};
inline daughter & operator--() {
ret_ptr -= _offset;
return static_cast<daughter &>(*this);
};
inline daughter & operator+=(const UInt n) {
ret_ptr += _offset * n;
return static_cast<daughter &>(*this);
}
inline daughter & operator-=(const UInt n) {
ret_ptr -= _offset * n;
return static_cast<daughter &>(*this);
}
inline proxy operator[](const UInt n) {
ret->values = ret_ptr + n * _offset;
return proxy(*ret);
}
inline const_proxy operator[](const UInt n) const {
ret->values = ret_ptr + n * _offset;
return const_proxy(*ret);
}
inline bool operator==(const iterator_internal & other) const {
return this->ret_ptr == other.ret_ptr;
}
inline bool operator!=(const iterator_internal & other) const {
return this->ret_ptr != other.ret_ptr;
}
inline bool operator<(const iterator_internal & other) const {
return this->ret_ptr < other.ret_ptr;
}
inline bool operator<=(const iterator_internal & other) const {
return this->ret_ptr <= other.ret_ptr;
}
inline bool operator>(const iterator_internal & other) const {
return this->ret_ptr > other.ret_ptr;
}
inline bool operator>=(const iterator_internal & other) const {
return this->ret_ptr >= other.ret_ptr;
}
inline daughter operator+(difference_type n) {
daughter tmp(static_cast<daughter &>(*this));
tmp += n;
return tmp;
}
inline daughter operator-(difference_type n) {
daughter tmp(static_cast<daughter &>(*this));
tmp -= n;
return tmp;
}
inline difference_type operator-(const iterator_internal & b) {
return (this->ret_ptr - b.ret_ptr) / _offset;
}
inline pointer_type data() const { return ret_ptr; }
inline difference_type offset() const { return _offset; }
protected:
UInt _offset{0};
pointer_type initial{nullptr};
std::unique_ptr<internal_value_type> ret{nullptr};
pointer_type ret_ptr{nullptr};
};
/* -------------------------------------------------------------------------- */
/**
* Specialization for scalar types
*/
template <class T, bool is_scal>
template <class R, class daughter, class IR>
class Array<T, is_scal>::iterator_internal<R, daughter, IR, false> {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
iterator_internal(pointer data = nullptr) : ret(data), initial(data){};
iterator_internal(const iterator_internal & it) = default;
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) = default;
UInt getCurrentIndex() { return (this->ret - this->initial); };
inline reference operator*() { return *ret; };
inline const_reference operator*() const { return *ret; };
inline pointer operator->() { return ret; };
inline daughter & operator++() {
++ret;
return static_cast<daughter &>(*this);
};
inline daughter & operator--() {
--ret;
return static_cast<daughter &>(*this);
};
inline daughter & operator+=(const UInt n) {
ret += n;
return static_cast<daughter &>(*this);
}
inline daughter & operator-=(const UInt n) {
ret -= n;
return static_cast<daughter &>(*this);
}
inline reference operator[](const UInt n) { return ret[n]; }
inline bool operator==(const iterator_internal & other) const {
return ret == other.ret;
}
inline bool operator!=(const iterator_internal & other) const {
return ret != other.ret;
}
inline bool operator<(const iterator_internal & other) const {
return ret < other.ret;
}
inline bool operator<=(const iterator_internal & other) const {
return ret <= other.ret;
}
inline bool operator>(const iterator_internal & other) const {
return ret > other.ret;
}
inline bool operator>=(const iterator_internal & other) const {
return ret >= other.ret;
}
inline daughter operator-(difference_type n) { return daughter(ret - n); }
inline daughter operator+(difference_type n) { return daughter(ret + n); }
inline difference_type operator-(const iterator_internal & b) {
return ret - b.ret;
}
inline pointer data() const { return ret; }
protected:
pointer ret{nullptr};
pointer initial{nullptr};
};
/* -------------------------------------------------------------------------- */
/* Begin/End functions implementation */
/* -------------------------------------------------------------------------- */
namespace detail {
template <class Tuple, size_t... Is>
constexpr auto take_front_impl(Tuple && t, std::index_sequence<Is...>) {
return std::make_tuple(std::get<Is>(std::forward<Tuple>(t))...);
}
template <size_t N, class Tuple> constexpr auto take_front(Tuple && t) {
return take_front_impl(std::forward<Tuple>(t),
std::make_index_sequence<N>{});
}
template <typename... V> constexpr auto product_all(V &&... v) {
std::common_type_t<int, V...> result = 1;
(void)std::initializer_list<int>{(result *= v, 0)...};
return result;
}
template <typename... T> std::string to_string_all(T &&... t) {
if (sizeof...(T) == 0)
return "";
std::stringstream ss;
bool noComma = true;
ss << "(";
(void)std::initializer_list<bool>{
(ss << (noComma ? "" : ", ") << t, noComma = false)...};
ss << ")";
return ss.str();
}
template <std::size_t N> struct InstantiationHelper {
template <typename type, typename T, typename... Ns>
static auto instantiate(T && data, Ns... ns) {
return std::make_unique<type>(data, ns...);
}
};
template <> struct InstantiationHelper<0> {
template <typename type, typename T> static auto instantiate(T && data) {
return data;
}
};
template <typename Arr, typename T, typename... Ns>
decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) {
using type = IteratorHelper_t<sizeof...(Ns) - 1, T>;
using array_type = std::decay_t<Arr>;
using iterator =
std::conditional_t<std::is_const<std::remove_reference_t<Arr>>::value,
typename array_type::template const_iterator<type>,
typename array_type::template iterator<type>>;
static_assert(sizeof...(Ns), "You should provide a least one size");
if (array.getNbComponent() * array.size() !=
product_all(std::forward<Ns>(ns)...)) {
AKANTU_CUSTOM_EXCEPTION_INFO(
debug::ArrayException(),
"The iterator on "
<< debug::demangle(typeid(Arr).name())
<< to_string_all(array.size(), array.getNbComponent())
<< "is not compatible with the type "
<< debug::demangle(typeid(type).name()) << to_string_all(ns...));
}
auto && wrapped = aka::apply(
[&](auto... n) {
return InstantiationHelper<sizeof...(n)>::template instantiate<type>(
data, n...);
},
take_front<sizeof...(Ns) - 1>(std::make_tuple(ns...)));
return iterator(std::move(wrapped));
}
} // namespace detail
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) {
- return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
+ return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...,
+ this->size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) {
- return detail::get_iterator(*this, values + nb_component * size_,
- std::forward<Ns>(ns)..., size_);
+ return detail::get_iterator(*this,
+ this->values + this->nb_component * this->size_,
+ std::forward<Ns>(ns)..., this->size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) const {
- return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
+ return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...,
+ this->size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) const {
- return detail::get_iterator(*this, values + nb_component * size_,
- std::forward<Ns>(ns)..., size_);
+ return detail::get_iterator(*this,
+ this->values + this->nb_component * this->size_,
+ std::forward<Ns>(ns)..., this->size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) {
- return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
+ return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) {
return detail::get_iterator(
- *this, values + detail::product_all(std::forward<Ns>(ns)...),
+ *this, this->values + detail::product_all(std::forward<Ns>(ns)...),
std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) const {
- return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
+ return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) const {
return detail::get_iterator(
- *this, values + detail::product_all(std::forward<Ns>(ns)...),
+ *this, this->values + detail::product_all(std::forward<Ns>(ns)...),
std::forward<Ns>(ns)...);
}
/* -------------------------------------------------------------------------- */
/* Views */
/* -------------------------------------------------------------------------- */
namespace detail {
template <typename Array, typename... Ns> class ArrayView {
using tuple = std::tuple<Ns...>;
public:
ArrayView(Array && array, Ns... ns)
: array(std::forward<Array>(array)), sizes(std::move(ns)...) {}
ArrayView(ArrayView && array_view) = default;
ArrayView & operator=(const ArrayView & array_view) = default;
ArrayView & operator=(ArrayView && array_view) = default;
decltype(auto) begin() {
return aka::apply(
[&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
}
decltype(auto) begin() const {
return aka::apply(
[&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
}
decltype(auto) end() {
return aka::apply(
[&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
}
decltype(auto) end() const {
return aka::apply(
[&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
}
decltype(auto) size() const {
return std::get<std::tuple_size<tuple>::value - 1>(sizes);
}
decltype(auto) dims() const { return std::tuple_size<tuple>::value - 1; }
private:
Array array;
tuple sizes;
};
} // namespace detail
/* -------------------------------------------------------------------------- */
template <typename Array, typename... Ns>
decltype(auto) make_view(Array && array, Ns... ns) {
static_assert(aka::conjunction<std::is_integral<std::decay_t<Ns>>...>::value,
"Ns should be integral types");
auto size = std::forward<decltype(array)>(array).size() *
std::forward<decltype(array)>(array).getNbComponent() /
detail::product_all(ns...);
return detail::ArrayView<Array, std::common_type_t<size_t, Ns>...,
std::common_type_t<size_t, decltype(size)>>(
std::forward<Array>(array), std::move(ns)..., size);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::const_iterator
: public iterator_internal<const R, Array<T, is_scal>::const_iterator<R>,
R> {
public:
using parent = iterator_internal<const R, const_iterator, R>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
const_iterator() : parent(){};
// const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}
// const_iterator(pointer warped) : parent(warped) {}
// const_iterator(const parent & it) : parent(it) {}
const_iterator(const const_iterator & it) = default;
const_iterator(const_iterator && it) = default;
template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
const_iterator(P * data) : parent(data) {}
- template <typename UP_P,
- typename =
- std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
+ template <typename UP_P, typename = std::enable_if_t<
+ is_tensor<typename UP_P::element_type>::value>>
const_iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
const_iterator & operator=(const const_iterator & it) = default;
};
template <class T, class R, bool is_tensor_ = is_tensor<R>::value>
struct ConstConverterIteratorHelper {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
return const_iterator(std::unique_ptr<R>(new R(*it, false)));
}
};
template <class T, class R> struct ConstConverterIteratorHelper<T, R, false> {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
return const_iterator(it.data());
}
};
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::iterator
: public iterator_internal<R, Array<T, is_scal>::iterator<R>> {
public:
using parent = iterator_internal<R, iterator>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
iterator() : parent(){};
iterator(const iterator & it) = default;
iterator(iterator && it) = default;
template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
iterator(P * data) : parent(data) {}
- template <typename UP_P,
- typename =
- std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
+ template <typename UP_P, typename = std::enable_if_t<
+ is_tensor<typename UP_P::element_type>::value>>
iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
iterator & operator=(const iterator & it) = default;
operator const_iterator<R>() {
return ConstConverterIteratorHelper<T, R>::convert(*this);
}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
inline typename Array<T, is_scal>::template iterator<R>
Array<T, is_scal>::erase(const iterator<R> & it) {
T * curr = it.data();
- UInt pos = (curr - values) / nb_component;
+ UInt pos = (curr - this->values) / this->nb_component;
erase(pos);
iterator<R> rit = it;
return --rit;
}
#endif
} // namespace akantu
#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */
diff --git a/src/common/aka_bbox.hh b/src/common/aka_bbox.hh
index b0bc5710f..9bc4183e8 100644
--- a/src/common/aka_bbox.hh
+++ b/src/common/aka_bbox.hh
@@ -1,261 +1,265 @@
/**
* @file aka_bbox.hh
*
* @author Nicolas Richart
*
* @date creation Mon Feb 12 2018
*
* @brief A simple bounding box class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_iterators.hh"
#include "aka_types.hh"
#include "communicator.hh"
/* -------------------------------------------------------------------------- */
+#include <map>
+/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_BBOX_HH__
#define __AKANTU_AKA_BBOX_HH__
namespace akantu {
class BBox {
public:
+ BBox() = default;
+
BBox(UInt spatial_dimension)
: dim(spatial_dimension),
lower_bounds(spatial_dimension, std::numeric_limits<Real>::max()),
upper_bounds(spatial_dimension, -std::numeric_limits<Real>::max()) {}
BBox(const BBox & other)
: dim(other.dim), empty{false}, lower_bounds(other.lower_bounds),
upper_bounds(other.upper_bounds) {}
BBox & operator=(const BBox & other) {
if (this != &other) {
this->dim = other.dim;
this->lower_bounds = other.lower_bounds;
this->upper_bounds = other.upper_bounds;
this->empty = other.empty;
}
return *this;
}
inline BBox & operator+=(const Vector<Real> & position) {
AKANTU_DEBUG_ASSERT(
this->dim == position.size(),
"You are adding a point of a wrong dimension to the bounding box");
this->empty = false;
for (auto s : arange(dim)) {
lower_bounds(s) = std::min(lower_bounds(s), position(s));
upper_bounds(s) = std::max(upper_bounds(s), position(s));
}
return *this;
}
/* ------------------------------------------------------------------------ */
inline bool intersects(const BBox & other,
const SpatialDirection & direction) const {
AKANTU_DEBUG_ASSERT(
this->dim == other.dim,
"You are intersecting bounding boxes of different dimensions");
return Math::intersects(lower_bounds(direction), upper_bounds(direction),
other.lower_bounds(direction),
other.upper_bounds(direction));
}
inline bool intersects(const BBox & other) const {
if (this->empty or other.empty)
return false;
bool intersects_ = true;
for (auto s : arange(this->dim)) {
intersects_ &= this->intersects(other, SpatialDirection(s));
}
return intersects_;
}
/* ------------------------------------------------------------------------ */
inline BBox intersection(const BBox & other) const {
AKANTU_DEBUG_ASSERT(
this->dim == other.dim,
"You are intersecting bounding boxes of different dimensions");
BBox intersection_(this->dim);
intersection_.empty = not this->intersects(other);
if (intersection_.empty)
return intersection_;
for (auto s : arange(this->dim)) {
// is lower point in range ?
bool point1 = Math::is_in_range(other.lower_bounds(s), lower_bounds(s),
upper_bounds(s));
// is upper point in range ?
bool point2 = Math::is_in_range(other.upper_bounds(s), lower_bounds(s),
upper_bounds(s));
if (point1 and not point2) {
// |-----------| this (i)
// |-----------| other(i)
// 1 2
intersection_.lower_bounds(s) = other.lower_bounds(s);
intersection_.upper_bounds(s) = upper_bounds(s);
} else if (point1 && point2) {
// |-----------------| this (i)
// |-----------| other(i)
// 1 2
intersection_.lower_bounds(s) = other.lower_bounds(s);
intersection_.upper_bounds(s) = other.upper_bounds(s);
} else if (!point1 && point2) {
// |-----------| this (i)
// |-----------| other(i)
// 1 2
intersection_.lower_bounds(s) = this->lower_bounds(s);
intersection_.upper_bounds(s) = other.upper_bounds(s);
} else {
// |-----------| this (i)
// |-----------------| other(i)
// 1 2
intersection_.lower_bounds(s) = this->lower_bounds(s);
intersection_.upper_bounds(s) = this->upper_bounds(s);
}
}
return intersection_;
}
/* ------------------------------------------------------------------------ */
inline bool contains(const Vector<Real> & point) const {
return (point >= lower_bounds) and (point <= upper_bounds);
}
/* ------------------------------------------------------------------------ */
inline void reset() {
lower_bounds.set(std::numeric_limits<Real>::max());
upper_bounds.set(-std::numeric_limits<Real>::max());
}
/* ------------------------------------------------------------------------ */
const Vector<Real> & getLowerBounds() const { return lower_bounds; }
const Vector<Real> & getUpperBounds() const { return upper_bounds; }
Vector<Real> & getLowerBounds() { return lower_bounds; }
Vector<Real> & getUpperBounds() { return upper_bounds; }
/* ------------------------------------------------------------------------ */
inline Real size(const SpatialDirection & direction) const {
return upper_bounds(direction) - lower_bounds(direction);
}
Vector<Real> size() const {
Vector<Real> size_(dim);
for (auto s : arange(this->dim)) {
size_(s) = this->size(SpatialDirection(s));
}
return size_;
}
inline operator bool() const { return not empty; }
/* ------------------------------------------------------------------------ */
BBox allSum(const Communicator & communicator) const {
Matrix<Real> reduce_bounds(dim, 2);
Vector<Real>(reduce_bounds(0)) = lower_bounds;
Vector<Real>(reduce_bounds(1)) = Real(-1.) * upper_bounds;
communicator.allReduce(reduce_bounds, SynchronizerOperation::_min);
BBox global(dim);
global.lower_bounds = Vector<Real>(reduce_bounds(0));
global.upper_bounds = Real(-1.) * Vector<Real>(reduce_bounds(1));
global.empty = false;
return global;
}
std::vector<BBox> allGather(const Communicator & communicator) const {
auto prank = communicator.whoAmI();
auto nb_proc = communicator.getNbProc();
Array<Real> bboxes_data(nb_proc, dim * 2 + 1);
auto * base = bboxes_data.storage() + prank * (2 * dim + 1);
Vector<Real>(base + dim * 0, dim) = lower_bounds;
Vector<Real>(base + dim * 1, dim) = upper_bounds;
base[dim * 2] = empty ? 1. : 0.; // ugly trick
communicator.allGather(bboxes_data);
std::vector<BBox> bboxes;
bboxes.reserve(nb_proc);
for (auto p : arange(nb_proc)) {
bboxes.emplace_back(dim);
auto & bbox = bboxes.back();
auto * base = bboxes_data.storage() + p * (2 * dim + 1);
bbox.lower_bounds = Vector<Real>(base + dim * 0, dim);
bbox.upper_bounds = Vector<Real>(base + dim * 1, dim);
bbox.empty = base[dim * 2] == 1. ? true : false;
}
return bboxes;
}
std::map<UInt, BBox> intersection(const BBox & other,
const Communicator & communicator) {
// todo: change for a custom reduction algorithm
auto other_bboxes = other.allGather(communicator);
std::map<UInt, BBox> intersections;
- for (auto & bbox : enumerate(other_bboxes)) {
+ for (const auto & bbox : enumerate(other_bboxes)) {
auto && tmp = this->intersection(std::get<1>(bbox));
if (tmp) {
intersections[std::get<0>(bbox)] = tmp;
}
}
return intersections;
}
void printself(std::ostream & stream) const {
stream << "BBox[";
if (not empty) {
stream << lower_bounds << " - " << upper_bounds;
}
stream << "]";
}
protected:
- UInt dim;
+ UInt dim{0};
bool empty{true};
Vector<Real> lower_bounds;
Vector<Real> upper_bounds;
};
inline std::ostream & operator<<(std::ostream & stream, const BBox & bbox) {
bbox.printself(stream);
return stream;
}
} // akantu
#endif /* __AKANTU_AKA_BBOX_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc
index 65cd472b9..6c576cf54 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc
@@ -1,551 +1,562 @@
/**
* @file material_cohesive.cc
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Seyedeh Mohadeseh Taheri Mousavi <mohadeseh.taherimousavi@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Feb 22 2012
* @date last modification: Mon Feb 19 2018
*
* @brief Specialization of the material class for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive.hh"
#include "aka_random_generator.hh"
#include "dof_synchronizer.hh"
#include "fe_engine_template.hh"
#include "integrator_gauss.hh"
#include "shape_cohesive.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id)
: Material(model, id),
facet_filter("facet_filter", id, this->getMemoryID()),
fem_cohesive(
model.getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine")),
reversible_energy("reversible_energy", *this),
total_energy("total_energy", *this), opening("opening", *this),
tractions("tractions", *this),
contact_tractions("contact_tractions", *this),
contact_opening("contact_opening", *this), delta_max("delta max", *this),
use_previous_delta_max(false), use_previous_opening(false),
damage("damage", *this), sigma_c("sigma_c", *this),
normal(0, spatial_dimension, "normal") {
AKANTU_DEBUG_IN();
this->model = dynamic_cast<SolidMechanicsModelCohesive *>(&model);
this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable,
"Critical stress");
this->registerParam("delta_c", delta_c, Real(0.),
_pat_parsable | _pat_readable, "Critical displacement");
this->element_filter.initialize(this->model->getMesh(),
_spatial_dimension = spatial_dimension,
_element_kind = _ek_cohesive);
// this->model->getMesh().initElementTypeMapArray(
// this->element_filter, 1, spatial_dimension, false, _ek_cohesive);
if (this->model->getIsExtrinsic())
this->facet_filter.initialize(this->model->getMeshFacets(),
_spatial_dimension = spatial_dimension - 1,
_element_kind = _ek_regular);
// this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1,
// spatial_dimension -
// 1);
this->reversible_energy.initialize(1);
this->total_energy.initialize(1);
this->tractions.initialize(spatial_dimension);
this->tractions.initializeHistory();
this->contact_tractions.initialize(spatial_dimension);
this->contact_opening.initialize(spatial_dimension);
this->opening.initialize(spatial_dimension);
this->opening.initializeHistory();
this->delta_max.initialize(1);
this->damage.initialize(1);
if (this->model->getIsExtrinsic())
this->sigma_c.initialize(1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MaterialCohesive::~MaterialCohesive() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
if (this->use_previous_delta_max)
this->delta_max.initializeHistory();
if (this->use_previous_opening)
this->opening.initializeHistory();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::assembleInternalForces(GhostType ghost_type) {
AKANTU_DEBUG_IN();
#if defined(AKANTU_DEBUG_TOOLS)
debug::element_manager.printData(debug::_dm_material_cohesive,
"Cohesive Tractions", tractions);
#endif
auto & internal_force = const_cast<Array<Real> &>(model->getInternalForce());
for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type,
_ek_cohesive)) {
auto & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (nb_element == 0)
continue;
const auto & shapes = fem_cohesive.getShapes(type, ghost_type);
auto & traction = tractions(type, ghost_type);
auto & contact_traction = contact_tractions(type, ghost_type);
UInt size_of_shapes = shapes.getNbComponent();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points =
fem_cohesive.getNbIntegrationPoints(type, ghost_type);
/// compute @f$t_i N_a@f$
Array<Real> * traction_cpy = new Array<Real>(
nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes);
auto traction_it = traction.begin(spatial_dimension, 1);
auto contact_traction_it = contact_traction.begin(spatial_dimension, 1);
auto shapes_filtered_begin = shapes.begin(1, size_of_shapes);
auto traction_cpy_it =
traction_cpy->begin(spatial_dimension, size_of_shapes);
Matrix<Real> traction_tmp(spatial_dimension, 1);
for (UInt el = 0; el < nb_element; ++el) {
UInt current_quad = elem_filter(el) * nb_quadrature_points;
for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it,
++contact_traction_it, ++current_quad, ++traction_cpy_it) {
const Matrix<Real> & shapes_filtered =
shapes_filtered_begin[current_quad];
traction_tmp.copy(*traction_it);
traction_tmp += *contact_traction_it;
traction_cpy_it->mul<false, false>(traction_tmp, shapes_filtered);
}
}
/**
* compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t
* \mathbf{t}_q \overline w_q J_q@f$
*/
- Array<Real> * int_t_N = new Array<Real>(
+ Array<Real> * partial_int_t_N = new Array<Real>(
nb_element, spatial_dimension * size_of_shapes, "int_t_N");
- fem_cohesive.integrate(*traction_cpy, *int_t_N,
+ fem_cohesive.integrate(*traction_cpy, *partial_int_t_N,
spatial_dimension * size_of_shapes, type, ghost_type,
elem_filter);
delete traction_cpy;
- int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent());
+ Array<Real> * int_t_N = new Array<Real>(
+ nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N");
Real * int_t_N_val = int_t_N->storage();
+ Real * partial_int_t_N_val = partial_int_t_N->storage();
for (UInt el = 0; el < nb_element; ++el) {
+ std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension,
+ int_t_N_val);
+ std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension,
+ int_t_N_val + size_of_shapes * spatial_dimension);
+
for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n)
int_t_N_val[n] *= -1.;
+
int_t_N_val += nb_nodes_per_element * spatial_dimension;
+ partial_int_t_N_val += size_of_shapes * spatial_dimension;
}
+ delete partial_int_t_N;
+
/// assemble
model->getDOFManager().assembleElementalArrayLocalArray(
*int_t_N, internal_force, type, ghost_type, 1, elem_filter);
delete int_t_N;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) {
AKANTU_DEBUG_IN();
for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type,
_ek_cohesive)) {
UInt nb_quadrature_points =
fem_cohesive.getNbIntegrationPoints(type, ghost_type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Array<Real> & shapes = fem_cohesive.getShapes(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (!nb_element)
continue;
UInt size_of_shapes = shapes.getNbComponent();
Array<Real> * shapes_filtered = new Array<Real>(
nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes");
Real * shapes_filtered_val = shapes_filtered->storage();
UInt * elem_filter_val = elem_filter.storage();
for (UInt el = 0; el < nb_element; ++el) {
- auto shapes_val =
- shapes.storage() +
- elem_filter_val[el] * size_of_shapes * nb_quadrature_points;
+ auto shapes_val = shapes.storage() + elem_filter_val[el] *
+ size_of_shapes *
+ nb_quadrature_points;
memcpy(shapes_filtered_val, shapes_val,
size_of_shapes * nb_quadrature_points * sizeof(Real));
shapes_filtered_val += size_of_shapes * nb_quadrature_points;
}
Matrix<Real> A(spatial_dimension * size_of_shapes,
spatial_dimension * nb_nodes_per_element);
for (UInt i = 0; i < spatial_dimension * size_of_shapes; ++i) {
A(i, i) = 1;
A(i, i + spatial_dimension * size_of_shapes) = -1;
}
/// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}}
/// @f$
Array<Real> * tangent_stiffness_matrix = new Array<Real>(
nb_element * nb_quadrature_points,
spatial_dimension * spatial_dimension, "tangent_stiffness_matrix");
// Array<Real> * normal = new Array<Real>(nb_element *
// nb_quadrature_points, spatial_dimension, "normal");
normal.resize(nb_quadrature_points);
computeNormal(model->getCurrentPosition(), normal, type, ghost_type);
/// compute openings @f$\mathbf{\delta}@f$
// computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
// ghost_type);
tangent_stiffness_matrix->clear();
computeTangentTraction(type, *tangent_stiffness_matrix, normal, ghost_type);
// delete normal;
UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element *
spatial_dimension * nb_nodes_per_element;
Array<Real> * at_nt_d_n_a = new Array<Real>(
nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A");
Array<Real>::iterator<Vector<Real>> shapes_filt_it =
shapes_filtered->begin(size_of_shapes);
Array<Real>::matrix_iterator D_it =
tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator At_Nt_D_N_A_it =
at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element,
spatial_dimension * nb_nodes_per_element);
Array<Real>::matrix_iterator At_Nt_D_N_A_end =
at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element,
spatial_dimension * nb_nodes_per_element);
Matrix<Real> N(spatial_dimension, spatial_dimension * size_of_shapes);
Matrix<Real> N_A(spatial_dimension,
spatial_dimension * nb_nodes_per_element);
Matrix<Real> D_N_A(spatial_dimension,
spatial_dimension * nb_nodes_per_element);
for (; At_Nt_D_N_A_it != At_Nt_D_N_A_end;
++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) {
N.clear();
/**
* store the shapes in voigt notations matrix @f$\mathbf{N} =
* \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\
* 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$
**/
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt n = 0; n < size_of_shapes; ++n)
N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n);
/**
* compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T
* \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}}
*{\partial{\delta}}
* \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$
**/
N_A.mul<false, false>(N, A);
D_N_A.mul<false, false>(*D_it, N_A);
(*At_Nt_D_N_A_it).mul<true, false>(D_N_A, N_A);
}
delete tangent_stiffness_matrix;
delete shapes_filtered;
Array<Real> * K_e = new Array<Real>(nb_element, size_at_nt_d_n_a, "K_e");
fem_cohesive.integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, type,
ghost_type, elem_filter);
delete at_nt_d_n_a;
model->getDOFManager().assembleElementalMatricesToMatrix(
"K", "displacement", *K_e, type, ghost_type, _unsymmetric, elem_filter);
delete K_e;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- *
* Compute traction from displacements
*
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void MaterialCohesive::computeTraction(GhostType ghost_type) {
AKANTU_DEBUG_IN();
#if defined(AKANTU_DEBUG_TOOLS)
debug::element_manager.printData(debug::_dm_material_cohesive,
"Cohesive Openings", opening);
#endif
for (auto & type : element_filter.elementTypes(spatial_dimension, ghost_type,
_ek_cohesive)) {
Array<UInt> & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (nb_element == 0)
continue;
UInt nb_quadrature_points =
nb_element * fem_cohesive.getNbIntegrationPoints(type, ghost_type);
normal.resize(nb_quadrature_points);
/// compute normals @f$\mathbf{n}@f$
computeNormal(model->getCurrentPosition(), normal, type, ghost_type);
/// compute openings @f$\mathbf{\delta}@f$
computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
ghost_type);
/// compute traction @f$\mathbf{t}@f$
computeTraction(normal, type, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::computeNormal(const Array<Real> & position,
Array<Real> & normal, ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & fem_cohesive =
this->model->getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine");
#define COMPUTE_NORMAL(type) \
fem_cohesive.getShapeFunctions() \
.computeNormalsOnIntegrationPoints<type, CohesiveReduceFunctionMean>( \
position, normal, ghost_type, element_filter(type, ghost_type));
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL);
#undef COMPUTE_NORMAL
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::computeOpening(const Array<Real> & displacement,
Array<Real> & opening, ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & fem_cohesive =
this->model->getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine");
#define COMPUTE_OPENING(type) \
fem_cohesive.getShapeFunctions() \
.interpolateOnIntegrationPoints<type, CohesiveReduceFunctionOpening>( \
displacement, opening, spatial_dimension, ghost_type, \
element_filter(type, ghost_type));
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING);
#undef COMPUTE_OPENING
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::updateEnergies(ElementType type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
if (Mesh::getKind(type) != _ek_cohesive)
return;
Vector<Real> b(spatial_dimension);
Vector<Real> h(spatial_dimension);
auto erev = reversible_energy(type, ghost_type).begin();
auto etot = total_energy(type, ghost_type).begin();
auto traction_it = tractions(type, ghost_type).begin(spatial_dimension);
auto traction_old_it =
tractions.previous(type, ghost_type).begin(spatial_dimension);
auto opening_it = opening(type, ghost_type).begin(spatial_dimension);
auto opening_old_it =
opening.previous(type, ghost_type).begin(spatial_dimension);
auto traction_end = tractions(type, ghost_type).end(spatial_dimension);
/// loop on each quadrature point
for (; traction_it != traction_end; ++traction_it, ++traction_old_it,
++opening_it, ++opening_old_it, ++erev,
++etot) {
/// trapezoidal integration
b = *opening_it;
b -= *opening_old_it;
h = *traction_old_it;
h += *traction_it;
*etot += .5 * b.dot(h);
*erev = .5 * traction_it->dot(*opening_it);
}
/// update old values
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getReversibleEnergy() {
AKANTU_DEBUG_IN();
Real erev = 0.;
/// integrate reversible energy for each type of elements
for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost,
_ek_cohesive)) {
erev +=
fem_cohesive.integrate(reversible_energy(type, _not_ghost), type,
_not_ghost, element_filter(type, _not_ghost));
}
AKANTU_DEBUG_OUT();
return erev;
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getDissipatedEnergy() {
AKANTU_DEBUG_IN();
Real edis = 0.;
/// integrate dissipated energy for each type of elements
for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost,
_ek_cohesive)) {
Array<Real> dissipated_energy(total_energy(type, _not_ghost));
dissipated_energy -= reversible_energy(type, _not_ghost);
edis += fem_cohesive.integrate(dissipated_energy, type, _not_ghost,
element_filter(type, _not_ghost));
}
AKANTU_DEBUG_OUT();
return edis;
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getContactEnergy() {
AKANTU_DEBUG_IN();
Real econ = 0.;
/// integrate contact energy for each type of elements
for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost,
_ek_cohesive)) {
auto & el_filter = element_filter(type, _not_ghost);
UInt nb_quad_per_el = fem_cohesive.getNbIntegrationPoints(type, _not_ghost);
UInt nb_quad_points = el_filter.size() * nb_quad_per_el;
Array<Real> contact_energy(nb_quad_points);
auto contact_traction_it =
contact_tractions(type, _not_ghost).begin(spatial_dimension);
auto contact_opening_it =
contact_opening(type, _not_ghost).begin(spatial_dimension);
/// loop on each quadrature point
for (UInt q = 0; q < nb_quad_points;
++contact_traction_it, ++contact_opening_it, ++q) {
contact_energy(q) = .5 * contact_traction_it->dot(*contact_opening_it);
}
econ += fem_cohesive.integrate(contact_energy, type, _not_ghost, el_filter);
}
AKANTU_DEBUG_OUT();
return econ;
}
/* -------------------------------------------------------------------------- */
Real MaterialCohesive::getEnergy(const std::string & type) {
AKANTU_DEBUG_IN();
if (type == "reversible")
return getReversibleEnergy();
else if (type == "dissipated")
return getDissipatedEnergy();
else if (type == "cohesive contact")
return getContactEnergy();
AKANTU_DEBUG_OUT();
return 0.;
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/model/time_step_solvers/time_step_solver_default.cc b/src/model/time_step_solvers/time_step_solver_default.cc
index cec055b5a..b115e433b 100644
--- a/src/model/time_step_solvers/time_step_solver_default.cc
+++ b/src/model/time_step_solvers/time_step_solver_default.cc
@@ -1,310 +1,302 @@
/**
* @file time_step_solver_default.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Sep 15 2015
* @date last modification: Wed Feb 21 2018
*
* @brief Default implementation of the time step solver
*
* @section LICENSE
*
* Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "time_step_solver_default.hh"
#include "dof_manager_default.hh"
#include "integration_scheme_1st_order.hh"
#include "integration_scheme_2nd_order.hh"
#include "mesh.hh"
#include "non_linear_solver.hh"
#include "pseudo_time.hh"
#include "sparse_matrix_aij.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
TimeStepSolverDefault::TimeStepSolverDefault(
DOFManagerDefault & dof_manager, const TimeStepSolverType & type,
NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id)
: TimeStepSolver(dof_manager, type, non_linear_solver, id, memory_id),
dof_manager(dof_manager), is_mass_lumped(false) {
switch (type) {
case _tsst_dynamic:
break;
case _tsst_dynamic_lumped:
this->is_mass_lumped = true;
break;
case _tsst_static:
/// initialize a static time solver for callback dofs
break;
default:
AKANTU_TO_IMPLEMENT();
}
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::setIntegrationScheme(
const ID & dof_id, const IntegrationSchemeType & type,
IntegrationScheme::SolutionType solution_type) {
if (this->integration_schemes.find(dof_id) !=
this->integration_schemes.end()) {
AKANTU_EXCEPTION("Their DOFs "
<< dof_id
<< " have already an integration scheme associated");
}
std::unique_ptr<IntegrationScheme> integration_scheme;
if (this->is_mass_lumped) {
switch (type) {
case _ist_forward_euler: {
integration_scheme = std::make_unique<ForwardEuler>(dof_manager, dof_id);
break;
}
case _ist_central_difference: {
integration_scheme =
std::make_unique<CentralDifference>(dof_manager, dof_id);
break;
}
default:
AKANTU_EXCEPTION(
"This integration scheme cannot be used in lumped dynamic");
}
} else {
switch (type) {
case _ist_pseudo_time: {
integration_scheme = std::make_unique<PseudoTime>(dof_manager, dof_id);
break;
}
case _ist_forward_euler: {
integration_scheme = std::make_unique<ForwardEuler>(dof_manager, dof_id);
break;
}
case _ist_trapezoidal_rule_1: {
integration_scheme =
std::make_unique<TrapezoidalRule1>(dof_manager, dof_id);
break;
}
case _ist_backward_euler: {
integration_scheme = std::make_unique<BackwardEuler>(dof_manager, dof_id);
break;
}
case _ist_central_difference: {
integration_scheme =
std::make_unique<CentralDifference>(dof_manager, dof_id);
break;
}
case _ist_fox_goodwin: {
integration_scheme = std::make_unique<FoxGoodwin>(dof_manager, dof_id);
break;
}
case _ist_trapezoidal_rule_2: {
integration_scheme =
std::make_unique<TrapezoidalRule2>(dof_manager, dof_id);
break;
}
case _ist_linear_acceleration: {
integration_scheme =
std::make_unique<LinearAceleration>(dof_manager, dof_id);
break;
}
case _ist_generalized_trapezoidal: {
integration_scheme =
std::make_unique<GeneralizedTrapezoidal>(dof_manager, dof_id);
break;
}
case _ist_newmark_beta:
integration_scheme = std::make_unique<NewmarkBeta>(dof_manager, dof_id);
break;
}
}
AKANTU_DEBUG_ASSERT(integration_scheme,
"No integration scheme was found for the provided types");
auto && matrices_names = integration_scheme->getNeededMatrixList();
for (auto && name : matrices_names) {
needed_matrices.insert({name, _mt_not_defined});
}
this->integration_schemes[dof_id] = std::move(integration_scheme);
this->solution_types[dof_id] = solution_type;
this->integration_schemes_owner.insert(dof_id);
}
/* -------------------------------------------------------------------------- */
bool TimeStepSolverDefault::hasIntegrationScheme(const ID & dof_id) const {
return this->integration_schemes.find(dof_id) !=
this->integration_schemes.end();
}
/* -------------------------------------------------------------------------- */
TimeStepSolverDefault::~TimeStepSolverDefault() = default;
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::solveStep(SolverCallback & solver_callback) {
this->solver_callback = &solver_callback;
this->solver_callback->beforeSolveStep();
this->non_linear_solver.solve(*this);
this->solver_callback->afterSolveStep();
this->solver_callback = nullptr;
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::predictor() {
- AKANTU_DEBUG_IN();
-
TimeStepSolver::predictor();
- for (auto & pair : this->integration_schemes) {
- auto & dof_id = pair.first;
+ for (auto && pair : this->integration_schemes) {
+ const auto & dof_id = pair.first;
auto & integration_scheme = pair.second;
if (this->dof_manager.hasPreviousDOFs(dof_id)) {
this->dof_manager.savePreviousDOFs(dof_id);
}
/// integrator predictor
integration_scheme->predictor(this->time_step);
}
-
- AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::corrector() {
AKANTU_DEBUG_IN();
TimeStepSolver::corrector();
for (auto & pair : this->integration_schemes) {
auto & dof_id = pair.first;
auto & integration_scheme = pair.second;
const auto & solution_type = this->solution_types[dof_id];
integration_scheme->corrector(solution_type, this->time_step);
/// computing the increment of dof if needed
if (this->dof_manager.hasDOFsIncrement(dof_id)) {
if (!this->dof_manager.hasPreviousDOFs(dof_id)) {
AKANTU_DEBUG_WARNING("In order to compute the increment of "
<< dof_id << " a 'previous' has to be registered");
continue;
}
Array<Real> & increment = this->dof_manager.getDOFsIncrement(dof_id);
Array<Real> & previous = this->dof_manager.getPreviousDOFs(dof_id);
UInt dof_array_comp = this->dof_manager.getDOFs(dof_id).getNbComponent();
auto prev_dof_it = previous.begin(dof_array_comp);
auto incr_it = increment.begin(dof_array_comp);
auto incr_end = increment.end(dof_array_comp);
increment.copy(this->dof_manager.getDOFs(dof_id));
for (; incr_it != incr_end; ++incr_it, ++prev_dof_it) {
*incr_it -= *prev_dof_it;
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::assembleMatrix(const ID & matrix_id) {
AKANTU_DEBUG_IN();
TimeStepSolver::assembleMatrix(matrix_id);
if (matrix_id != "J")
return;
for (auto & pair : this->integration_schemes) {
auto & dof_id = pair.first;
auto & integration_scheme = pair.second;
const auto & solution_type = this->solution_types[dof_id];
integration_scheme->assembleJacobian(solution_type, this->time_step);
}
this->dof_manager.applyBoundary("J");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::assembleResidual() {
- AKANTU_DEBUG_IN();
-
if (this->needed_matrices.find("M") != needed_matrices.end()) {
if (this->is_mass_lumped) {
this->assembleLumpedMatrix("M");
} else {
this->assembleMatrix("M");
}
}
TimeStepSolver::assembleResidual();
- for (auto & pair : this->integration_schemes) {
+ for (auto && pair : this->integration_schemes) {
auto & integration_scheme = pair.second;
integration_scheme->assembleResidual(this->is_mass_lumped);
}
-
- AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void TimeStepSolverDefault::assembleResidual(const ID & residual_part) {
AKANTU_DEBUG_IN();
if (this->needed_matrices.find("M") != needed_matrices.end()) {
if (this->is_mass_lumped) {
this->assembleLumpedMatrix("M");
} else {
this->assembleMatrix("M");
}
}
if (residual_part != "inertial") {
TimeStepSolver::assembleResidual(residual_part);
}
if (residual_part == "inertial") {
for (auto & pair : this->integration_schemes) {
auto & integration_scheme = pair.second;
integration_scheme->assembleResidual(this->is_mass_lumped);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/test/test_common/test_array.cc b/test/test_common/test_array.cc
index 6b587b845..081874ad9 100644
--- a/test/test_common/test_array.cc
+++ b/test/test_common/test_array.cc
@@ -1,285 +1,289 @@
/**
* @file test_array.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Nov 09 2017
* @date last modification: Fri Jan 26 2018
*
* @brief Test the arry class
*
* @section LICENSE
*
* Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "aka_types.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
#include <memory>
#include <typeindex>
#include <typeinfo>
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
class NonTrivial {
public:
NonTrivial() = default;
NonTrivial(int a) : a(a){};
bool operator==(const NonTrivial & rhs) { return a == rhs.a; }
int a{0};
};
bool operator==(const int & a, const NonTrivial & rhs) { return a == rhs.a; }
std::ostream & operator<<(std::ostream & stream, const NonTrivial & _this) {
stream << _this.a;
return stream;
}
/* -------------------------------------------------------------------------- */
using TestTypes = ::testing::Types<Real, UInt, NonTrivial>;
/* -------------------------------------------------------------------------- */
::testing::AssertionResult AssertType(const char * /*a_expr*/,
const char * /*b_expr*/,
const std::type_info & a,
const std::type_info & b) {
if (std::type_index(a) == std::type_index(b))
return ::testing::AssertionSuccess();
return ::testing::AssertionFailure()
<< debug::demangle(a.name()) << " != " << debug::demangle(b.name())
<< ") are different";
}
/* -------------------------------------------------------------------------- */
template <typename T> class ArrayConstructor : public ::testing::Test {
protected:
using type = T;
void SetUp() override { type_str = debug::demangle(typeid(T).name()); }
template <typename... P> decltype(auto) construct(P &&... params) {
return std::make_unique<Array<T>>(std::forward<P>(params)...);
}
protected:
std::string type_str;
};
TYPED_TEST_CASE(ArrayConstructor, TestTypes);
TYPED_TEST(ArrayConstructor, ConstructDefault1) {
auto array = this->construct();
EXPECT_EQ(0, array->size());
EXPECT_EQ(1, array->getNbComponent());
EXPECT_STREQ("", array->getID().c_str());
}
TYPED_TEST(ArrayConstructor, ConstructDefault2) {
auto array = this->construct(1000);
EXPECT_EQ(1000, array->size());
EXPECT_EQ(1, array->getNbComponent());
EXPECT_STREQ("", array->getID().c_str());
}
TYPED_TEST(ArrayConstructor, ConstructDefault3) {
auto array = this->construct(1000, 10);
EXPECT_EQ(1000, array->size());
EXPECT_EQ(10, array->getNbComponent());
EXPECT_STREQ("", array->getID().c_str());
}
TYPED_TEST(ArrayConstructor, ConstructDefault4) {
auto array = this->construct(1000, 10, "test");
EXPECT_EQ(1000, array->size());
EXPECT_EQ(10, array->getNbComponent());
EXPECT_STREQ("test", array->getID().c_str());
}
TYPED_TEST(ArrayConstructor, ConstructDefault5) {
auto array = this->construct(1000, 10, 1);
EXPECT_EQ(1000, array->size());
EXPECT_EQ(10, array->getNbComponent());
EXPECT_EQ(1, array->operator()(10, 6));
EXPECT_STREQ("", array->getID().c_str());
}
-TYPED_TEST(ArrayConstructor, ConstructDefault6) {
- typename TestFixture::type defaultv[2] = {0, 1};
+// TYPED_TEST(ArrayConstructor, ConstructDefault6) {
+// typename TestFixture::type defaultv[2] = {0, 1};
- auto array = this->construct(1000, 2, defaultv);
- EXPECT_EQ(1000, array->size());
- EXPECT_EQ(2, array->getNbComponent());
- EXPECT_EQ(1, array->operator()(10, 1));
- EXPECT_EQ(0, array->operator()(603, 0));
- EXPECT_STREQ("", array->getID().c_str());
-}
+// auto array = this->construct(1000, 2, defaultv);
+// EXPECT_EQ(1000, array->size());
+// EXPECT_EQ(2, array->getNbComponent());
+// EXPECT_EQ(1, array->operator()(10, 1));
+// EXPECT_EQ(0, array->operator()(603, 0));
+// EXPECT_STREQ("", array->getID().c_str());
+// }
/* -------------------------------------------------------------------------- */
template <typename T> class ArrayFixture : public ArrayConstructor<T> {
public:
void SetUp() override {
ArrayConstructor<T>::SetUp();
array = this->construct(1000, 10);
}
void TearDown() override { array.reset(nullptr); }
protected:
std::unique_ptr<Array<T>> array;
};
TYPED_TEST_CASE(ArrayFixture, TestTypes);
TYPED_TEST(ArrayFixture, Copy) {
Array<typename TestFixture::type> copy(*this->array);
EXPECT_EQ(1000, copy.size());
EXPECT_EQ(10, copy.getNbComponent());
EXPECT_NE(this->array->storage(), copy.storage());
}
TYPED_TEST(ArrayFixture, Set) {
auto & arr = *(this->array);
arr.set(12);
EXPECT_EQ(12, arr(156, 5));
EXPECT_EQ(12, arr(520, 7));
EXPECT_EQ(12, arr(999, 9));
}
TYPED_TEST(ArrayFixture, Resize) {
auto & arr = *(this->array);
- auto ptr = arr.storage();
+ auto * ptr = arr.storage();
arr.resize(0);
EXPECT_EQ(0, arr.size());
- EXPECT_EQ(ptr, arr.storage());
- EXPECT_EQ(1000, arr.getAllocatedSize());
+ EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
+ EXPECT_LE(0, arr.getAllocatedSize());
arr.resize(3000);
EXPECT_EQ(3000, arr.size());
- EXPECT_EQ(3000, arr.getAllocatedSize());
+ EXPECT_LE(3000, arr.getAllocatedSize());
+
+ ptr = arr.storage();
arr.resize(0);
EXPECT_EQ(0, arr.size());
- EXPECT_EQ(nullptr, arr.storage());
- EXPECT_EQ(0, arr.getAllocatedSize());
+ EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
+ EXPECT_LE(0, arr.getAllocatedSize());
}
TYPED_TEST(ArrayFixture, PushBack) {
auto & arr = *(this->array);
- auto ptr = arr.storage();
+ auto * ptr = arr.storage();
arr.resize(0);
EXPECT_EQ(0, arr.size());
- EXPECT_EQ(ptr, arr.storage());
- EXPECT_EQ(1000, arr.getAllocatedSize());
+ EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
+ EXPECT_LE(0, arr.getAllocatedSize());
arr.resize(3000);
EXPECT_EQ(3000, arr.size());
- EXPECT_EQ(3000, arr.getAllocatedSize());
+ EXPECT_LE(3000, arr.getAllocatedSize());
+
+ ptr = arr.storage();
arr.resize(0);
EXPECT_EQ(0, arr.size());
- EXPECT_EQ(nullptr, arr.storage());
- EXPECT_EQ(0, arr.getAllocatedSize());
+ EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
+ EXPECT_LE(0, arr.getAllocatedSize());
}
TYPED_TEST(ArrayFixture, ViewVector) {
auto && view = make_view(*this->array, 10);
EXPECT_NO_THROW(view.begin());
{
auto it = view.begin();
EXPECT_EQ(10, it->size());
EXPECT_PRED_FORMAT2(AssertType, typeid(*it),
typeid(Vector<typename TestFixture::type>));
EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]),
typeid(VectorProxy<typename TestFixture::type>));
}
}
TYPED_TEST(ArrayFixture, ViewMatrix) {
{
auto && view = make_view(*this->array, 2, 5);
EXPECT_NO_THROW(view.begin());
{
auto it = view.begin();
EXPECT_EQ(10, it->size());
EXPECT_EQ(2, it->size(0));
EXPECT_EQ(5, it->size(1));
EXPECT_PRED_FORMAT2(AssertType, typeid(*it),
typeid(Matrix<typename TestFixture::type>));
EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]),
typeid(MatrixProxy<typename TestFixture::type>));
}
}
}
TYPED_TEST(ArrayFixture, ViewVectorWrong) {
auto && view = make_view(*this->array, 11);
EXPECT_THROW(view.begin(), debug::ArrayException);
}
TYPED_TEST(ArrayFixture, ViewMatrixWrong) {
auto && view = make_view(*this->array, 3, 7);
EXPECT_THROW(view.begin(), debug::ArrayException);
}
TYPED_TEST(ArrayFixture, ViewMatrixIter) {
std::size_t count = 0;
for (auto && mat : make_view(*this->array, 10, 10)) {
EXPECT_EQ(100, mat.size());
EXPECT_EQ(10, mat.size(0));
EXPECT_EQ(10, mat.size(1));
EXPECT_PRED_FORMAT2(AssertType, typeid(mat),
typeid(Matrix<typename TestFixture::type>));
++count;
}
EXPECT_EQ(100, count);
}
TYPED_TEST(ArrayFixture, ConstViewVector) {
const auto & carray = *this->array;
auto && view = make_view(carray, 10);
EXPECT_NO_THROW(view.begin());
{
auto it = view.begin();
EXPECT_EQ(10, it->size());
EXPECT_PRED_FORMAT2(AssertType, typeid(*it),
typeid(Vector<typename TestFixture::type>));
EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]),
typeid(VectorProxy<typename TestFixture::type>));
}
}
} // namespace
diff --git a/test/test_model/test_model_solver/test_model_solver.cc b/test/test_model/test_model_solver/test_model_solver.cc
index bbca4d32c..9ccc544b0 100644
--- a/test/test_model/test_model_solver/test_model_solver.cc
+++ b/test/test_model/test_model_solver/test_model_solver.cc
@@ -1,148 +1,170 @@
/**
* @file test_model_solver.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Apr 13 2016
* @date last modification: Thu Feb 01 2018
*
* @brief Test default dof manager
*
* @section LICENSE
*
* Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_random_generator.hh"
#include "dof_manager.hh"
#include "dof_synchronizer.hh"
#include "mesh.hh"
#include "mesh_accessor.hh"
#include "model_solver.hh"
#include "non_linear_solver.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include "test_model_solver_my_model.hh"
/* -------------------------------------------------------------------------- */
#include <cmath>
/* -------------------------------------------------------------------------- */
using namespace akantu;
static void genMesh(Mesh & mesh, UInt nb_nodes);
static void printResults(MyModel & model, UInt nb_nodes);
Real F = -10;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize(argc, argv);
UInt prank = Communicator::getStaticCommunicator().whoAmI();
std::cout << std::setprecision(7);
UInt global_nb_nodes = 100;
Mesh mesh(1);
RandomGenerator<UInt>::seed(1);
if (prank == 0) {
genMesh(mesh, global_nb_nodes);
}
// std::cout << prank << RandGenerator<Real>::seed() << std::endl;
mesh.distribute();
MyModel model(F, mesh, false);
model.getNewSolver("static", _tsst_static, _nls_newton_raphson);
model.setIntegrationScheme("static", "disp", _ist_pseudo_time);
NonLinearSolver & solver = model.getDOFManager().getNonLinearSolver("static");
solver.set("max_iterations", 2);
model.solveStep();
printResults(model, global_nb_nodes);
finalize();
return EXIT_SUCCESS;
}
/* -------------------------------------------------------------------------- */
void genMesh(Mesh & mesh, UInt nb_nodes) {
MeshAccessor mesh_accessor(mesh);
Array<Real> & nodes = mesh_accessor.getNodes();
Array<UInt> & conn = mesh_accessor.getConnectivity(_segment_2);
nodes.resize(nb_nodes);
mesh_accessor.setNbGlobalNodes(nb_nodes);
for (UInt n = 0; n < nb_nodes; ++n) {
nodes(n, _x) = n * (1. / (nb_nodes - 1));
}
conn.resize(nb_nodes - 1);
for (UInt n = 0; n < nb_nodes - 1; ++n) {
conn(n, 0) = n;
conn(n, 1) = n + 1;
}
+ mesh_accessor.makeReady();
}
/* -------------------------------------------------------------------------- */
void printResults(MyModel & model, UInt nb_nodes) {
- UInt prank = model.mesh.getCommunicator().whoAmI();
- auto & sync = dynamic_cast<DOFManagerDefault &>(model.getDOFManager())
- .getSynchronizer();
-
- if (prank == 0) {
- Array<Real> global_displacement(nb_nodes);
- Array<Real> global_forces(nb_nodes);
- Array<bool> global_blocked(nb_nodes);
-
- sync.gather(model.forces, global_forces);
- sync.gather(model.displacement, global_displacement);
- sync.gather(model.blocked, global_blocked);
-
- auto force_it = global_forces.begin();
- auto disp_it = global_displacement.begin();
- auto blocked_it = global_blocked.begin();
+ if (model.mesh.isDistributed()) {
+ UInt prank = model.mesh.getCommunicator().whoAmI();
+ auto & sync = dynamic_cast<DOFManagerDefault &>(model.getDOFManager())
+ .getSynchronizer();
+
+ if (prank == 0) {
+ Array<Real> global_displacement(nb_nodes);
+ Array<Real> global_forces(nb_nodes);
+ Array<bool> global_blocked(nb_nodes);
+
+ sync.gather(model.forces, global_forces);
+ sync.gather(model.displacement, global_displacement);
+ sync.gather(model.blocked, global_blocked);
+
+ auto force_it = global_forces.begin();
+ auto disp_it = global_displacement.begin();
+ auto blocked_it = global_blocked.begin();
+
+ std::cout << "node"
+ << ", " << std::setw(8) << "disp"
+ << ", " << std::setw(8) << "force"
+ << ", " << std::setw(8) << "blocked" << std::endl;
+
+ UInt node = 0;
+ for (; disp_it != global_displacement.end();
+ ++disp_it, ++force_it, ++blocked_it, ++node) {
+ std::cout << node << ", " << std::setw(8) << *disp_it << ", "
+ << std::setw(8) << *force_it << ", " << std::setw(8)
+ << *blocked_it << std::endl;
+
+ std::cout << std::flush;
+ }
+ } else {
+ sync.gather(model.forces);
+ sync.gather(model.displacement);
+ sync.gather(model.blocked);
+ }
+ } else {
+ auto force_it = model.forces.begin();
+ auto disp_it = model.displacement.begin();
+ auto blocked_it = model.blocked.begin();
std::cout << "node"
- << ", " << std::setw(8) << "disp"
- << ", " << std::setw(8) << "force"
- << ", " << std::setw(8) << "blocked" << std::endl;
+ << ", " << std::setw(8) << "disp"
+ << ", " << std::setw(8) << "force"
+ << ", " << std::setw(8) << "blocked" << std::endl;
UInt node = 0;
- for (; disp_it != global_displacement.end();
- ++disp_it, ++force_it, ++blocked_it, ++node) {
- std::cout << node << ", " << std::setw(8) << *disp_it << ", "
- << std::setw(8) << *force_it << ", " << std::setw(8)
- << *blocked_it << std::endl;
-
- std::cout << std::flush;
- }
- } else {
- sync.gather(model.forces);
- sync.gather(model.displacement);
- sync.gather(model.blocked);
+ for (; disp_it != model.displacement.end();
+ ++disp_it, ++force_it, ++blocked_it, ++node) {
+ std::cout << node << ", " << std::setw(8) << *disp_it << ", "
+ << std::setw(8) << *force_it << ", " << std::setw(8)
+ << *blocked_it << std::endl;
+
+ std::cout << std::flush;
+ }
}
}