diff --git a/src/common/field.hh b/src/common/field.hh index f6ed3c2..0dd38d6 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,413 +1,429 @@ /** * 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; 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 const T* get_ptr_to_entry(const size_t&& index) const; + inline const T& get_ref_to_entry(const size_t&& index) const; 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 + const T* TypedFieldBase:: + get_ptr_to_entry(const size_t&& index) const { + return &this->array[std::move(index)](0, 0); + } + + /* ---------------------------------------------------------------------- */ + template + const T& TypedFieldBase:: + get_ref_to_entry(const size_t && index) const { + 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_collection_global.hh b/src/common/field_collection_global.hh index 8289924..aa158d2 100644 --- a/src/common/field_collection_global.hh +++ b/src/common/field_collection_global.hh @@ -1,157 +1,165 @@ /** * file field_collection_global.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for global 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_COLLECTION_GLOBAL_H #define FIELD_COLLECTION_GLOBAL_H #include "common/field_collection_base.hh" #include "common/ccoord_operations.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ //! DimS spatial dimension (dimension of problem) //! DimM material_dimension (dimension of constitutive law) template class GlobalFieldCollection: public FieldCollectionBase> { public: using Parent = FieldCollectionBase >; using Ccoord = typename Parent::Ccoord; using Field_p = typename Parent::Field_p; //! Default constructor GlobalFieldCollection(); //! Copy constructor GlobalFieldCollection(const GlobalFieldCollection &other) = delete; //! Move constructor GlobalFieldCollection (GlobalFieldCollection &&other) noexcept = delete; //! Destructor virtual ~GlobalFieldCollection() noexcept = default; //! Copy assignment operator GlobalFieldCollection& operator=(const GlobalFieldCollection &other) = delete; //! Move assignment operator GlobalFieldCollection& operator=(GlobalFieldCollection &&other) noexcept = delete; /** allocate memory, etc. At this point, the collection is informed aboud the size and shape of the domain (through the sizes parameter). The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' any field of a different size is wrong TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(Ccoord sizes); //! return the pixel sizes inline const Ccoord & get_sizes() const; //! returns the linear index corresponding to cell coordinates - inline size_t get_index(Ccoord && ccoord) const; + template + inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; static constexpr inline Dim_t spatial_dim() {return DimS;} static constexpr inline Dim_t material_dim() {return DimM;} protected: //! number of discretisation cells in each of the DimS spatial directions Ccoord sizes{0}; private: }; /* ---------------------------------------------------------------------- */ template GlobalFieldCollection::GlobalFieldCollection() :Parent(){} /* ---------------------------------------------------------------------- */ template void GlobalFieldCollection:: initialise(Ccoord sizes) { this->size_ = std::accumulate(sizes.begin(), sizes.end(), 1, std::multiplies()); this->sizes = sizes; std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! return the pixel sizes template const typename GlobalFieldCollection::Ccoord & GlobalFieldCollection::get_sizes() const { return this->sizes; } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename GlobalFieldCollection::Ccoord GlobalFieldCollection::get_ccoord(size_t index) const { return CcoordOps::get_ccoord(this->get_sizes(), std::move(index)); } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template + template size_t - GlobalFieldCollection::get_index(Ccoord && ccoord) const { - return CcoordOps::get_index(this->get_sizes(), std::move(ccoord)); + GlobalFieldCollection::get_index(CcoordRef && ccoord) const { + static_assert(std::is_same< + Ccoord, + std::remove_const_t< + std::remove_reference_t>>::value, + "can only be called with values or references of Ccoord"); + return CcoordOps::get_index(this->get_sizes(), + std::forward(ccoord)); } } // muSpectre #endif /* FIELD_COLLECTION_GLOBAL_H */ diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh index 3edefaf..2c52a25 100644 --- a/src/common/field_collection_local.hh +++ b/src/common/field_collection_local.hh @@ -1,159 +1,166 @@ /** * file field_collection_local.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for local 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_COLLECTION_LOCAL_H #define FIELD_COLLECTION_LOCAL_H #include "common/field_collection_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ //! DimS spatial dimension (dimension of problem) //! DimM material_dimension (dimension of constitutive law) template class LocalFieldCollection: public FieldCollectionBase> { public: using Parent = FieldCollectionBase>; using Ccoord = typename Parent::Ccoord; using Field_p = typename Parent::Field_p; //! Default constructor LocalFieldCollection(); //! Copy constructor LocalFieldCollection(const LocalFieldCollection &other) = delete; //! Move constructor LocalFieldCollection(LocalFieldCollection &&other) noexcept = delete; //! Destructor virtual ~LocalFieldCollection() noexcept = default; //! Copy assignment operator LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete; //! Move assignment operator LocalFieldCollection& operator=(LocalFieldCollection &&other) noexcept = delete; //! add a pixel/voxel to the field collection inline void add_pixel(const Ccoord & local_ccoord); /** allocate memory, etc. at this point, the field collection knows how many entries it should have from the size of the coords containes (which grows by one every time add_pixel is called. The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' any field of a different size is wrong TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(); //! returns the linear index corresponding to cell coordinates - inline size_t get_index(Ccoord && ccoord) const; + template + inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; protected: //! container of pixel coords for non-global collections std::vector ccoords{}; //! container of indices for non-global collections (slow!) std::map indices{}; private: }; /* ---------------------------------------------------------------------- */ template LocalFieldCollection::LocalFieldCollection() :Parent(){} /* ---------------------------------------------------------------------- */ template void LocalFieldCollection:: add_pixel(const Ccoord & local_ccoord) { if (this->is_initialised) { throw FieldCollectionError ("once a field collection has been initialised, you can't add new " "pixels."); } this->indices[local_ccoord] = this->ccoords.size(); this->ccoords.push_back(local_ccoord); } /* ---------------------------------------------------------------------- */ template void LocalFieldCollection:: initialise() { this->size_ = this->ccoords.size(); std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template + template size_t - LocalFieldCollection::get_index(Ccoord && ccoord) const { - return this->indices[std::move(ccoord)]; + LocalFieldCollection::get_index(CcoordRef && ccoord) const { + static_assert(std::is_same< + Ccoord, + std::remove_const_t< + std::remove_reference_t>>::value, + "can only be called with values or references of Ccoord"); + return this->indices[std::forward(ccoord)]; } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename LocalFieldCollection::Ccoord LocalFieldCollection::get_ccoord(size_t index) const { return this->ccoords[std::move(index)]; } } // muSpectre #endif /* FIELD_COLLECTION_LOCAL_H */ diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh index 60197f3..79ac830 100644 --- a/src/common/field_map_base.hh +++ b/src/common/field_map_base.hh @@ -1,553 +1,575 @@ /** * 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 "field_collection_base.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_nc = TypedFieldBase ; using TypedField = std::conditional_t; using Field = typename TypedField::Parent; using size_type = std::size_t; - + using pointer = std::conditional_t; //! Default constructor FieldMap() = delete; template FieldMap(std::enable_if_t field); template FieldMap(std::enable_if_t field); template FieldMap(TypedFieldBase & field); //! Copy constructor - FieldMap(const FieldMap &other) = delete; + FieldMap(const FieldMap &other) = default; //! Move constructor - FieldMap(FieldMap &&other) noexcept = delete; + FieldMap(FieldMap &&other) noexcept = default; //! Destructor virtual ~FieldMap() noexcept = default; //! Copy assignment operator FieldMap& operator=(const FieldMap &other) = delete; //! Move assignment operator FieldMap& operator=(FieldMap &&other) noexcept = delete; //! 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 { static_assert(!((ConstIter==false) && (ConstField==true)), "You can't have a non-const iterator over a const " "field"); 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; using reference = typename FullyTypedFieldMap::reference; //! 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 (ConstIter) { 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 pointer get_ptr_to_entry(size_t index); + inline const T* get_ptr_to_entry(size_t index) const; inline T& get_ref_to_entry(size_t index); + inline const T& get_ref_to_entry(size_t index) const; const FieldCollection & collection; TypedField & field; private: }; } // internal namespace internal { /* ---------------------------------------------------------------------- */ template template FieldMap:: FieldMap(std::enable_if_t 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(std::enable_if_t 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()) { std::string err{"Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' of type '" + this->field.get_stored_typeid().name() + "'"}; throw FieldInterpretationError (err); } //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 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) { + typename FieldMap::pointer + FieldMap:: + get_ptr_to_entry(size_t index) { + return this->field.get_ptr_to_entry(std::move(index)); + } + + /* ---------------------------------------------------------------------- */ + template + const T* + FieldMap:: + get_ptr_to_entry(size_t index) const { return this->field.get_ptr_to_entry(std::move(index)); } /* ---------------------------------------------------------------------- */ template T& - FieldMap::get_ref_to_entry(size_t index) { + FieldMap:: + get_ref_to_entry(size_t index) { + return this->field.get_ref_to_entry(std::move(index)); + } + + /* ---------------------------------------------------------------------- */ + template + const T& + FieldMap:: + get_ref_to_entry(size_t index) const { 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 ceb0798..2401f91 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,274 +1,275 @@ /** * 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 { }; template<> struct NameWrapper { static std::string field_info_root() {return "Array";} }; template<> struct NameWrapper { static std::string field_info_root() {return "Matrix";} }; template<> struct NameWrapper { static std::string field_info_root() {return "T4Matrix";} }; /* ---------------------------------------------------------------------- */ template class MatrixLikeFieldMap: public FieldMap { public: using Parent = FieldMap; using Ccoord = Ccoord_t; using value_type = EigenArray; using const_reference = const value_type; using reference = std::conditional_t; // since it's a resource handle using size_type = typename Parent::size_type; using pointer = std::unique_ptr; using TypedField = typename Parent::TypedField; using Field = typename TypedField::Parent; using const_iterator= typename Parent::template iterator; using iterator = std::conditional_t< ConstField, 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). */ template MatrixLikeFieldMap(std::enable_if_t field); template MatrixLikeFieldMap(std::enable_if_t 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; + MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = default; //! Move constructor - MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = delete; + MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = default; //! 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 template MatrixLikeFieldMap:: MatrixLikeFieldMap(std::enable_if_t field) :Parent(field) { static_assert((isntConst != ConstField), "let the compiler deduce isntConst, this is a SFINAE " "parameter"); this->check_compatibility(); } /* ---------------------------------------------------------------------- */ template template MatrixLikeFieldMap:: MatrixLikeFieldMap(std::enable_if_t field) :Parent(field) { static_assert((isConst == ConstField), "let the compiler deduce isntConst, this is a SFINAE " "parameter"); 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); + size_t index{}; + 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 >, Eigen::Map>>, internal::Map_t::Matrix, ConstField>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using T4MatrixFieldMap = internal::MatrixLikeFieldMap , internal::Map_t::T4Matrix, MapConst>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen array map as iterate template using ArrayFieldMap = internal::MatrixLikeFieldMap >, Eigen::Map>>, internal::Map_t::Array, ConstField>; } // muSpectre #endif /* FIELD_MAP_MATRIXLIKE_H */ diff --git a/src/common/field_map_scalar.hh b/src/common/field_map_scalar.hh index 4ad0897..a9a52b1 100644 --- a/src/common/field_map_scalar.hh +++ b/src/common/field_map_scalar.hh @@ -1,165 +1,165 @@ /** * file field_map_scalar.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief maps over scalar 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_SCALAR_H #define FIELD_MAP_SCALAR_H #include "common/field_map_base.hh" namespace muSpectre { template class ScalarFieldMap : public internal::FieldMap { public: using parent = internal::FieldMap; using Ccoord = Ccoord_t; using value_type = T; using const_reference = const value_type &; using reference = std::conditional_t; using size_type = typename parent::size_type; using pointer = T*; using Field = typename parent::Field; using const_iterator= typename parent::template iterator; using iterator = std::conditional_t >; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor ScalarFieldMap() = delete; template ScalarFieldMap(std::enable_if_t field); template ScalarFieldMap(std::enable_if_t field); //! Copy constructor - ScalarFieldMap(const ScalarFieldMap &other) = delete; + ScalarFieldMap(const ScalarFieldMap &other) = default; //! Move constructor - ScalarFieldMap(ScalarFieldMap &&other) noexcept = delete; + ScalarFieldMap(ScalarFieldMap &&other) noexcept = default; //! Destructor virtual ~ScalarFieldMap() noexcept = default; //! Copy assignment operator ScalarFieldMap& operator=(const ScalarFieldMap &other) = delete; //! Move assignment operator ScalarFieldMap& operator=(ScalarFieldMap &&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); const static std::string field_info_root; private: }; /* ---------------------------------------------------------------------- */ template template ScalarFieldMap:: ScalarFieldMap(std::enable_if_t field) :parent(field) { static_assert((isntConst != ConstField), "let the compiler deduce isntConst, this is a SFINAE " "parameter"); this->check_compatibility(); } /* ---------------------------------------------------------------------- */ template template ScalarFieldMap:: ScalarFieldMap(std::enable_if_t field) :parent(field) { static_assert((isConst == ConstField), "let the compiler deduce isntConst, this is a SFINAE " "parameter"); this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string ScalarFieldMap::info_string() const { std::stringstream info; info << "Scalar(" << typeid(T).name() << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ template const std::string ScalarFieldMap::field_info_root{ "Scalar"}; /* ---------------------------------------------------------------------- */ //! member access template template ref_t ScalarFieldMap::operator[](size_type index) { return this->get_ref_to_entry(std::move(index)); } /* ---------------------------------------------------------------------- */ //! member access template typename ScalarFieldMap::reference ScalarFieldMap::operator[](const Ccoord& ccoord) { auto && index = this->collection.get_index(std::move(ccoord)); return this->get_ref_to_entry(std::move(index)); } } // muSpectre #endif /* FIELD_MAP_SCALAR_H */ diff --git a/src/common/field_map_tensor.hh b/src/common/field_map_tensor.hh index 3a9fc87..1d1e1cf 100644 --- a/src/common/field_map_tensor.hh +++ b/src/common/field_map_tensor.hh @@ -1,195 +1,195 @@ /** * 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, ConstField> { public: using Parent = internal::FieldMap; using Ccoord = Ccoord_t; using Sizes = typename SizesByOrder::Sizes; using T_t = Eigen::TensorFixedSize; using value_type = Eigen::TensorMap; using const_reference = Eigen::TensorMap; using reference = std::conditional_t; // since it's a resource handle using size_type = typename Parent::size_type; using pointer = std::unique_ptr; using TypedField = typename Parent::TypedField; using Field = typename TypedField::Parent; using const_iterator = typename Parent::template iterator; using iterator = std::conditional_t< ConstField, 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; template TensorFieldMap(std::enable_if_t field); template TensorFieldMap(std::enable_if_t field); //! Copy constructor - TensorFieldMap(const TensorFieldMap &other) = delete; + TensorFieldMap(const TensorFieldMap &other) = default; //! Move constructor - TensorFieldMap(TensorFieldMap &&other) noexcept = delete; + TensorFieldMap(TensorFieldMap &&other) noexcept = default; //! 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 template TensorFieldMap:: TensorFieldMap(std::enable_if_t field) :Parent(field) { static_assert((isntConst != ConstField), "let the compiler deduce isntConst, this is a SFINAE " "parameter"); this->check_compatibility(); } /* ---------------------------------------------------------------------- */ template template TensorFieldMap:: TensorFieldMap(std::enable_if_t field) :Parent(field) { static_assert((isConst == ConstField), "let the compiler deduce isntConst, this is a SFINAE " "parameter"); 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 e67ab10..f80ac63 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 7f8ae79..48cf11c 100644 --- a/src/materials/material_hyper_elastic1.cc +++ b/src/materials/material_hyper_elastic1.cc @@ -1,70 +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) { + throw 12; 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_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 9133940..fac451e 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,519 +1,657 @@ /** * 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 -#include #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common//utilities.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 NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; using GFieldCollection_t = typename Parent::GFieldCollection_t; 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; using DefaultInternalVariables = std::tuple<>; //! 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); + 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); + 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; - /** - * get a combined range that can be used to iterate over all fields - */ - inline decltype(auto) get_zipped_fields(const StrainField_t & F, - StressField_t & P); - inline decltype(auto) get_zipped_fields(const StrainField_t & F, - StressField_t & P, - TangentField_t & K); - template - inline decltype(auto) get_zipped_fields_worker(Args && ... args); - /** * inheriting classes with internal variables need to overload this function */ decltype(auto) get_internals() { + // the default material has no internal variables + return typename Material::InternalVariables{};} + decltype(auto) get_internals() const { + // the default material has no internal variables return typename Material::InternalVariables{};} 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>; - // TODO: CLEAN THIS - /* These lambdas are 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_small_strain = [this] - (auto && F, auto && P, auto && K, auto && ... internal_variables) { + (auto && Strains, auto && Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, Material::strain_measure)}; + + auto & this_mat = static_cast(*this); + + // Transformation gradient is first in the strains tuple + auto && F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent moduli - std::tie(P, K) = static_cast(*this).evaluate_stress_tangent(std::move(strain), internal_variables...); + Stresses = + std::apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress_tangent(std::move(strain), + internals...);}, + internal_variables); }; auto constitutive_law_finite_strain = [this] - (auto && F, auto && P, auto && K, auto && ... internal_variables) { + (auto && Strains, auto && Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, Material::strain_measure)}; + auto & this_mat = static_cast(*this); + // Transformation gradient is first in the strains tuple + auto && F = std::get<0>(Strains); 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...); + + // TODO: Figure this out: I can't std::move(internals...), + // because if there are no internals, compilation fails with "no + // matching function for call to ‘move()’'. These are tuples of + // lvalue references, so it shouldn't be too bad, but still + // irksome. + + // return value contains a tuple of rvalue_refs to both stress + // and tangent moduli + auto && stress_tgt = + std::apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress_tangent(std::move(strain), + internals...);}, + internal_variables); auto && stress = std::get<0>(stress_tgt); auto && tangent = std::get<1>(stress_tgt); - std::tie(P, K) = MatTB::PK1_stress + Stresses = MatTB::PK1_stress (F, stress, tangent); }; - auto it{this->get_zipped_fields(F, P, K)}; - for (auto && arglist: 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). + iterable_proxy fields{*this, F, P, K}; + for (auto && arglist: fields) { + /** + * arglist is a tuple of three tuples containing only Lvalue + * references (see value_tye in the class definition of + * iterable_proxy::iterator). Tuples contain strains, stresses + * and internal variables, respectively, + */ + //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this + static_assert(std::is_same(std::get<0>(arglist)))>>::value, + "Type mismatch for strain reference, check iterator " + "value_type"); + static_assert(std::is_same(std::get<1>(arglist)))>>::value, + "Type mismatch for stress reference, check iterator" + "value_type"); + static_assert(std::is_same(std::get<1>(arglist)))>>::value, + "Type mismatch for tangent reference, check iterator" + "value_type"); switch (Form) { case Formulation::small_strain: { - std::apply(constitutive_law_small_strain, asStdTuple(arglist)); + std::apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { - std::apply(constitutive_law_finite_strain, asStdTuple(arglist)); + std::apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template - decltype(auto) - MaterialMuSpectre:: - get_zipped_fields(const StrainField_t & F, StressField_t & P) { - typename Material::StrainMap_t Fmap(F);//TODO: think about whether F and P should be std::forward<>()'ed - typename Material::StressMap_t Pmap(P); + template + void MaterialMuSpectre:: + compute_stresses_worker(const StrainField_t & F, + StressField_t & P){ - return this->get_zipped_fields_worker(std::move(Fmap), - std::move(Pmap)); - } + /* These lambdas are executed for every integration point. - /* ---------------------------------------------------------------------- */ - template - decltype(auto) - MaterialMuSpectre:: - get_zipped_fields(const StrainField_t & F, StressField_t & P, - TangentField_t & K) { - typename Material::StrainMap_t Fmap(F); - typename Material::StressMap_t Pmap(P); - typename Material::TangentMap_t Kmap(K); - - return this->get_zipped_fields_worker(std::move(Fmap), - std::move(Pmap), - std::move(Kmap)); - } + F contains the transformation gradient for finite strain calculations and + the infinitesimal strain tensor in small strain problems - /* ---------------------------------------------------------------------- */ - template - template - decltype(auto) - MaterialMuSpectre:: - get_zipped_fields_worker(Args && ... args) { - return std::apply - ([this, &args...] (auto & ... internal_fields) { - return boost::combine(args..., internal_fields...);}, - this->get_internals()); - } + The internal_variables tuple contains whatever internal variables + Material declared (e.g., eigenstrain, strain rate, etc.) + */ + auto constitutive_law_small_strain = [this] + (auto && Strains, auto && Stresses, auto && internal_variables) { + constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; + constexpr StrainMeasure expected_strain_m{ + get_formulation_strain_type(Form, Material::strain_measure)}; + auto & this_mat = static_cast(*this); + // Transformation gradient is first in the strains tuple + auto && F = std::get<0>(Strains); + auto && strain = MatTB::convert_strain(F); + // return value contains a tuple of rvalue_refs to both stress and tangent moduli + auto && sigma = std::get<0>(Stresses); + sigma = + std::apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress(std::move(strain), + internals...);}, + internal_variables); + }; - /** - * Struct with a single variadic function only used to determine the - * tuple type to store the various sub-iterators in - */ - template::NeedTangent Needtgt> - strict tuple_type { - using NeedTangent = - typename MaterialMuSpectre::NeedTangent; - template < + auto constitutive_law_finite_strain = [this] + (auto && Strains, auto && Stresses, auto && internal_variables) { + constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; + constexpr StrainMeasure expected_strain_m{ + get_formulation_strain_type(Form, Material::strain_measure)}; + auto & this_mat = static_cast(*this); + + // Transformation gradient is first in the strains tuple + auto && F = std::get<0>(Strains); + auto && strain = MatTB::convert_strain(F); + + // TODO: Figure this out: I can't std::move(internals...), + // because if there are no internals, compilation fails with "no + // matching function for call to ‘move()’'. These are tuples of + // lvalue references, so it shouldn't be too bad, but still + // irksome. + + // return value contains a tuple of rvalue_refs to both stress + // and tangent moduli + auto && stress = + std::apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress(std::move(strain), + internals...);}, + internal_variables); + auto && P = get<0>(Stresses); + P = MatTB::PK1_stress + (F, stress); + }; + + iterable_proxy fields{*this, F, P}; + + for (auto && arglist: fields) { + /** + * arglist is a tuple of three tuples containing only Lvalue + * references (see value_tye in the class definition of + * iterable_proxy::iterator). Tuples contain strains, stresses + * and internal variables, respectively, + */ + + //auto && stress_tgt = std::get<0>(tuples); + //auto && inputs = std::get<1>(tuples);TODO:clean this + static_assert(std::is_same(std::get<0>(arglist)))>>::value, + "Type mismatch for strain reference, check iterator " + "value_type"); + static_assert(std::is_same(std::get<1>(arglist)))>>::value, + "Type mismatch for stress reference, check iterator" + "value_type"); + + switch (Form) { + case Formulation::small_strain: { + std::apply(constitutive_law_small_strain, std::move(arglist)); + break; + } + case Formulation::finite_strain: { + std::apply(constitutive_law_finite_strain, std::move(arglist)); + break; + } + } + } } - + + /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template - template ::NeedTangent NeedTgt> + template class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; using NeedTangent = typename MaterialMuSpectre::NeedTangent; /** 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. **/ template iterable_proxy(const MaterialMuSpectre & mat, const StrainField_t & F, StressField_t & P, std::enable_if_t & K); template iterable_proxy(const MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P); + using StrainMap_t = typename Material::StrainMap_t; + using StressMap_t = typename Material::StressMap_t; + using TangentMap_t = typename Material::TangentMap_t; using Strain_t = typename Material::Strain_t; using Stress_t = typename Material::Stress_t; using Tangent_t = typename Material::Tangent_t; + using InternalVariables = typename Material::InternalVariables; + using StressFieldTup = std::conditional_t + <(NeedTgt == NeedTangent::yes), + std::tuple, + std::tuple>; + using StressMapTup = std::conditional_t + <(NeedTgt == NeedTangent::yes), + std::tuple, + std::tuple>; + using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes), + std::tuple, + std::tuple>; + //! 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 InternalReferences = MatTB::ReferenceTuple_t; using value_type = - std::conditional_t - <(NeedTgt == NeedTangent::yes), - std::tuple, std::tuple>, - std::tuple, std::tuple>>; + std::tuple, Stress_tTup, InternalReferences>; 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 iterable_proxy & it, bool begin = true) - : it{it}, index{begin ? 0:it.mat.internal_fields.size()}{} + : it{it}, strain_map{it.strain_field}, + stress_map {it.stress_tup}, + index{begin ? 0:it.material.internal_fields.size()}{} //! 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 iterable_proxy & it; + StrainMap_t strain_map; + StressMapTup stress_map; size_t index; private: }; - iterator begin() {return iterator(*this);} - iterator end() {return iterator(*this, false);} + iterator begin() {return std::move(iterator(*this));} + iterator end() {return std::move(iterator(*this, false));} protected: - using StressTup = std::conditional_t - <(NeedTgt == NeedTangent::yes), - std::tuple, - std::tuple>; const MaterialMuSpectre & material; const StrainField_t & strain_field; - auto ZippedFields; - std::tupleStrainField_t - StressTup stress_tup; + InternalVariables & internals; + StressFieldTup stress_tup; private: }; + + /* ---------------------------------------------------------------------- */ + template + template + bool + MaterialMuSpectre::iterable_proxy::iterator:: + operator!=(const iterator & other) const { + return (this->index != other.index); + } + + /* ---------------------------------------------------------------------- */ + template + template + typename MaterialMuSpectre:: + template iterable_proxy:: + iterator & + MaterialMuSpectre::iterable_proxy::iterator:: + operator++() { + this->index++; + return *this; + } + + template + template + typename MaterialMuSpectre:: + template iterable_proxy::iterator:: + value_type + MaterialMuSpectre::iterable_proxy::iterator:: + operator*() { + + const Ccoord_t pixel{ + this->it.material.internal_fields.get_ccoord(this->index)}; + auto && strain = std::make_tuple(this->strain_map[pixel]); + + auto && stresses = + std::apply([&pixel] (auto && ... stress_tgt) { + return std::make_tuple(stress_tgt[pixel]...);}, + this->stress_map); + const auto & internal = this->it.material.get_internals(); + const auto id{index}; + auto && internals = + std::apply([&internal, id] (auto && ... internals) { + return std::make_tuple(internals[id]...);}, + internal); + return std::make_tuple(std::move(strain), + std::move(stresses), + std::move(internals)); + } + // ///* ---------------------------------------------------------------------- */ //template //template ::NeedTangent Tangent> //template //MaterialMuSpectre::iterable_proxy:: //iterable_proxy(const MaterialMuSpectre & mat, // const StrainField_t & F, // StressField_t & P, // std::enable_if_t & K) // : material{mat}, strain_field{F}, stress_tup(P, K), // ZippedFields{std::apply(boost::combine, // std::tuple_cat(std::tuple, // stress_tup, // get_internals()))} //{ // static_assert((Tangent == NeedTangent::yes), // "Explicitly choosing the template parameter 'DoNeedTgt' " // "breaks the SFINAE intention of thi code. Let it be deduced " // "by the compiler."); //} // ///* ---------------------------------------------------------------------- */ //template //template ::NeedTangent Tangent> //template //MaterialMuSpectre::iterable_proxy:: //iterable_proxy(const MaterialMuSpectre & mat, // const StrainField_t & F, // std::enable_if_t & P) // : material{mat}, strain_field{F}, stress_tup(P), // ZippedFields{std::apply(boost::combine, // std::tuple_cat(std::tuple, // stress_tup, // get_internals()))} { // static_assert((Tangent == NeedTangent::no), // "Explicitly choosing the template parameter 'DontNeedTgt' " // "breaks the SFINAE intention of thi code. Let it be deduced " // "by the compiler."); //} // ///* ---------------------------------------------------------------------- */ ////! 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 192a027..c2c246d 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,329 +1,375 @@ /** * 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 #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){} }; + /* ---------------------------------------------------------------------- */ + /** + * Flag used to designate whether the material should compute both stress + * and tangent moduli or only stress + */ + enum class NeedTangent { + yes, // compute both stress and tangent moduli + no}; // compute only stress + /** + * Struct with a single variadic function only used to determine the + * tuple type to store in the various sub-iterators or to return + */ + + template + struct StoredTuple { + template + decltype(auto) operator()(Args && ... args) { + return std::tie(args ... ); + } + }; + + /** + * struct used to determine the exact type of a tuple of references obtained + * when a bunch of iterators over fiel_maps are dereferenced and their + * results are concatenated into a tuple + */ + template + struct ReferenceTuple { + using type = std::tuple; + }; + + /** + * specialisation for tuples + */ + template <> + template + struct ReferenceTuple> { + using type = typename ReferenceTuple::type; + }; + + /** + * helper type for ReferenceTuple + */ + template + using ReferenceTuple_t = typename ReferenceTuple::type; + /* ---------------------------------------------------------------------- */ /** 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 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 { 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 { 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 { 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 { 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) { 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) { Kmap(i,m,j,n) += F(i,r)*get(C,r,m,n,s)*(F(j,s)); } } } } } } 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) { 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 && 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 */