diff --git a/src/common/T4_map_proxy.hh b/src/common/T4_map_proxy.hh index 596130d..8783f3f 100644 --- a/src/common/T4_map_proxy.hh +++ b/src/common/T4_map_proxy.hh @@ -1,152 +1,167 @@ /** * file T4_map_proxy.hh * * @author Till Junge * * @date 19 Nov 2017 * * @brief Map type to allow fourth-order tensor-like maps on 2D matrices * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ +#include + #include -#include "Eigen/src/Core/util/Constants.h" +#include #ifndef T4_MAP_PROXY_H #define T4_MAP_PROXY_H namespace muSpectre { /* ---------------------------------------------------------------------- */ /** Proxy class mapping a fourth-order tensor onto a 2D matrix (in order to avoid the use of Eigen::Tensor. This class is, however byte-compatible with Tensors (i.e., you can map this onto a tensor instead of a matrix) **/ - template > class T4Map: - public Eigen::MapBase> { + public Eigen::MapBase> + { public: typedef Eigen::MapBase Base; EIGEN_DENSE_PUBLIC_INTERFACE(T4Map); - using PlainObjectType = Eigen::Matrix; + using matrix_type = Eigen::Matrix; + using PlainObjectType = + std::conditional_t; + using ConstType = T4Map; using Base::colStride; using Base::rowStride; using Base::IsRowMajor; typedef typename Base::PointerType PointerType; typedef PointerType PointerArgType; EIGEN_DEVICE_FUNC inline PointerType cast_to_pointer_type(PointerArgType ptr) { return ptr; } EIGEN_DEVICE_FUNC inline Eigen::Index innerStride() const { return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1; } + template + inline T4Map & operator=(const Eigen::MatrixBase & other) { + this->map = other; + return *this; + } + EIGEN_DEVICE_FUNC inline Eigen::Index outerStride() const { return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer() : IsVectorAtCompileTime ? this->size() : int(Flags)&Eigen::RowMajorBit ? this->cols() : this->rows(); } /** Constructor in the fixed-size case. * * \param dataPtr pointer to the array to map * \param stride optional Stride object, passing the strides. */ EIGEN_DEVICE_FUNC explicit inline T4Map(PointerArgType dataPtr, const StrideType& stride = StrideType()) : Base(cast_to_pointer_type(dataPtr)), m_stride(stride), map(cast_to_pointer_type(dataPtr)) { PlainObjectType::Base::_check_template_params(); } EIGEN_INHERIT_ASSIGNMENT_OPERATORS(T4Map); /** My accessor to mimick tensorial access **/ inline const Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) const { constexpr auto myColStride{ (colStride() == 1) ? colStride(): colStride()/Dim}; constexpr auto myRowStride{ (rowStride() == 1) ? rowStride(): rowStride()/Dim}; return this->operator()(i * myRowStride + j * myColStride, k * myRowStride + l * myColStride); } inline Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) { const auto myColStride{ (colStride() == 1) ? colStride(): colStride()/Dim}; const auto myRowStride{ (rowStride() == 1) ? rowStride(): rowStride()/Dim}; return this->map(i * myRowStride + j * myColStride, k * myRowStride + l * myColStride); } protected: StrideType m_stride; Eigen::Map map; }; } // muSpectre namespace Eigen { //! forward declarations template struct traits; /* ---------------------------------------------------------------------- */ namespace internal { - template - struct traits > + struct traits > : public traits> { using PlainObjectType = Matrix; typedef traits TraitsBase; enum { InnerStrideAtCompileTime = StrideType::InnerStrideAtCompileTime == 0 ? int(PlainObjectType::InnerStrideAtCompileTime) : int(StrideType::InnerStrideAtCompileTime), OuterStrideAtCompileTime = StrideType::OuterStrideAtCompileTime == 0 ? int(PlainObjectType::OuterStrideAtCompileTime) : int(StrideType::OuterStrideAtCompileTime), Alignment = int(MapOptions)&int(AlignedMask), Flags0 = TraitsBase::Flags & (~NestByRefBit), Flags = is_lvalue::value ? int(Flags0) : (int(Flags0) & ~LvalueBit) }; private: enum { Options }; // Expressions don't have Options }; } // namespace internal } // namespace Eigen #endif /* T4_MAP_PROXY_H */ diff --git a/src/common/eigen_tools.hh b/src/common/eigen_tools.hh index 3ef5dbd..d9d4cee 100644 --- a/src/common/eigen_tools.hh +++ b/src/common/eigen_tools.hh @@ -1,94 +1,140 @@ /** * file eigen_tools.hh * * @author Till Junge * * @date 20 Sep 2017 * * @brief small tools to be used with Eigen * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef EIGEN_TOOLS_H #define EIGEN_TOOLS_H #include #include #include "common/common.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ //! Creates a Eigen::Sizes type for a Tensor defined by an order and dim namespace internal { template struct SizesByOrderHelper { using Sizes = typename SizesByOrderHelper::Sizes; }; template struct SizesByOrderHelper<0, dim, dims...> { using Sizes = Eigen::Sizes; }; } // internal template struct SizesByOrder { static_assert(order > 0, "works only for order greater than zero"); using Sizes = typename internal::SizesByOrderHelper::Sizes; }; /* ---------------------------------------------------------------------- */ //! Call a passed lambda with the unpacked sizes as arguments namespace internal { /* ---------------------------------------------------------------------- */ template struct CallSizesHelper { static decltype(auto) call(Fun_t && fun) { static_assert(order > 0, "can't handle empty sizes b)"); return CallSizesHelper::call (fun); } }; /* ---------------------------------------------------------------------- */ template struct CallSizesHelper<0, Fun_t, dim, args...> { static decltype(auto) call(Fun_t && fun) { return fun(args...); } }; } // internal template inline decltype(auto) call_sizes(Fun_t && fun) { static_assert(order > 1, "can't handle empty sizes"); return internal::CallSizesHelper:: call(std::forward(fun)); } + + /** + * Structure to determine whether an expression can be evaluated into a Matrix, Array, etc. and which helps determine compile-time size + */ + struct EigenCheck { + template + constexpr static bool isDense(const Eigen::DenseBase & /*mat*/) { + return true; + } + + template + constexpr static bool isMatrix(const Eigen::MatrixBase & /*mat*/) { + return true; + } + template + constexpr static bool isFixed(T && t) { + static_assert + (isDense(t), + "The type you check for fixed size is not a Eigen::Dense type"); + using bT = std::remove_reference_t; + return ((bT::RowsAtCompileTime != Eigen::Dynamic) && + + (bT::ColsAtCompileTime != Eigen::Dynamic)); + } + + template + constexpr static bool isSquare(T && t) { + static_assert + (isDense(t), + "The type you check for fixed size is not a Eigen::Dense type"); + using bT = std::remove_reference_t; + return (bT::RowsAtCompileTime == bT::ColsAtCompileTime); + } + + template + constexpr static int TensorDim(T && t) { + static_assert + (isMatrix(t), "The type of t is not understood as an Eigen::Matrix"); + static_assert(isFixed(t), "t's dimension is not known at compile time"); + static_assert(isSquare(t), "t's matrix isn't square"); + return std::remove_reference_t::RowsAtCompileTime; + } + }; + + + } // muSpectre #endif /* EIGEN_TOOLS_H */ diff --git a/src/common/field.hh b/src/common/field.hh index c443b8e..59e5be5 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,408 +1,412 @@ /** * file field.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief header-only implementation of a field for field collections * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_H #define FIELD_H #include #include #include #include #include #include #include #include #include namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ template class FieldBase { protected: //! constructor //! unique name (whithin Collection) //! number of components //! collection to which this field belongs (eg, material, system) FieldBase(std::string unique_name, size_t nb_components, FieldCollection & collection); public: + using collection_t = FieldCollection; //! Copy constructor FieldBase(const FieldBase &other) = delete; //! Move constructor FieldBase(FieldBase &&other) noexcept = delete; //! Destructor virtual ~FieldBase() noexcept = default; //! Copy assignment operator FieldBase& operator=(const FieldBase &other) = delete; //! Move assignment operator FieldBase& operator=(FieldBase &&other) noexcept = delete; /* ---------------------------------------------------------------------- */ //!Identifying accessors //! return field name inline const std::string & get_name() const; //! return field type //inline const Field_t & get_type() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; //! return my collection (for iterating) inline const size_t & get_nb_components() const; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const = 0; virtual size_t size() const = 0; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) virtual void set_zero() = 0; //! give access to collections friend FieldCollection; protected: /* ---------------------------------------------------------------------- */ //! allocate memory etc virtual void resize(size_t size) = 0; const std::string name; const size_t nb_components; const FieldCollection & collection; private: }; /* ---------------------------------------------------------------------- */ //! declaraton for friending template class FieldMap; /* ---------------------------------------------------------------------- */ template class TypedFieldBase: public FieldBase { friend class FieldMap; public: + constexpr static auto nb_components{NbComponents}; using Parent = FieldBase; + using Parent::collection_t; + using Scalar = T; using Base = Parent; //using storage_type = Eigen::Array; using StoredType = Eigen::Array; using StorageType = std::vector>; TypedFieldBase(std::string unique_name, FieldCollection& collection); virtual ~TypedFieldBase() = default; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) inline void set_zero() override final; inline void push_back(const StoredType & value); template< class... Args> inline void emplace_back(Args&&... args); size_t size() const override final; protected: inline T* get_ptr_to_entry(const size_t&& index); inline T& get_ref_to_entry(const size_t&& index); inline virtual void resize(size_t size) override final; StorageType array; }; } // internal /* ---------------------------------------------------------------------- */ template class TensorField: public internal::TypedFieldBase { public: using Parent = internal::TypedFieldBase; using Base = typename Parent::Base; using Field_p = typename FieldCollection::Field_p; using component_type = T; //! Copy constructor TensorField(const TensorField &other) = delete; //! Move constructor TensorField(TensorField &&other) noexcept = delete; //! Destructor virtual ~TensorField() noexcept = default; //! Copy assignment operator TensorField& operator=(const TensorField &other) = delete; //! Move assignment operator TensorField& operator=(TensorField &&other) noexcept = delete; //! accessors inline Dim_t get_order() const; inline Dim_t get_dim() const; //! factory function template friend typename FieldType::Base& make_field(std::string unique_name, CollectionType & collection, Args&&... args); protected: //! constructor protected! TensorField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ template class MatrixField: public internal::TypedFieldBase { public: using Parent = internal::TypedFieldBase; using Base = typename Parent::Base; using Field_p = std::unique_ptr>; using component_type = T; //! Copy constructor MatrixField(const MatrixField &other) = delete; //! Move constructor MatrixField(MatrixField &&other) noexcept = delete; //! Destructor virtual ~MatrixField() noexcept = default; //! Copy assignment operator MatrixField& operator=(const MatrixField &other) = delete; //! Move assignment operator MatrixField& operator=(MatrixField &&other) noexcept = delete; //! accessors inline Dim_t get_nb_row() const; inline Dim_t get_nb_col() const; //! factory function template friend typename FieldType::Base& make_field(std::string unique_name, CollectionType & collection, Args&&... args); protected: //! constructor protected! MatrixField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ //! convenience alias template using ScalarField = TensorField; /* ---------------------------------------------------------------------- */ // Implementations namespace internal { /* ---------------------------------------------------------------------- */ template FieldBase::FieldBase(std::string unique_name, size_t nb_components_, FieldCollection & collection_) :name(unique_name), nb_components(nb_components_), collection(collection_) {} /* ---------------------------------------------------------------------- */ template inline const std::string & FieldBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template inline const FieldCollection & FieldBase:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ template inline const size_t & FieldBase:: get_nb_components() const { return this->nb_components; } /* ---------------------------------------------------------------------- */ template TypedFieldBase:: TypedFieldBase(std::string unique_name, FieldCollection & collection) :FieldBase(unique_name, NbComponents, collection){ static_assert ((std::is_arithmetic::value || std::is_same::value), "Use TypedFieldBase for integer, real or complex scalars for T"); static_assert(NbComponents > 0, "Only fields with more than 0 components"); } /* ---------------------------------------------------------------------- */ //! return type_id of stored type template const std::type_info & TypedFieldBase:: get_stored_typeid() const { return typeid(T); } /* ---------------------------------------------------------------------- */ template void TypedFieldBase:: set_zero() { std::fill(this->array.begin(), this->array.end(), T{}); } /* ---------------------------------------------------------------------- */ template size_t TypedFieldBase:: size() const { return this->array.size(); } /* ---------------------------------------------------------------------- */ template T* TypedFieldBase:: get_ptr_to_entry(const size_t&& index) { return &this->array[std::move(index)](0, 0); } /* ---------------------------------------------------------------------- */ template T& TypedFieldBase:: get_ref_to_entry(const size_t && index) { return this->array[std::move(index)](0, 0); } /* ---------------------------------------------------------------------- */ template void TypedFieldBase:: resize(size_t size) { this->array.resize(size); } /* ---------------------------------------------------------------------- */ template void TypedFieldBase:: push_back(const StoredType & value) { this->array.push_back(value); } /* ---------------------------------------------------------------------- */ template template void TypedFieldBase:: emplace_back(Args&&... args) { this->array.emplace_back(std::move(args...)); } } // internal /* ---------------------------------------------------------------------- */ //! Factory function, guarantees that only fields get created that are //! properly registered and linked to a collection. template typename FieldType::Base & make_field(std::string unique_name, FieldCollection & collection, Args&&... args) { auto && ptr = std::unique_ptr{ new FieldType(unique_name, collection, args...)}; auto& retref = *ptr; collection.register_field(std::move(ptr)); return retref; } /* ---------------------------------------------------------------------- */ template TensorField:: TensorField(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_order() const { return order; } /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_dim() const { return dim; } /* ---------------------------------------------------------------------- */ template MatrixField:: MatrixField(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_col() const { return NbCol; } /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_row() const { return NbRow; } } // muSpectre #endif /* FIELD_H */ diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh index 391ac26..9c355a3 100644 --- a/src/common/field_map_base.hh +++ b/src/common/field_map_base.hh @@ -1,499 +1,529 @@ /** * file field_map.hh * * @author Till Junge * * @date 12 Sep 2017 * * @brief Defined a strongly defines proxy that iterates efficiently over a field * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_H #define FIELD_MAP_H #include #include #include #include #include #include "common/field.hh" #include "common/common.hh" namespace muSpectre { namespace internal { //----------------------------------------------------------------------------// //! little helper to automate creation of const maps without duplication template struct const_corrector { using type = typename T::reference; }; template struct const_corrector { using type = typename T::const_reference; }; template using const_corrector_t = typename const_corrector::type; //----------------------------------------------------------------------------// template class FieldMap { public: + constexpr static auto nb_components{NbComponents}; using TypedField = ::muSpectre::internal::TypedFieldBase ; using Field = typename TypedField::Parent; using size_type = std::size_t; //! Default constructor FieldMap() = delete; FieldMap(Field & field); + template + FieldMap(TypedFieldBase & field); + //! Copy constructor FieldMap(const FieldMap &other) = delete; //! Move constructor FieldMap(FieldMap &&other) noexcept = delete; //! Destructor virtual ~FieldMap() noexcept = default; //! Copy assignment operator FieldMap& operator=(const FieldMap &other) = delete; //! Move assignment operator FieldMap& operator=(FieldMap &&other) noexcept = delete; - //! give static access to nb_components - inline constexpr static size_t nb_components(); - //! give human-readable field map type virtual std::string info_string() const = 0; //! return field name inline const std::string & get_name() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; //! member access needs to be implemented by inheriting classes //inline value_type operator[](size_t index); //inline value_type operator[](Ccoord ccord); //! check compatibility (must be called by all inheriting classes at the //! end of their constructors inline void check_compatibility(); //! convenience call to collection's size method inline size_t size() const; + //! compile-time compatibility check + template + struct is_compatible; + template class iterator { public: using value_type = const_corrector_t; using pointer = typename FullyTypedFieldMap::pointer; using difference_type = std::ptrdiff_t; using iterator_category = std::random_access_iterator_tag; using Ccoord = typename FieldCollection::Ccoord; //! Default constructor iterator() = delete; //! constructor inline iterator(FullyTypedFieldMap & fieldmap, bool begin=true); //! constructor for random access inline iterator(FullyTypedFieldMap & fieldmap, size_t index); //! Copy constructor iterator(const iterator &other)= default; //! Move constructor iterator(iterator &&other) noexcept = default; //! Destructor virtual ~iterator() noexcept = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) noexcept = default; //! pre-increment inline iterator & operator++(); //! post-increment inline iterator operator++(int); //! dereference inline value_type operator*(); //! member of pointer inline pointer operator->(); //! pre-decrement inline iterator & operator--(); //! post-decrement inline iterator operator--(int); //! access subscripting inline value_type operator[](difference_type diff); //! equality inline bool operator==(const iterator & other) const; //! inequality inline bool operator!=(const iterator & other) const; //! div. comparisons inline bool operator<(const iterator & other) const; inline bool operator<=(const iterator & other) const; inline bool operator>(const iterator & other) const; inline bool operator>=(const iterator & other) const; //! additions, subtractions and corresponding assignments inline iterator operator+(difference_type diff) const; inline iterator operator-(difference_type diff) const; inline iterator& operator+=(difference_type diff); inline iterator& operator-=(difference_type diff); //! get pixel coordinates inline Ccoord get_ccoord() const; //! ostream operator (mainly for debug friend std::ostream & operator<<(std::ostream & os, const iterator& it) { if (isConst) { os << "const "; } os << "iterator on field '" << it.fieldmap.get_name() << "', entry " << it.index; return os; } protected: const FieldCollection & collection; FullyTypedFieldMap & fieldmap; size_t index; private: }; protected: inline T* get_ptr_to_entry(size_t index); inline T& get_ref_to_entry(size_t index); const FieldCollection & collection; TypedField & field; private: }; } // internal namespace internal { /* ---------------------------------------------------------------------- */ template FieldMap:: FieldMap(Field & field) :collection(field.get_collection()), field(static_cast(field)) { static_assert(NbComponents>0, "Only fields with more than 0 components allowed"); } + /* ---------------------------------------------------------------------- */ + template + template + FieldMap:: + FieldMap(TypedFieldBase & field) + :collection(field.get_collection()), field(static_cast(field)) { + static_assert(std::is_same::value, + "The field does not have the expected FieldCollection type"); + static_assert(std::is_same::value, + "The field does not have the expected Scalar type"); + static_assert((NbC == NbComponents), + "The field does not have the expected number of components"); + } + /* ---------------------------------------------------------------------- */ template void FieldMap:: check_compatibility() { if (typeid(T).hash_code() != this->field.get_stored_typeid().hash_code()) { throw FieldInterpretationError ("Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' of type '" + this->field.get_stored_typeid().name() + "'"); } //check size compatibility if (NbComponents != this->field.get_nb_components()) { throw FieldInterpretationError ("Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' with " + std::to_string(this->field.get_nb_components()) + " components"); } } /* ---------------------------------------------------------------------- */ template size_t FieldMap:: size() const { return this->collection.size(); } /* ---------------------------------------------------------------------- */ - template - constexpr size_t - FieldMap:: - nb_components() { - return NbComponents; - } + template + template + struct FieldMap::is_compatible { + constexpr static bool explain() { + static_assert + (std::is_same::value, + "The field does not have the expected FieldCollection type"); + static_assert + (std::is_same::value, + "The // FIXME: eld does not have the expected Scalar type"); + static_assert((Field::nb_components == NbComponents), + "The field does not have the expected number of components"); + //The static asserts wouldn't pass in the incompatible case, so this is it + return true; + } + constexpr static bool value{std::is_base_of::value}; + }; /* ---------------------------------------------------------------------- */ template const std::string & FieldMap:: get_name() const { return this->field.get_name(); } /* ---------------------------------------------------------------------- */ template const FieldCollection & FieldMap:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ // Iterator implementations //! constructor template template FieldMap::iterator:: iterator(FullyTypedFieldMap & fieldmap, bool begin) :collection(fieldmap.get_collection()), fieldmap(fieldmap), index(begin ? 0 : fieldmap.field.size()) {} /* ---------------------------------------------------------------------- */ //! constructor for random access template template FieldMap::iterator:: iterator(FullyTypedFieldMap & fieldmap, size_t index) :collection(fieldmap.collection), fieldmap(fieldmap), index(index) {} /* ---------------------------------------------------------------------- */ //! pre-increment template template FieldMap::iterator & FieldMap::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ //! post-increment template template FieldMap::iterator FieldMap::iterator:: operator++(int) { iterator current = *this; this->index++; return current; } /* ---------------------------------------------------------------------- */ //! dereference template template typename FieldMap::template iterator::value_type FieldMap::iterator:: operator*() { return this->fieldmap.template operator[](this->index); } /* ---------------------------------------------------------------------- */ //! member of pointer template template typename FullyTypedFieldMap::pointer FieldMap::iterator:: operator->() { return this->fieldmap.ptr_to_val_t(this->index); } /* ---------------------------------------------------------------------- */ //! pre-decrement template template FieldMap::iterator & FieldMap::iterator:: operator--() { this->index--; return *this; } /* ---------------------------------------------------------------------- */ //! post-decrement template template FieldMap::iterator FieldMap::iterator:: operator--(int) { iterator current = *this; this->index--; return current; } /* ---------------------------------------------------------------------- */ //! Access subscripting template template typename FieldMap::template iterator::value_type FieldMap::iterator:: operator[](difference_type diff) { return this->fieldmap[this->index+diff]; } /* ---------------------------------------------------------------------- */ //! equality template template bool FieldMap::iterator:: operator==(const iterator & other) const { return (this->index == other.index && &this->fieldmap == &other.fieldmap); } /* ---------------------------------------------------------------------- */ //! inquality template template bool FieldMap::iterator:: operator!=(const iterator & other) const { return !(*this == other); } /* ---------------------------------------------------------------------- */ //! div. comparisons template template bool FieldMap::iterator:: operator<(const iterator & other) const { return (this->index < other.index); } template template bool FieldMap::iterator:: operator<=(const iterator & other) const { return (this->index <= other.index); } template template bool FieldMap::iterator:: operator>(const iterator & other) const { return (this->index > other.index); } template template bool FieldMap::iterator:: operator>=(const iterator & other) const { return (this->index >= other.index); } /* ---------------------------------------------------------------------- */ //! additions, subtractions and corresponding assignments template template FieldMap::iterator FieldMap::iterator:: operator+(difference_type diff) const { return iterator(this->fieldmap, this->index + diff); } template template FieldMap::iterator FieldMap::iterator:: operator-(difference_type diff) const { return iterator(this->fieldmap, this->index - diff); } template template FieldMap::iterator & FieldMap::iterator:: operator+=(difference_type diff) { this->index += diff; return *this; } template template FieldMap::iterator & FieldMap::iterator:: operator-=(difference_type diff) { this->index -= diff; return *this; } /* ---------------------------------------------------------------------- */ //! get pixel coordinates template template typename FieldCollection::Ccoord FieldMap::iterator:: get_ccoord() const { return this->collection.get_ccoord(this->index); } ////----------------------------------------------------------------------------// //template //std::ostream & operator << //(std::ostream &os, // const typename FieldMap:: // template iterator & it) { // os << "iterator on field '" // << it.field.get_name() // << "', entry " << it.index; // return os; //} /* ---------------------------------------------------------------------- */ template T* FieldMap::get_ptr_to_entry(size_t index) { return this->field.get_ptr_to_entry(std::move(index)); } /* ---------------------------------------------------------------------- */ template T& FieldMap::get_ref_to_entry(size_t index) { return this->field.get_ref_to_entry(std::move(index)); } } // internal } // muSpectre #endif /* FIELD_MAP_H */ diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh index cad2979..055b53f 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,204 +1,227 @@ /** * file field_map_matrixlike.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Eigen-Matrix and -Array maps over strongly typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_MATRIXLIKE_H #define FIELD_MAP_MATRIXLIKE_H #include #include #include "common/field_map_base.hh" #include "common/T4_map_proxy.hh" namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ enum class Map_t{Matrix, Array, T4Matrix}; template struct NameWrapper { const static std::string field_info_root; }; template<> const std::string NameWrapper::field_info_root{"Array"}; template<> const std::string NameWrapper::field_info_root{"Matrix"}; template<> const std::string NameWrapper::field_info_root{"T4Matrix"}; /* ---------------------------------------------------------------------- */ template class MatrixLikeFieldMap: public FieldMap { public: using Parent = FieldMap; using Ccoord = Ccoord_t; using value_type = EigenArray; using reference = value_type; // since it's a resource handle using const_reference = const value_type; using size_type = typename Parent::size_type; using pointer = std::unique_ptr; using TypedField = typename Parent::TypedField; using Field = typename TypedField::Parent; using iterator = typename Parent::template iterator; using const_iterator= typename Parent::template iterator; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor MatrixLikeFieldMap() = delete; + /** + * Constructor using a (non-typed) field. Compatibility is enforced at + * runtime. This should not be a performance concern, as this constructor + * will not be called in anny inner loops (if used as intended). + */ MatrixLikeFieldMap(Field & field); + /** + * Constructor using a typed field. Compatibility is enforced + * statically. It is not always possible to call this constructor, as the + * type of the field might not be known at compile time. + */ + template + MatrixLikeFieldMap(TypedFieldBase & field); + //! Copy constructor MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = delete; //! Move constructor MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = delete; //! Destructor virtual ~MatrixLikeFieldMap() noexcept = default; //! Copy assignment operator MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete; //! Move assignment operator MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) noexcept = delete; //! give human-readable field map type inline std::string info_string() const override final; //! member access template inline ref operator[](size_type index); inline reference operator[](const Ccoord& ccoord); //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} inline const_iterator cbegin(){return const_iterator(*this);} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);}; inline const_iterator cend(){return const_iterator(*this, false);} protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); const static std::string field_info_root; private: }; template const std::string MatrixLikeFieldMap:: field_info_root{NameWrapper::field_info_root}; /* ---------------------------------------------------------------------- */ template MatrixLikeFieldMap:: MatrixLikeFieldMap(Field & field) :Parent(field) { this->check_compatibility(); } + + /* ---------------------------------------------------------------------- */ + template + template + MatrixLikeFieldMap:: + MatrixLikeFieldMap(TypedFieldBase & field) + :Parent(field) { + } + /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string MatrixLikeFieldMap:: info_string() const { std::stringstream info; info << field_info_root << "(" << typeid(typename EigenArray::value_type).name() << ", " << EigenArray::RowsAtCompileTime << "x" << EigenArray::ColsAtCompileTime << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template template ref MatrixLikeFieldMap:: operator[](size_type index) { return ref(this->get_ptr_to_entry(index)); } template typename MatrixLikeFieldMap::reference MatrixLikeFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); return reference(this->get_ptr_to_entry(std::move(index))); } /* ---------------------------------------------------------------------- */ template typename MatrixLikeFieldMap::pointer MatrixLikeFieldMap:: ptr_to_val_t(size_type index) { return std::make_unique (this->get_ptr_to_entry(std::move(index))); } } // internal /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using MatrixFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Matrix>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate - template + template using T4MatrixFieldMap = internal::MatrixLikeFieldMap , + T4Map, internal::Map_t::T4Matrix>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen array map as iterate template using ArrayFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Array>; } // muSpectre #endif /* FIELD_MAP_MATRIXLIKE_H */ diff --git a/src/common/field_map_tensor.hh b/src/common/field_map_tensor.hh index ee27a3e..6c4204c 100644 --- a/src/common/field_map_tensor.hh +++ b/src/common/field_map_tensor.hh @@ -1,169 +1,169 @@ /** * file field_map_tensor.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Defines an Eigen-Tensor map over strongly typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_TENSOR_H #define FIELD_MAP_TENSOR_H #include #include #include "common/eigen_tools.hh" #include "common/field_map_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template class TensorFieldMap: public internal::FieldMap::Sizes::total_size > { public: using Parent = internal::FieldMap; + TensorFieldMap::nb_components>; using Ccoord = Ccoord_t; using Sizes = typename SizesByOrder::Sizes; using T_t = Eigen::TensorFixedSize; using value_type = Eigen::TensorMap; using reference = value_type; // since it's a resource handle using const_reference = Eigen::TensorMap; using size_type = typename Parent::size_type; using pointer = std::unique_ptr; using TypedField = typename Parent::TypedField; using Field = typename TypedField::Parent; using iterator = typename Parent::template iterator; using const_iterator = typename Parent::template iterator; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor TensorFieldMap() = delete; TensorFieldMap(Field & field); //! Copy constructor TensorFieldMap(const TensorFieldMap &other) = delete; //! Move constructor TensorFieldMap(TensorFieldMap &&other) noexcept = delete; //! Destructor virtual ~TensorFieldMap() noexcept = default; //! Copy assignment operator TensorFieldMap& operator=(const TensorFieldMap &other) = delete; //! Move assignment operator TensorFieldMap& operator=(TensorFieldMap &&other) noexcept = delete; //! give human-readable field map type inline std::string info_string() const override final; //! member access template inline ref_t operator[](size_type index); inline reference operator[](const Ccoord & ccoord); //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} inline const_iterator cbegin(){return const_iterator(*this);} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);} inline const_iterator cend(){return const_iterator(*this, false);} protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); private: }; /* ---------------------------------------------------------------------- */ template TensorFieldMap:: TensorFieldMap(Field & field) :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string TensorFieldMap::info_string() const { std::stringstream info; info << "Tensor(" << typeid(T).name() << ", " << order << "_o, " << dim << "_d)"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template template ref_t TensorFieldMap::operator[](size_type index) { auto && lambda = [this, &index](auto&&...tens_sizes) { return ref_t(this->get_ptr_to_entry(index), tens_sizes...); }; return call_sizes(lambda); } template typename TensorFieldMap::reference TensorFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); auto && lambda = [this, &index](auto&&...sizes) { return ref_t(this->get_ptr_to_entry(index), sizes...); }; return call_sizes(lambda); } /* ---------------------------------------------------------------------- */ //! for sad, legacy iterator use. Don't use unless you have to. template typename TensorFieldMap::pointer TensorFieldMap:: ptr_to_val_t(size_type index) { auto && lambda = [this, &index](auto&&... tens_sizes) { return std::make_unique(this->get_ptr_to_entry(index), tens_sizes...); }; return call_sizes(lambda); } } // muSpectre #endif /* FIELD_MAP_TENSOR_H */ diff --git a/src/materials/CMakeLists.txt b/src/materials/CMakeLists.txt index f80ac63..e67ab10 100644 --- a/src/materials/CMakeLists.txt +++ b/src/materials/CMakeLists.txt @@ -1,9 +1,9 @@ set (materials_SRC material_base.cc - material_hyper_elastic1.cc + #material_hyper_elastic1.cc #material_linear_elastic.cc ) add_library(materials ${materials_SRC}) target_link_libraries(materials Eigen3::Eigen common) diff --git a/src/materials/material_hyper_elastic1.cc b/src/materials/material_hyper_elastic1.cc index a4efbd4..7f8ae79 100644 --- a/src/materials/material_hyper_elastic1.cc +++ b/src/materials/material_hyper_elastic1.cc @@ -1,68 +1,70 @@ /** * file material_hyper_elastic1.cc * * @author Till Junge * * @date 14 Nov 2017 * * @brief Implementation for materialhyperelastic1 * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_hyper_elastic1.hh" #include "common/tensor_algebra.hh" #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template MaterialHyperElastic1::MaterialHyperElastic1(std::string name, Real young, Real poisson) :Parent(name), young{young}, poisson{poisson}, lambda{young*poisson/((1+poisson)*(1-2*poisson))}, mu{young/(2*(1+poisson))}, C{lambda*Tensors::outer(Tensors::I2(),Tensors::I2()) + 2*mu*Tensors::I4S()} {} /* ---------------------------------------------------------------------- */ template template decltype(auto) MaterialHyperElastic1::evaluate_stress(s_t && E) { return E.trace()*lambda * Strain_t::Identity() + 2*mu*E; } /* ---------------------------------------------------------------------- */ template template decltype(auto) MaterialHyperElastic1::evaluate_stress_tangent(s_t && E) { - return std::forward_as_tuple(this->evaluate_stress(std::move(E)), this->C); + return std::forward_as_tuple(this->evaluate_stress(std::move(E)), + Tangent_t(const_cast(this->C.data()))); + throw 12; } template class MaterialHyperElastic1<2, 2>; template class MaterialHyperElastic1<2, 3>; template class MaterialHyperElastic1<3, 3>; } // muSpectre diff --git a/src/materials/material_hyper_elastic1.hh b/src/materials/material_hyper_elastic1.hh index d587aaf..1e3ef6b 100644 --- a/src/materials/material_hyper_elastic1.hh +++ b/src/materials/material_hyper_elastic1.hh @@ -1,101 +1,101 @@ /** * file material_hyper_elastic1.hh * * @author Till Junge * * @date 13 Nov 2017 * * @brief Implementation for hyperelastic reference material like in de Geus * 2017. This follows the simplest and likely not most efficient * implementation (with exception of the Python law) * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_muSpectre_base.hh" #ifndef MATERIAL_HYPER_ELASTIC1_H #define MATERIAL_HYPER_ELASTIC1_H namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class MaterialHyperElastic1: public MaterialMuSpectre, DimS, DimM> { public: using Parent = MaterialMuSpectre; using GFieldCollection_t = typename Parent::GFieldCollection_t; // declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; // declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; // declare whether the derivative of stress with respect to strain is uniform constexpr static bool uniform_stiffness = true; // declare the type in which you wish to receive your strain measure using StrainMap_t = MatrixFieldMap; using StressMap_t = StrainMap_t; - using TangentMap_t = T4MatrixFieldMap; + using TangentMap_t = T4MatrixFieldMap; using Strain_t = typename StrainMap_t::const_reference; using Stress_t = typename StressMap_t::reference; - using Tangent_t = typename TangentMap_t::reference; + using Tangent_t = typename TangentMap_t::reference::ConstType; using Stiffness_t = Eigen::TensorFixedSize >; //! Default constructor MaterialHyperElastic1() = delete; //! Copy constructor MaterialHyperElastic1(const MaterialHyperElastic1 &other) = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialHyperElastic1(std::string name, Real young, Real poisson); //! Move constructor MaterialHyperElastic1(MaterialHyperElastic1 &&other) noexcept = delete; //! Destructor virtual ~MaterialHyperElastic1() noexcept = default; //! Copy assignment operator MaterialHyperElastic1& operator=(const MaterialHyperElastic1 &other) = delete; //! Move assignment operator MaterialHyperElastic1& operator=(MaterialHyperElastic1 &&other) noexcept = delete; template decltype(auto) evaluate_stress(s_t && E); template decltype(auto) evaluate_stress_tangent(s_t && E); const Tangent_t & get_stiffness() const; protected: const Real young, poisson, lambda, mu; const Stiffness_t C; private: }; } // muSpectre #endif /* MATERIAL_HYPER_ELASTIC1_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index e082852..3f97428 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,351 +1,383 @@ /** * file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include "stdexcept" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #ifndef MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H namespace muSpectre { //! 'Material' is a CRTP template class MaterialMuSpectre: public MaterialBase { public: enum class NeedTangent {yes, no}; using Parent = MaterialBase; using MFieldCollection_t = typename Parent::MFieldCollection_t; using StressField_t = typename Parent::StressField_t; using StrainField_t = typename Parent::StrainField_t; using TangentField_t = typename Parent::TangentField_t; + + //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) noexcept = delete; //! Destructor virtual ~MaterialMuSpectre() noexcept = default; //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) noexcept = delete; //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form) override final; //! computes stress and tangent modulus virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) override final; protected: //! computes stress with the formulation available at compile time template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P); //! computes stress with the formulation available at compile time template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate over //! that does or does not include tangent moduli template class iterable_proxy; private: }; + /* ---------------------------------------------------------------------- */ + template + MaterialMuSpectre:: + MaterialMuSpectre(std::string name):Parent(name) { + using stress_compatible = typename Material::StressMap_t:: + template is_compatible; + using strain_compatible = typename Material::StrainMap_t:: + template is_compatible; + using tangent_compatible = typename Material::TangentMap_t:: + template is_compatible; + + static_assert((stress_compatible::value && + stress_compatible::explain()), + "The material's declared stress map is not compatible " + "with the stress field. More info in previously shown " + "assert."); + + static_assert((strain_compatible::value && + strain_compatible::explain()), + "The material's declared strain map is not compatible " + "with the strain field. More info in previously shown " + "assert."); + + static_assert((tangent_compatible::value && + tangent_compatible::explain()), + "The material's declared tangent map is not compatible " + "with the tangent field. More info in previously shown " + "assert."); + } + /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses(const StrainField_t &F, StressField_t &P, muSpectre::Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, muSpectre::Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P, K); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P, K); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K){ //! in small strain problems, we always want Cauchy stress as a function //! of the infinitesimal strain tensor using IterableSmallStrain_t = iterable_proxy; using IterableFiniteStrain_t = iterable_proxy; using iterable = std::conditional_t <(Form==Formulation::small_strain), IterableSmallStrain_t, IterableFiniteStrain_t>; constexpr auto stored_strain_m = get_stored_strain_type(Form); constexpr auto expected_strain_m = get_formulation_strain_type(Form, Material::strain_measure); /* This lambda is executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ auto constitutive_law = [this](auto && F, auto && ... internal_variables) { auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent moduli auto && stress_tgt = static_cast(*this).evaluate_stress_tangent(std::move(strain), internal_variables...); if (Form == Formulation::small_strain) { return std::move(stress_tgt); } else { auto && stress = std::get<0>(stress_tgt); auto && tangent = std::get<1>(stress_tgt); return MatTB::PK1_stress (F, stress, tangent); } }; iterable it{*this}; for (auto && tuples: it) { // the iterator yields a pair of tuples. this first tuple contains // references to stress and stiffness in the global arrays, the second // contains references to the deformation gradient and internal variables // (some of them are const). auto && stress_tgt = std::get<0>(tuples); auto && inputs = std::get<1>(tuples); stress_tgt = std::apply(constitutive_law, std::move(inputs)); } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template template::NeedTangent NeedTgt> class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterable_proxy(const MaterialMuSpectre & mat); using Strain_t = typename Material::Strain_t; using Stress_t = typename Material::Stress_t; using Tangent_t = typename Material::Tangent_t; //! Copy constructor iterable_proxy(const iterable_proxy &other) = default; //! Move constructor iterable_proxy(iterable_proxy &&other) noexcept = default; //! Destructor virtual ~iterable_proxy() noexcept = default; //! Copy assignment operator iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator iterable_proxy& operator=(iterable_proxy &&other) noexcept = default; class iterator { public: using value_type = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>, std::tuple, std::tuple>>; using iterator_category = std::forward_iterator_tag; //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterator(const MaterialMuSpectre & mat, bool begin = true); //! Copy constructor iterator(const iterator &other) = default; //! Move constructor iterator(iterator &&other) noexcept = default; //! Destructor virtual ~iterator() noexcept = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) noexcept = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const MaterialMuSpectre & material; size_t index; private: }; iterator begin() {return iterator(*this);} iterator end() {return iterator(*this, false);} protected: const MaterialMuSpectre & material; private: }; /* ---------------------------------------------------------------------- */ template template::NeedTangent Tangent> MaterialMuSpectre::iterable_proxy:: iterable_proxy(const MaterialMuSpectre & mat) : material{mat} {} /* ---------------------------------------------------------------------- */ template template::NeedTangent Tangent> MaterialMuSpectre::iterable_proxy:: iterator::iterator(const MaterialMuSpectre & mat, bool begin) : material{mat}, index{begin ? 0:mat.internal_fields.size()} {} /* ---------------------------------------------------------------------- */ //! pre-increment template template::NeedTangent Tangent> typename MaterialMuSpectre:: template iterable_proxy::iterator & MaterialMuSpectre:: template iterable_proxy::iterator:: operator++() { ++this->index; return *this; } /* ---------------------------------------------------------------------- */ //! dereference template template::NeedTangent Tangent> typename MaterialMuSpectre:: template iterable_proxy::iterator::value_type MaterialMuSpectre:: template iterable_proxy::iterator:: operator*() { static_assert(Tangent, "Template specialisation failed"); return std::make_tuple (this->material.internal_fields); } } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh index dce28e4..af55928 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,306 +1,328 @@ /** * file materials_toolbox.hh * * @author Till Junge * * @date 02 Nov 2017 * * @brief collection of common continuum mechanics tools * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MATERIALS_TOOLBOX_H #define MATERIALS_TOOLBOX_H +#include #include #include #include #include #include "common/common.hh" #include "common/tensor_algebra.hh" +#include "common/eigen_tools.hh" +#include "common/T4_map_proxy.hh" namespace muSpectre { namespace MatTB { class MaterialsToolboxError:public std::runtime_error{ public: explicit MaterialsToolboxError(const std::string& what) :std::runtime_error(what){} explicit MaterialsToolboxError(const char * what) :std::runtime_error(what){} }; /* ---------------------------------------------------------------------- */ /** Structure for functions returning one strain measure as a function of another **/ namespace internal { template struct ConvertStrain { template inline static decltype(auto) compute(Strain_t&& input) { // transparent case, in which no conversion is required: // just a perfect forwarding static_assert ((In == Out), "This particular strain conversion is not implemented"); return std::forward(input); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Green-Lagrange strain from the transformation gradient **/ template <> template decltype(auto) ConvertStrain:: compute(Strain_t && F) { return (F.transpose()*F - Strain_t::PlainObject::Identity()); } } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning one strain measure as a function of //! another template decltype(auto) convert_strain(Strain_t && strain) { return internal::ConvertStrain::compute(std::move(strain)); }; /* ---------------------------------------------------------------------- */ /** Structure for functions returning PK1 stress from other stress measures **/ namespace internal { - template + template struct PK1_stress { template inline static decltype(auto) compute(Strain_t && /*strain*/, Stress_t && /*stress*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } template inline static decltype(auto) compute(Strain_t && /*strain*/, Stress_t && /*stress*/, Tangent_t && /*stiffness*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress **/ - template - struct PK1_stress: - public PK1_stress + struct PK1_stress: + public PK1_stress { template inline static decltype(auto) compute(Strain_t && /*dummy*/, Stress_t && P) { return std::forward(P); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress *and* stiffness is given with respect to the transformation gradient **/ - template <> - struct PK1_stress: - public PK1_stress + struct PK1_stress: + public PK1_stress { + using Parent = PK1_stress; + using Parent::compute; + template decltype(auto) compute(Strain_t && /*dummy*/, Stress_t && P, Tangent_t && K) { return std::forward_as_tuple(P, K); } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) */ - template - struct PK1_stress: - public PK1_stress + struct PK1_stress: + public PK1_stress { template inline static decltype(auto) compute(Strain_t && F, Stress_t && S) { return F*S; } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) derived * with respect to Green-Lagrange strain */ - template <> - struct PK1_stress: - public PK1_stress + struct PK1_stress: + public PK1_stress { - + using Parent = PK1_stress; + using Parent::compute; + template inline static decltype(auto) compute(Strain_t && F, Stress_t && S, Tangent_t && C) { using T4 = typename Tangent_t::PlainObject; + using Tmap = T4Map; T4 K; + Tmap Kmap{K.data()}; K.setZero(); constexpr int dim{Strain_t::ColsAtCompileTime}; for (int i = 0; i < dim; ++i) { for (int m = 0; m < dim; ++m) { for (int n = 0; n < dim; ++n) { - K(i,m,i,n) += S(m,n); + Kmap(i,m,i,n) += S(m,n); for (int j = 0; j < dim; ++j) { for (int r = 0; r < dim; ++r) { for (int s = 0; s < dim; ++s) { - K(i,m,j,n) += F(i,r)*C(r,m,n,s)*(F(j,s)); + Kmap(i,m,j,n) += F(i,r)*C(r,m,n,s)*(F(j,s)); } } } } } } - return K; + auto && P = compute(std::forward(F), + std::forward(S)); + return std::forward_as_tuple(std::move(P), K); } }; } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) { - return internal::PK1_stress::compute(std::move(stress), - std::move(strain)); + constexpr Dim_t dim{EigenCheck::TensorDim(strain)}; + static_assert((dim == EigenCheck::TensorDim(stress)), + "Stress and strain tensors have differing dimensions"); + return internal::PK1_stress::compute + (std::move(stress), std::move(strain)); }; /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template - decltype(auto) PK1_stress(Strain_t && strain, - Stress_t && stress, - Tangent_t && stiffness) { - return internal::PK1_stress::compute(std::move(stress), - std::move(strain), - std::move(stiffness)); + decltype(auto) PK1_stress(Strain_t && strain, + Stress_t && stress, + Tangent_t && tangent) { + constexpr Dim_t dim{EigenCheck::TensorDim(strain)}; + static_assert((dim == EigenCheck::TensorDim(stress)), + "Stress and strain tensors have differing dimensions"); + static_assert((dim*dim == EigenCheck::TensorDim(tangent)), + "Stress and tangent tensors have differing dimensions"); + return internal::PK1_stress::compute + (std::move(stress), std::move(strain), std::move(tangent)); }; /* ---------------------------------------------------------------------- */ /** Structure for functions returning ∂strain/∂F (where F is the transformation * gradient. (e.g. for a material law expressed as PK2 stress S as * function of Green-Lagrange strain E, this would return F^T ⊗̲ I) (⊗_ * refers to Curnier's notation) WARNING: this function assumes that your * stress measure has the two minor symmetries due to stress equilibrium * and the major symmetry due to the existance of an elastic potential * W(E). Do not rely on this function if you are using exotic stress * measures **/ namespace internal { template struct StrainDerivative { template inline static decltype(auto) compute (Grad_t && grad) { // the following test always fails to generate a compile-time error static_assert((StrainM == StrainMeasure::Almansi) && (StrainM == StrainMeasure::Biot), "The requested StrainDerivative calculation is not " "implemented. You either made a programming mistake or " "need to implement it as a specialisation of this " "function. See " "StrainDerivative for an " "example."); } }; template <> template decltype(auto) StrainDerivative:: compute(Grad_t && grad) { } } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning ∂strain/∂F template auto StrainDerivative = [](auto && F) decltype(auto) { // the following test always fails to generate a compile-time error static_assert((StrainM == StrainMeasure::Almansi) && (StrainM == StrainMeasure::Biot), "The requested StrainDerivative calculation is not " "implemented. You either made a programming mistake or " "need to implement it as a specialisation of this " "function. See " "StrainDerivative for an " "example."); }; /* ---------------------------------------------------------------------- */ template<> auto StrainDerivative = [] (auto && F) -> decltype(auto) { constexpr size_t order{2}; constexpr Dim_t dim{F.Dimensions[0]}; return Tensors::outer_under (F.shuffle(std::array{1,0}), Tensors::I2()); }; } // MatTB } // muSpectre #endif /* MATERIALS_TOOLBOX_H */ diff --git a/tests/test_t4_map.cc b/tests/test_t4_map.cc index 6007092..8e7941f 100644 --- a/tests/test_t4_map.cc +++ b/tests/test_t4_map.cc @@ -1,90 +1,99 @@ /** * file test_t4_map.cc * * @author Till Junge * * @date 20 Nov 2017 * * @brief Test the fourth-order map on second-order tensor implementation * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include #include "common/common.hh" #include "tests.hh" #include "common/T4_map_proxy.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(T4map_tests); /** * Test fixture for construction of T4Map for the time being, symmetry is not * exploited */ template struct T4_fixture { T4_fixture():matrix{}, tensor(matrix.data()){} EIGEN_MAKE_ALIGNED_OPERATOR_NEW; using M4 = Eigen::Matrix; using T4 = T4Map; constexpr static Dim_t dim{Dim}; M4 matrix; T4 tensor; }; using fix_collection = boost::mpl::list, T4_fixture>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, fix_collection, F) { BOOST_CHECK_EQUAL(F::tensor.cols(), F::dim*F::dim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(write_access_test, F, fix_collection, F) { auto & t4 = F::tensor; constexpr Dim_t dim{F::dim}; Eigen::TensorFixedSize> t4c; Eigen::Map t4c_map(t4c.data()); for (Dim_t i = 0; i < F::dim; ++i) { for (Dim_t j = 0; j < F::dim; ++j) { for (Dim_t k = 0; k < F::dim; ++k) { for (Dim_t l = 0; l < F::dim; ++l) { t4(i,j,k,l) = 1000*(i+1) + 100*(j+1) + 10*(k+1) + l+1; t4c(i,j,k,l) = 1000*(i+1) + 100*(j+1) + 10*(k+1) + l+1; } } } } for (Dim_t i = 0; i < ipow(dim,4); ++i) { BOOST_CHECK_EQUAL(F::matrix.data()[i], t4c.data()[i]); } } + BOOST_FIXTURE_TEST_CASE_TEMPLATE(assign_matrix_test, F, fix_collection, F) { + decltype(F::matrix) matrix; + matrix.setRandom(); + F::tensor = matrix; + for (Dim_t i = 0; i < ipow(F::dim,4); ++i) { + BOOST_CHECK_EQUAL(F::matrix.data()[i], matrix.data()[i]); + } + } + BOOST_AUTO_TEST_SUITE_END(); } // muSpectre