diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh index c0cda98d0..88e02987b 100644 --- a/src/common/aka_array_tmpl.hh +++ b/src/common/aka_array_tmpl.hh @@ -1,1182 +1,1182 @@ /** * @file aka_array_tmpl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Jul 15 2010 * @date last modification: Fri Feb 26 2021 * * @brief Inline functions of the classes Array and ArrayBase * * * @section LICENSE * * Copyright (©) 2010-2021 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 . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" // NOLINT /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ //#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__ //#define __AKANTU_AKA_ARRAY_TMPL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { namespace debug { struct ArrayException : public Exception {}; } // namespace debug /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(Int size, Int nb_component, const ID & id) : ArrayBase(id) { allocate(size, nb_component); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(Int size, Int nb_component, const_reference value, const ID & id) : ArrayBase(id) { allocate(size, nb_component, value); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::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(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer( const std::vector & vect) { this->data_storage = vect; this->size_ = vect.size(); this->nb_component = 1; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer & ArrayDataLayer::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 void ArrayDataLayer::allocate(Int new_size, Int nb_component) { this->nb_component = nb_component; this->resize(new_size); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(Int new_size, Int nb_component, const T & val) { this->nb_component = nb_component; this->resize(new_size, val); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(Int new_size) { this->data_storage.resize(new_size * this->nb_component); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(Int 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 void ArrayDataLayer::reserve(Int size, Int new_size) { if (new_size != -1) { this->data_storage.resize(new_size * this->nb_component); } 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 inline void ArrayDataLayer::push_back(const T & value) { this->data_storage.push_back(value); this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ /** * append a matrix or a vector to the array * @param new_elem a reference to a Matrix or Vector */ template template inline void ArrayDataLayer::push_back( const Eigen::MatrixBase & 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 << ")."); for (Idx i = 0; i < new_elem.size(); ++i) { this->data_storage.push_back(new_elem.array()[i]); } this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ template inline Int ArrayDataLayer::getAllocatedSize() const { return this->data_storage.capacity() / this->nb_component; } /* -------------------------------------------------------------------------- */ template inline Int ArrayDataLayer::getMemorySize() const { return this->data_storage.capacity() * sizeof(T); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class ArrayDataLayer : public ArrayBase { public: using value_type = T; using reference = value_type &; using pointer_type = value_type *; using size_type = typename ArrayBase::size_type; using const_reference = const value_type &; public: ~ArrayDataLayer() override { deallocate(); } /// Allocation of a new vector ArrayDataLayer(Int size = 0, Int nb_component = 1, const ID & id = "") : ArrayBase(id) { allocate(size, nb_component); } /// Allocation of a new vector with a default value ArrayDataLayer(Int size, Int 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.data(), this->size_ * this->nb_component, values); } /// Copy constructor (deep copy) explicit ArrayDataLayer(const std::vector & vect) { allocate(vect.size(), 1); std::copy_n(vect.data(), this->size_ * this->nb_component, values); } // copy operator inline ArrayDataLayer & operator=(const ArrayDataLayer & other) { if (this != &other) { allocate(other.size(), other.getNbComponent()); std::copy_n(other.data(), this->size_ * this->nb_component, values); } return *this; } // move constructor inline ArrayDataLayer(ArrayDataLayer && other) noexcept = default; // move assign inline ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default; protected: // deallocate the memory virtual void deallocate() { // NOLINTNEXTLINE(cppcoreguidelines-owning-memory, // cppcoreguidelines-no-malloc) free(this->values); } // allocate the memory virtual inline void allocate(Int size, Int nb_component) { if (size != 0) { // malloc can return a non NULL pointer in case size is 0 this->values = static_cast( // NOLINT std::malloc(nb_component * size * sizeof(T))); // NOLINT } 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 virtual inline void allocate(Int size, Int 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); } /// append a Vector or a Matrix template inline void push_back(const Eigen::MatrixBase & 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); make_view(*this, nb_component).begin()[this->size_ - 1].array() = new_elem.array(); } /// changes the allocated size but not the size virtual void reserve(Int size, Int new_size = Int(-1)) { auto tmp_size = this->size_; if (new_size != Int(-1)) { tmp_size = new_size; } this->resize(size); this->size_ = std::min(this->size_, tmp_size); } /// change the size of the Array virtual void resize(Int size) { if (size * this->nb_component == 0) { free(values); // NOLINT: cppcoreguidelines-no-malloc values = nullptr; this->allocated_size = 0; } else { if (this->values == nullptr) { this->allocate(size, this->nb_component); return; } Int diff = size - allocated_size; Int size_to_allocate = (std::abs(diff) > AKANTU_MIN_ALLOCATION) ? size : (diff > 0) ? allocated_size + AKANTU_MIN_ALLOCATION : allocated_size; if (size_to_allocate == allocated_size) { // otherwhy the reserve + push_back might fail... this->size_ = size; return; } auto * tmp_ptr = reinterpret_cast( // NOLINT 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 virtual void resize(Int size, const T & val) { Int tmp_size = this->size_; this->resize(size); if (size > tmp_size) { // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) std::fill_n(values + this->nb_component * tmp_size, (size - tmp_size) * this->nb_component, val); } } /// get the amount of space allocated in bytes inline size_type getMemorySize() const final { return this->allocated_size * this->nb_component * sizeof(T); } /// Get the real size allocated in memory inline Int getAllocatedSize() const { return this->allocated_size; } /// give the address of the memory allocated for this vector [[deprecated("use data instead to be stl compatible")]] T * storage() const { return values; }; T * data() const { return values; }; protected: /// allocation type agnostic data access T * values{nullptr}; Int allocated_size{0}; }; /* -------------------------------------------------------------------------- */ template inline auto Array::operator()(Int i, Int 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]; } /* -------------------------------------------------------------------------- */ template inline auto Array::operator()(Int i, Int 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 << "\""); // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) return this->values[i * this->nb_component + j]; } template inline auto Array::operator[](Int 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]; } /* -------------------------------------------------------------------------- */ template inline auto Array::operator[](Int 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 inline void Array::erase(Idx i) { AKANTU_DEBUG_IN(); 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 (Idx j = 0; j < this->nb_component; ++j) { this->values[i * this->nb_component + j] = this->values[(this->size_ - 1) * this->nb_component + j]; } } this->resize(this->size_ - 1); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ template template inline auto Array::erase(const view_iterator & it) { auto && curr = it.data(); Idx pos = (curr - this->values) / this->nb_component; erase(pos); view_iterator rit = it; return --rit; } /* -------------------------------------------------------------------------- */ /** * 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 Array & Array::operator-=(const Array & vect) { AKANTU_DEBUG_ASSERT((this->size_ == vect.size_) && (this->nb_component == vect.nb_component), "The too array don't have the same sizes"); T * a = this->values; T * b = vect.data(); for (Idx 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 Array & Array::operator+=(const Array & vect) { AKANTU_DEBUG_ASSERT((this->size_ == vect.size()) && (this->nb_component == vect.nb_component), "The too array don't have the same sizes"); T * a = this->values; T * b = vect.data(); for (Idx 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 */ template auto Array::operator*=(const T & alpha) -> Array & { T * a = this->values; for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a++ *= alpha; } return *this; } /* ------------------------------------------------------------------------- */ /** * Compare this array element by element to another. * @param other array to compare to * @return true it all element are equal and arrays have * the same shape, else false */ template bool Array::operator==(const Array & other) const { bool equal = this->nb_component == other.nb_component && this->size_ == other.size_ && this->id == other.id; if (not equal) { return false; } if (this->values == other.data()) { return true; } return std::equal(this->values, this->values + this->size_ * this->nb_component, other.data()); } /* ------------------------------------------------------------------------ */ template bool Array::operator!=(const Array & other) const { return !operator==(other); } /* ------------------------------------------------------------------------ */ /** * set all tuples of the array to a given vector or matrix * @param vm Matrix or Vector to fill the array with */ template template ::value> *> inline void Array::set(const C & elem) { AKANTU_DEBUG_ASSERT( this->nb_component == elem.array().size(), "The size of the object does not match the number of components"); for (auto && v : make_view(*this, this->nb_component)) { - v = elem.array(); + v = elem.reshaped(this->nb_component, 1); } } /* ------------------------------------------------------------------------ */ template void Array::append(const Array & other) { AKANTU_DEBUG_ASSERT( this->nb_component == other.nb_component, "Cannot append an array with a different number of component"); Idx old_size = this->size_; this->resize(this->size_ + other.size()); T * tmp = this->values + this->nb_component * old_size; std::copy_n(other.data(), other.size() * this->nb_component, tmp); } /* ------------------------------------------------------------------------ */ /* Functions Array */ /* ------------------------------------------------------------------------ */ template Array::Array(Int size, Int nb_component, const ID & id) : parent(size, nb_component, id) {} template <> inline Array::Array(Int size, Int nb_component, const ID & id) : parent(size, nb_component, "", id) {} /* ------------------------------------------------------------------------ */ template Array::Array(Int size, Int nb_component, const_reference value, const ID & id) : parent(size, nb_component, value, id) {} /* ------------------------------------------------------------------------ */ template Array::Array(const Array & vect, const ID & id) : parent(vect, id) {} /* ------------------------------------------------------------------------ */ template auto Array::operator=(const Array & other) -> Array & { AKANTU_DEBUG_WARNING("You are copying the array " << this->id << " are you sure it is on purpose"); if (&other == this) { return *this; } parent::operator=(other); return *this; } /* ------------------------------------------------------------------------ */ template Array::Array(const std::vector & vect) : parent(vect) {} /* ------------------------------------------------------------------------ */ template Array::~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 Idx Array::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 : Idx(-1); } /* ------------------------------------------------------------------------ */ template template ::value> *> inline Idx Array::find(const V & elem) { AKANTU_DEBUG_ASSERT(elem.size() == this->nb_component, "Cannot find an element with a wrong size (" << elem.size() << ") != " << this->nb_component); auto && view = make_view(*this, elem.size()); auto begin = view.begin(); auto end = view.end(); auto it = std::find(begin, end, elem); return (it != end) ? it - begin : Idx(-1); } /* ------------------------------------------------------------------------ */ /** * 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. * @param no_sanity_check turns off all checkes */ template void Array::copy(const Array & other, bool no_sanity_check) { AKANTU_DEBUG_IN(); if (not no_sanity_check and (other.nb_component != this->nb_component)) { AKANTU_ERROR("The two arrays do not have the same " "number of components"); } this->resize((other.size_ * other.nb_component) / this->nb_component); std::copy_n(other.data(), this->size_ * this->nb_component, this->values); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ template class ArrayPrintHelper { public: template static void print_content(const Array & vect, std::ostream & stream, int indent) { std::string space(indent, AKANTU_INDENT); stream << space << " + values : {"; for (Idx i = 0; i < vect.size(); ++i) { stream << "{"; for (Idx 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 { public: template static void print_content(__attribute__((unused)) const Array & vect, __attribute__((unused)) std::ostream & stream, __attribute__((unused)) int indent) {} }; /* ------------------------------------------------------------------------ */ template void Array::printself(std::ostream & stream, int indent) const { std::string space(indent, 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->getAllocatedSize() << std::endl; stream << space << " + memory size : " << printMemorySize(this->getMemorySize()) << std::endl; if (not AKANTU_DEBUG_LEVEL_IS_TEST()) { stream << space << " + address : " << std::hex << this->values << std::dec << std::endl; } stream.precision(prec); stream.flags(ff); if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) { ArrayPrintHelper::value>::print_content( *this, stream, indent); } stream << space << "]" << std::endl; } /* ------------------------------------------------------------------------ */ /* ArrayFilter */ /* ------------------------------------------------------------------------ */ template class ArrayFilter { public: /// const iterator template class const_iterator { public: Idx getCurrentIndex() { return (*this->filter_it * this->nb_item_per_elem + this->sub_element_counter); } inline const_iterator() = default; inline const_iterator(original_iterator origin_it, filter_iterator filter_it, Int nb_item_per_elem) : origin_it(std::move(origin_it)), filter_it(std::move(filter_it)), nb_item_per_elem(nb_item_per_elem), sub_element_counter(0){}; inline bool operator!=(const_iterator & other) const { return !((*this) == other); } inline bool operator==(const_iterator & other) const { return (this->origin_it == other.origin_it && this->filter_it == other.filter_it && this->sub_element_counter == other.sub_element_counter); } inline bool operator!=(const const_iterator & other) const { return !((*this) == other); } inline bool operator==(const const_iterator & other) const { return (this->origin_it == other.origin_it && this->filter_it == other.filter_it && this->sub_element_counter == other.sub_element_counter); } inline const_iterator & operator++() { ++sub_element_counter; if (sub_element_counter == nb_item_per_elem) { sub_element_counter = 0; ++filter_it; } return *this; }; inline decltype(auto) operator*() { return origin_it[nb_item_per_elem * (*filter_it) + sub_element_counter]; }; private: original_iterator origin_it; filter_iterator filter_it; /// the number of item per element Int nb_item_per_elem; /// counter for every sub element group Int sub_element_counter; }; using vector_iterator = void; // iterator>; using array_type = Array; using const_vector_iterator = const_iterator::const_scalar_iterator>; using value_type = typename array_type::value_type; private: /* ---------------------------------------------------------------------- */ /* Constructors/Destructors */ /* ---------------------------------------------------------------------- */ public: ArrayFilter(const Array & array, const Array & filter, Int nb_item_per_elem) : array(array), filter(filter), nb_item_per_elem(nb_item_per_elem){}; decltype(auto) begin_reinterpret(Int n, Int new_size) const { Int new_nb_item_per_elem = this->nb_item_per_elem; if (new_size != 0 && n != 0) new_nb_item_per_elem = this->array.getNbComponent() * this->filter.size() * this->nb_item_per_elem / (n * new_size); return const_vector_iterator(make_view(this->array, n).begin(), this->filter.begin(), new_nb_item_per_elem); }; decltype(auto) end_reinterpret(Int n, Int new_size) const { Int new_nb_item_per_elem = this->nb_item_per_elem; if (new_size != 0 && n != 0) new_nb_item_per_elem = this->array.getNbComponent() * this->filter.size() * this->nb_item_per_elem / (n * new_size); return const_vector_iterator(make_view(this->array, n).begin(), this->filter.end(), new_nb_item_per_elem); }; // vector_iterator begin_reinterpret(Int, Int) { throw; }; // vector_iterator end_reinterpret(Int, Int) { throw; }; /// return the size of the filtered array which is the filter size Int size() const { return this->filter.size() * this->nb_item_per_elem; }; /// the number of components of the filtered array Int getNbComponent() const { return this->array.getNbComponent(); }; /// tells if the container is empty bool empty() const [[gnu::warn_unused_result]] { return (size() == 0); } /* ---------------------------------------------------------------------- */ /* Class Members */ /* ---------------------------------------------------------------------- */ private: /// reference to array of data const Array & array; /// reference to the filter used to select elements const Array & filter; /// the number of item per element Int nb_item_per_elem; }; /* ------------------------------------------------------------------------ */ /* Begin/End functions implementation */ /* ------------------------------------------------------------------------ */ namespace detail { template constexpr auto take_front_impl(Tuple && t, std::index_sequence /*idxs*/) { return std::make_tuple(std::get(std::forward(t))...); } template constexpr auto take_front(Tuple && t) { return take_front_impl(std::forward(t), std::make_index_sequence{}); } template std::string to_string_all(T &&... t) { if (sizeof...(T) == 0) { return ""; } std::stringstream ss; bool noComma = true; ss << "("; (void)std::initializer_list{ (ss << (noComma ? "" : ", ") << t, noComma = false)...}; ss << ")"; return ss.str(); } template struct InstantiationHelper { template static auto instantiate(T && data, Ns... ns) { return std::make_unique(data, ns...); } }; template <> struct InstantiationHelper<0> { template static auto instantiate(T && data) { return data; } }; template decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) { const auto is_const_arr = std::is_const>::value; using type = ViewIteratorHelper_t; using iterator = std::conditional_t, view_iterator>; static_assert(sizeof...(Ns), "You should provide a least one size"); if (array.getNbComponent() * array.size() != Int(product_all(std::forward(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...)); } return aka::apply([&](auto... n) { return iterator(data, n...); }, take_front(std::make_tuple(ns...))); } } // namespace detail /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ template template inline auto Array::begin(Ns &&... ns) { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::end(Ns &&... ns) { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::begin(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::end(Ns &&... ns) const { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::cbegin(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::cend(Ns &&... ns) const { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::begin_reinterpret(Ns &&... ns) { return detail::get_iterator(*this, this->values, std::forward(ns)...); } template template inline auto Array::end_reinterpret(Ns &&... ns) { return detail::get_iterator( *this, this->values + detail::product_all(std::forward(ns)...), std::forward(ns)...); } template template inline auto Array::begin_reinterpret(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)...); } template template inline auto Array::end_reinterpret(Ns &&... ns) const { return detail::get_iterator( *this, this->values + detail::product_all(std::forward(ns)...), std::forward(ns)...); } /* ------------------------------------------------------------------------ */ /* Views */ /* ------------------------------------------------------------------------ */ namespace detail { template class ArrayView { using tuple = std::tuple; public: using size_type = Idx; ~ArrayView() = default; constexpr ArrayView(Array && array, Ns... ns) noexcept : array(array), sizes(std::move(ns)...) {} constexpr ArrayView(const ArrayView & array_view) = default; constexpr ArrayView(ArrayView && array_view) noexcept = default; constexpr ArrayView & operator=(const ArrayView & array_view) = default; constexpr ArrayView & operator=(ArrayView && array_view) noexcept = default; auto begin() { return aka::apply( [&](auto &&... ns) { return detail::get_iterator(array.get(), array.get().data(), std::forward(ns)...); }, sizes); } auto begin() const { return aka::apply( [&](auto &&... ns) { return detail::get_iterator(array.get(), array.get().data(), std::forward(ns)...); }, sizes); } auto end() { return aka::apply( [&](auto &&... ns) { return detail::get_iterator( array.get(), array.get().data() + detail::product_all(std::forward(ns)...), std::forward(ns)...); }, sizes); } auto end() const { return aka::apply( [&](auto &&... ns) { return detail::get_iterator( array.get(), array.get().data() + detail::product_all(std::forward(ns)...), std::forward(ns)...); }, sizes); } auto cbegin() const { return this->begin(); } auto cend() const { return this->end(); } constexpr auto size() const { return std::get::value - 1>(sizes); } constexpr auto dims() const { return std::tuple_size::value - 1; } private: std::reference_wrapper> array; tuple sizes; }; /* ---------------------------------------------------------------------- */ template class ArrayView &, Ns...> { using tuple = std::tuple; public: constexpr ArrayView(const ArrayFilter & array, Ns... ns) : array(array), sizes(std::move(ns)...) {} constexpr ArrayView(const ArrayView & array_view) = default; constexpr ArrayView(ArrayView && array_view) = default; constexpr ArrayView & operator=(const ArrayView & array_view) = default; constexpr ArrayView & operator=(ArrayView && array_view) = default; auto begin() const { return aka::apply( [&](auto &&... ns) { return array.get().begin_reinterpret( std::forward(ns)...); }, sizes); } auto end() const { return aka::apply( [&](auto &&... ns) { return array.get().end_reinterpret( std::forward(ns)...); }, sizes); } auto cbegin() const { return this->begin(); } auto cend() const { return this->end(); } constexpr auto size() const { return std::get::value - 1>(sizes); } constexpr auto dims() const { return std::tuple_size::value - 1; } private: std::reference_wrapper> array; tuple sizes; }; } // namespace detail /* ------------------------------------------------------------------------ */ template decltype(auto) make_view(Array && array, const Ns... ns) { static_assert(aka::conjunction>...>::value, "Ns should be integral types"); AKANTU_DEBUG_ASSERT((detail::product_all(ns...) != 0), "You must specify non zero dimensions"); auto size = std::forward(array).size() * std::forward(array).getNbComponent() / detail::product_all(ns...); return detail::ArrayView..., std::common_type_t>( std::forward(array), std::move(ns)..., size); } } // namespace akantu //#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */ diff --git a/src/common/aka_types_eigen_matrix_base_plugin.hh b/src/common/aka_types_eigen_matrix_base_plugin.hh index 883186d6a..6b17a3673 100644 --- a/src/common/aka_types_eigen_matrix_base_plugin.hh +++ b/src/common/aka_types_eigen_matrix_base_plugin.hh @@ -1,169 +1,169 @@ using size_type = Index; EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void zero(); template * = nullptr> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto operator()(Index c) { +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index c) { auto & d = this->derived(); return d.col(c); } template * = nullptr> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto operator()(Index c) const { +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index c) const { const auto & d = this->derived(); return d.col(c); } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index c) { return Base::operator()(c); } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index c) const { return Base::operator()(c); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index i, Index j) { return Base::operator()(i, j); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index i, Index j) const { return Base::operator()(i, j); } template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin(); template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end(); template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin() const; template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end() const; template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and not ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin(); template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and not ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end(); template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin() const; template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end() const; // clang-format off [[deprecated("use data instead to be stl compatible")]] Scalar * storage() { return this->data(); } [[deprecated("use data instead to be stl compatible")]] const Scalar * storage() const { return this->data(); } // clang-format on EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size() const { return this->rows() * this->cols(); } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size(Index i) const { AKANTU_DEBUG_ASSERT(i < 2, "This tensor has only " << 2 << " dimensions, not " << (i + 1)); return (i == 0) ? this->rows() : this->cols(); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void set(const Scalar & t) { this->fill(t); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eye(const Scalar & t = Scalar()) { (*this) = t * Matrix::Identity(this->rows(), this->cols()); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void clear() { this->fill(Scalar()); }; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto distance(const MatrixBase & other) const { return (*this - other).norm(); } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar doubleDot(const MatrixBase & other) const { Scalar sum = 0; eigen_assert(rows() == cols() and rows() == other.rows() and cols() == other.cols()); this->cwiseProduct(other).sum(); } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eig(const MatrixBase & other) const; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eig(const MatrixBase & values, const MatrixBase & vectors, bool sort = std::is_floating_point::value) const; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eigh(const MatrixBase & other) const; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eigh(const MatrixBase & values, const MatrixBase & vectors, bool sort = true) const; template inline bool operator<=(const MatrixBase & v) const { return this->isMuchSmallerThan(v); } template inline bool operator>=(const MatrixBase & v) const { return v.isMuchSmallerThan(*this); } template inline bool operator<(const MatrixBase & v) const { return (*this <= v) and (*this != v); } template inline bool operator>(const MatrixBase & v) const { return (*this >= v) and (*this != v); } diff --git a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh index 33c0d2ac1..e62bf6482 100644 --- a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh +++ b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh @@ -1,118 +1,115 @@ /** * @file element_class_segment_2_inline_impl.hh * * @author Emil Gallyamov * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jul 16 2010 * @date last modification: Fri Dec 11 2020 * * @brief Specialization of the element_class class for the type _segment_2 * * * @section LICENSE * * Copyright (©) 2010-2021 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 . * */ /** * @verbatim q --x--------|--------x---> x -1 0 1 @endverbatim * * @f{eqnarray*}{ * w_1(x) &=& 1/2(1 - x) \\ * w_2(x) &=& 1/2(1 + x) * @f} * * @f{eqnarray*}{ * x_{q} &=& 0 * @f} */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_segment_2, _gt_segment_2, _itp_lagrange_segment_2, _ek_regular, 1, _git_segment, 1); /* -------------------------------------------------------------------------- */ template <> template ::value> *> inline void InterpolationElement<_itp_lagrange_segment_2>::computeShapes( const Eigen::MatrixBase & natural_coords, Eigen::MatrixBase & N) { /// natural coordinate Real c = natural_coords(0); /// shape functions N(0) = 0.5 * (1 - c); N(1) = 0.5 * (1 + c); } /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_segment_2>::computeDNDS( const Eigen::MatrixBase & /*natural_coords*/, Eigen::MatrixBase & dnds) { /// dN1/de dnds(0, 0) = -.5; /// dN2/de dnds(0, 1) = .5; } /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_segment_2>::computeD2NDS2( const vector_type & /*natural_coords*/, matrix_type & d2nds2) { d2nds2.zero(); } /* -------------------------------------------------------------------------- */ template <> template inline Real GeometricalElement<_gt_segment_2>::getInradius( const Eigen::MatrixBase & coord) { - auto a{coord(0)}; - auto b{coord(1)}; - - return (b - a).norm(); + return (coord.col(1) - coord.col(0)).norm(); } /* -------------------------------------------------------------------------- */ template <> template inline void GeometricalElement<_gt_segment_2>::getNormal( const Eigen::MatrixBase & coord, Eigen::MatrixBase & normal) { AKANTU_DEBUG_ASSERT(normal.size() == 2, "The normal is only uniquely defined in 2D"); Math::normal(coord.col(0) - coord.col(1), normal); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/fe_engine/fe_engine_template_tmpl_field.hh b/src/fe_engine/fe_engine_template_tmpl_field.hh index cee4a13a3..0b65b7f92 100644 --- a/src/fe_engine/fe_engine_template_tmpl_field.hh +++ b/src/fe_engine/fe_engine_template_tmpl_field.hh @@ -1,420 +1,420 @@ /** * @file fe_engine_template_tmpl_field.hh * * @author Nicolas Richart * * @date creation: Wed Aug 09 2017 * @date last modification: Sat Mar 13 2021 * * @brief implementation of the assemble field s functions * * * @section LICENSE * * Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ //#include "fe_engine_template.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_ #define AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* Matrix lumping functions */ /* -------------------------------------------------------------------------- */ namespace fe_engine { namespace details { namespace { template void fillField(const Functor & field_funct, Array & field, Int nb_element, Int nb_integration_points, ElementType type, GhostType ghost_type) { auto nb_degree_of_freedom = field.getNbComponent(); field.resize(nb_integration_points * nb_element); Matrix mat(nb_degree_of_freedom, nb_integration_points); Element el{type, 0, ghost_type}; for (auto && data : enumerate(make_view(field, nb_degree_of_freedom, nb_integration_points))) { el.element = std::get<0>(data); - mat.fill(0.); + mat.zero(); field_funct(mat, el); - mat = std::get<1>(data); + std::get<1>(data) = mat; } } } // namespace } // namespace details } // namespace fe_engine #define DISPATCH_HELPER(Func, ...) \ auto && call = [&](auto && enum_type) { \ constexpr ElementType type = std::decay_t::value; \ this->Func(__VA_ARGS__); \ }; \ tuple_dispatch>(call, type); /* -------------------------------------------------------------------------- */ template