diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh index 6aa1db7..a1f32ff 100644 --- a/src/common/field_map_base.hh +++ b/src/common/field_map_base.hh @@ -1,556 +1,572 @@ /** * 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) = default; //! Move constructor 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 const_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*(); + //! dereference + inline const_value_type operator*() const; //! 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 pointer get_ptr_to_entry(size_t index); inline const T* get_ptr_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 + typename FieldMap:: + template iterator::value_type FieldMap::iterator:: operator*() { return this->fieldmap.template operator[](this->index); } + /* ---------------------------------------------------------------------- */ + //! dereference + template + template + typename FieldMap:: + template iterator::const_value_type + FieldMap::iterator:: + operator*() const { + 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 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)); } } // internal } // muSpectre #endif /* FIELD_MAP_H */ diff --git a/src/materials/material_hyper_elastic1.hh b/src/materials/material_hyper_elastic1.hh index ca98bf5..b8656f5 100644 --- a/src/materials/material_hyper_elastic1.hh +++ b/src/materials/material_hyper_elastic1.hh @@ -1,133 +1,133 @@ /** * 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 { template class MaterialHyperElastic1; template <> template struct MaterialMuSpectre_traits>: public MaterialMuSpectre_traits { using Parent = MaterialMuSpectre_traits; using InternalVariables = typename Parent::DefaultInternalVariables; }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template - class MaterialHyperElastic1: public MaterialMuSpectre, DimS, DimM> + class MaterialHyperElastic1: + public MaterialMuSpectre, DimS, DimM> { public: using Parent = MaterialMuSpectre; using NeedTangent = typename Parent::NeedTangent; 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 = MatrixFieldMap; 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 Stiffness_t = Eigen::TensorFixedSize >; using InternalVariables = typename Parent::DefaultInternalVariables; //! 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 inline decltype(auto) evaluate_stress(s_t && E); template inline 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: }; /* ---------------------------------------------------------------------- */ 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()))); } } // muSpectre #endif /* MATERIAL_HYPER_ELASTIC1_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 802fa76..65bc683 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,691 +1,704 @@ /** * 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 "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 { template struct MaterialMuSpectre_traits { }; template class MaterialMuSpectre; template <> struct MaterialMuSpectre_traits { using DefaultInternalVariables = std::tuple<>; }; //! 'Material' is a CRTP template class MaterialMuSpectre: public MaterialBase { public: 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<>; using traits = MaterialMuSpectre_traits; //! 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; //* allocate memory, etc virtual void initialise(bool stiffness = false); using Parent::compute_stresses; using Parent::compute_stresses_tangent; //! 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; /** * 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{};} typename traits::InternalVariables internal_variables{}; 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:: initialise(bool /*stiffness*/) { this->internal_fields.initialise(); } /* ---------------------------------------------------------------------- */ 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){ /* 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.) */ + using Strains_t = std::tuple; + using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] - (auto && Strains, auto && Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t 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 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 && Strains, auto && Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t 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_tgt = + 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); + //auto && stress = std::get<0>(stress_tgt); + //auto && tangent = std::get<1>(stress_tgt); + Eigen::Matrix stress = std::get<0>(stress_tgt); + Eigen::Matrix tangent = std::get<1>(stress_tgt); Stresses = MatTB::PK1_stress - (F, stress, tangent); + (std::move(F), std::move(stress), std::move(tangent)); + std::cout << "bla" << std::endl; }; 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, std::move(arglist)); break; } case Formulation::finite_strain: { std::apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P){ /* 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.) */ + + using Strains_t = std::tuple; + using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] - (auto && Strains, auto && Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t 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); }; auto constitutive_law_finite_strain = [this] - (auto && Strains, auto && Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t && 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)}; + size_t StrainsSize = sizeof(Strains_t); + std::cout << "Size of strain type = " << StrainsSize << ", map_size = " + << sizeof(typename Material::Strain_t) << " ptr_size = " + << sizeof(typename Material::Strain_t::PointerType) << std::endl; 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 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) :material(mat), strain_field(F), stress_tup{P,K}, internals(material.internal_variables){}; template iterable_proxy(const MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) :material(mat), strain_field(F), stress_tup{P}, internals(material.internal_variables){}; 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::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}, 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 std::move(iterator(*this));} iterator end() {return std::move(iterator(*this, false));} protected: const MaterialMuSpectre & material; const StrainField_t & strain_field; StressFieldTup stress_tup; const InternalVariables & internals; 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/tests/test_material_hyper_elastic1.cc b/tests/test_material_hyper_elastic1.cc index de2cd2d..f81d46f 100644 --- a/tests/test_material_hyper_elastic1.cc +++ b/tests/test_material_hyper_elastic1.cc @@ -1,196 +1,216 @@ /** * file test_material_hyper_elastic1.cc * * @author Till Junge * * @date 28 Nov 2017 * * @brief Tests for the large-strain, objective Hooke's law, implemented in * the convenient strategy (i.e., using MaterialMuSpectre) * * @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 "materials/material_hyper_elastic1.hh" #include "tests.hh" #include "common/test_goodies.hh" #include "common/field_collection.hh" #include "common/iterators.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_hyper_elastic_1); template struct MaterialFixture { using Mat_t = MaterialHyperElastic1; constexpr static Real lambda{2}, mu{1.5}; - MaterialFixture():mat("Name", lambda, mu){}; + constexpr static Real young{mu*(3*lambda + 2*mu)/(lambda + mu)}; + constexpr static Real poisson{lambda/(2*(lambda + mu))}; + MaterialFixture():mat("Name", young, poisson){}; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; Mat_t mat; }; using mat_list = boost::mpl::list, MaterialFixture, MaterialFixture>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_add_pixel, Fix, mat_list, Fix) { auto & mat{Fix::mat}; constexpr Dim_t sdim{Fix::sdim}; testGoodies::RandRange rng;; const Dim_t nb_pixel{7}, box_size{17}; using Ccoord = Ccoord_t; for (Dim_t i = 0; i < nb_pixel; ++i) { Ccoord c; for (Dim_t j = 0; j < sdim; ++j) { c[j] = rng.randval(0, box_size); } BOOST_CHECK_NO_THROW(mat.add_pixel(c)); } BOOST_CHECK_NO_THROW(mat.initialise()); } template struct MaterialFixtureFilled: public MaterialFixture { using Mat_t = typename MaterialFixture::Mat_t; constexpr static Dim_t box_size{3}; MaterialFixtureFilled():MaterialFixture(){ using Ccoord = Ccoord_t; Ccoord cube{CcoordOps::get_cube(box_size)}; CcoordOps::Pixels pixels(cube); for (auto pixel: pixels) { this->mat.add_pixel(pixel); } this->mat.initialise(); }; }; using mat_fill = boost::mpl::list, MaterialFixtureFilled, MaterialFixtureFilled>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_law, Fix, mat_fill, Fix) { constexpr auto cube{CcoordOps::get_cube(Fix::box_size)}; auto & mat{Fix::mat}; using FC_t = FieldCollection; FC_t globalfields; auto & F{make_field ("Transformation Gradient", globalfields)}; auto & P1 = make_field ("Nominal Stress1", globalfields); // to be computed alone auto & P2 = make_field ("Nominal Stress2", globalfields); // to be computed with tangent auto & K = make_field ("Tangent Moduli", globalfields); // to be computed with tangent auto & Pr = make_field ("Nominal Stress reference", globalfields); auto & Kr = make_field ("Tangent Moduli reference", globalfields); // to be computed with tangent globalfields.initialise(cube); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); { // block to contain not-constant gradient map typename Fix::Mat_t::StressMap_t grad_map (globalfields["Transformation Gradient"]); for (auto F_: grad_map) { F_.setRandom(); } grad_map[0] = grad_map[0].Identity(); // identifiable gradients for debug grad_map[1] = 1.2*grad_map[1].Identity(); // ditto - size_t counter{0}; - for (auto F_: grad_map) { - std::cout << "F(" << counter++ << ")" << std::endl << F_ << std::endl; - } } //compute stresses using material mat.compute_stresses(globalfields["Transformation Gradient"], globalfields["Nominal Stress1"], Formulation::finite_strain); //compute stresses and tangent moduli using material BOOST_CHECK_THROW (mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Nominal Stress2"], Formulation::finite_strain), std::runtime_error); - //mat.compute_stresses_tangent(globalfields["Transformation Gradient"], - // globalfields["Nominal Stress2"], - // globalfields["Tangent Moduli"], - // Formulation::finite_strain); + mat.compute_stresses_tangent(globalfields["Transformation Gradient"], + globalfields["Nominal Stress2"], + globalfields["Tangent Moduli"], + Formulation::finite_strain); - // auto zip_it = akantu::zip - // (typename Fix::Mat_t::StrainMap_t(globalfields["Transformation Gradient"]), - // typename Fix::Mat_t::StressMap_t(globalfields["Nominal Stress reference"]), - // typename Fix::Mat_t::TangentMap_t(globalfields["Tangent Moduli reference"])); - // for (auto tup: zip_it) { typename Fix::Mat_t::StrainMap_t Fmap(globalfields["Transformation Gradient"]); - typename Fix::Mat_t::StressMap_t Pmap(globalfields["Nominal Stress reference"]); - typename Fix::Mat_t::TangentMap_t Kmap(globalfields["Tangent Moduli reference"]); - for (size_t i = 0; i < globalfields.size(); ++i) { - const typename Fix::Mat_t::Strain_t F_ = Fmap[i];//std::get<0>(tup); - std::cout << "F" << std::endl << F_ << std::endl; - typename Fix::Mat_t::Stress_t P_ = Pmap[i];//std::get<1>(tup); - typename Fix::Mat_t::Tangent_t K_ = Kmap[i];//std::get<2>(tup); + typename Fix::Mat_t::StressMap_t Pmap_ref(globalfields["Nominal Stress reference"]); + typename Fix::Mat_t::TangentMap_t Kmap_ref(globalfields["Tangent Moduli reference"]); + + for (auto tup: boost::combine(Fmap, Pmap_ref, Kmap_ref)) { + auto F_ = boost::get<0>(tup); + auto P_ = boost::get<1>(tup); + auto K_ = boost::get<2>(tup); std::tie(P_,K_) = testGoodies::objective_hooke_explicit (Fix::lambda, Fix::mu, F_); } + typename Fix::Mat_t::StressMap_t Pmap_1(globalfields["Nominal Stress1"]); + for (auto tup: boost::combine(Pmap_ref, Pmap_1)) { + auto P_r = boost::get<0>(tup); + auto P_1 = boost::get<1>(tup); + Real error = (P_r - P_1).norm(); + BOOST_CHECK_LT(error, tol); + } + + typename Fix::Mat_t::StressMap_t Pmap_2(globalfields["Nominal Stress2"]); + typename Fix::Mat_t::TangentMap_t Kmap(globalfields["Tangent Moduli"]); + for (auto tup: boost::combine(Pmap_ref, Pmap_2, Kmap_ref, Kmap)) { + auto P_r = boost::get<0>(tup); + auto P = boost::get<1>(tup); + Real error = (P_r - P).norm(); + std::cout << "P_r =" << std::endl << P_r << std::endl; + std::cout << "P =" << std::endl << P << std::endl; + BOOST_CHECK_LT(error, tol); + + auto K_r = boost::get<2>(tup); + auto K = boost::get<3>(tup); + std::cout << "K_r =" << std::endl << K_r << std::endl; + std::cout << "K =" << std::endl << K << std::endl; + error = (K_r - K).norm(); + BOOST_CHECK_LT(error, tol); + } + } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre