Homec4science
Diffusion µSpectre (bugfix/T2542) Browse

muSpectre/src/materials/material_muSpectre_base.hh
fb3a30f7d8f0bugfix/T2542

Last Change
Raw File

material_muSpectre_base.hh
View Options

/**
* file material_muSpectre_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @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 LICENSE
*
* Copyright © 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 MATERIAL_MUSPECTRE_BASE_H
#define MATERIAL_MUSPECTRE_BASE_H
#include "common/common.hh"
#include "materials/material_base.hh"
#include "materials/materials_toolbox.hh"
#include "common/field_collection.hh"
#include "common/field.hh"
#include "common//utilities.hh"
#include <tuple>
#include <type_traits>
#include <iterator>
#include <stdexcept>
namespace muSpectre {
/**
* Forward declaration for factory function
*/
template <Dim_t DimS, Dim_t DimM>
class SystemBase;
template <class Material>
struct MaterialMuSpectre_traits {
};
template <class Material, Dim_t DimS, Dim_t DimM>
class MaterialMuSpectre;
template <>
struct MaterialMuSpectre_traits<void> {
using DefaultInternalVariables = std::tuple<>;
};
//! 'Material' is a CRTP
template <class Material, Dim_t DimS, Dim_t DimM>
class MaterialMuSpectre: public MaterialBase<DimS, DimM>
{
public:
using NeedTangent = MatTB::NeedTangent;
using Parent = MaterialBase<DimS, DimM>;
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<Material>;
//! Default constructor
MaterialMuSpectre() = delete;
//! Construct by name
MaterialMuSpectre(std::string name);
//! Copy constructor
MaterialMuSpectre(const MaterialMuSpectre &other) = delete;
//! Move constructor
MaterialMuSpectre(MaterialMuSpectre &&other) = delete;
//! Destructor
virtual ~MaterialMuSpectre() = default;
//! Factory
template <class... ConstructorArgs>
static Material & make(SystemBase<DimS, DimM> & sys,
ConstructorArgs &&... args);
//! Copy assignment operator
MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete;
//! Move assignment operator
MaterialMuSpectre& operator=(MaterialMuSpectre &&other) = delete;
//* allocate memory, etc
virtual void initialise(bool stiffness = false) override final;
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
//! __attribute__ required by g++-6 and g++-7 because of this bug:
//! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
template <Formulation Form>
inline void compute_stresses_worker(const StrainField_t & F,
StressField_t & P)
__attribute__ ((visibility ("default")));
//! computes stress with the formulation available at compile time
//! __attribute__ required by g++-6 and g++-7 because of this bug:
//! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
template <Formulation Form>
inline void compute_stresses_worker(const StrainField_t & F,
StressField_t & P,
TangentField_t & K)
__attribute__ ((visibility ("default")));
//! 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<NeedTangent need_tgt = NeedTangent::no>
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{};
bool is_initialised{false};
private:
};
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
MaterialMuSpectre<Material, DimS, DimM>::
MaterialMuSpectre(std::string name)
:Parent(name) {
using stress_compatible = typename Material::StressMap_t::
template is_compatible<StressField_t>;
using strain_compatible = typename Material::StrainMap_t::
template is_compatible<StrainField_t>;
using tangent_compatible = typename Material::TangentMap_t::
template is_compatible<TangentField_t>;
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 <class Material, Dim_t DimS, Dim_t DimM>
template <class... ConstructorArgs>
Material & MaterialMuSpectre<Material, DimS, DimM>::
make(SystemBase<DimS, DimM> & sys,
ConstructorArgs && ... args) {
auto mat = std::make_unique<Material>(args...);
auto & mat_ref = *mat;
sys.add_material(std::move(mat));
return mat_ref;
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
initialise(bool /*stiffness*/) {
if (!this->is_initialised) {
this->internal_fields.initialise();
this->is_initialised = true;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses(const StrainField_t &F, StressField_t &P,
Formulation form) {
switch (form) {
case Formulation::finite_strain: {
this->template compute_stresses_worker<Formulation::finite_strain>(F, P);
break;
}
case Formulation::small_strain: {
this->template compute_stresses_worker<Formulation::small_strain>(F, P);
break;
}
default:
throw std::runtime_error("Unknown formulation");
break;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses_tangent(const StrainField_t & F, StressField_t & P,
TangentField_t & K,
Formulation form) {
switch (form) {
case Formulation::finite_strain: {
this->template compute_stresses_worker<Formulation::finite_strain>(F, P, K);
break;
}
case Formulation::small_strain: {
this->template compute_stresses_worker<Formulation::small_strain>(F, P, K);
break;
}
default:
throw std::runtime_error("Unknown formulation");
break;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <Formulation Form>
void MaterialMuSpectre<Material, DimS, DimM>::
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<typename Material::Strain_t>;
using Stresses_t = std::tuple<typename Material::Stress_t,
typename Material::Tangent_t>;
auto constitutive_law_small_strain = [this]
(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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && F = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
// return value contains a tuple of rvalue_refs to both stress and tangent moduli
Stresses =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
};
auto constitutive_law_finite_strain = [this]
(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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && grad = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
// 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 =
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);
Stresses = MatTB::PK1_stress<Material::stress_measure, Material::strain_measure>
(std::move(grad),
std::move(stress),
std::move(tangent));
};
iterable_proxy<NeedTangent::yes> 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<typename Material::Strain_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<0>(arglist)))>>::value,
"Type mismatch for strain reference, check iterator "
"value_type");
static_assert(std::is_same<typename Material::Stress_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<1>(arglist)))>>::value,
"Type mismatch for stress reference, check iterator"
"value_type");
static_assert(std::is_same<typename Material::Tangent_t,
std::remove_reference_t<
decltype(std::get<1>(std::get<1>(arglist)))>>::value,
"Type mismatch for tangent reference, check iterator"
"value_type");
switch (Form) {
case Formulation::small_strain: {
apply(constitutive_law_small_strain, std::move(arglist));
break;
}
case Formulation::finite_strain: {
apply(constitutive_law_finite_strain, std::move(arglist));
break;
}
}
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <Formulation Form>
void MaterialMuSpectre<Material, DimS, DimM>::
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<typename Material::Strain_t>;
using Stresses_t = std::tuple<typename Material::Stress_t>;
auto constitutive_law_small_strain = [this]
(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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && F = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
// return value contains a tuple of rvalue_refs to both stress and tangent moduli
auto && sigma = std::get<0>(Stresses);
sigma =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress(std::move(strain),
internals...);},
internal_variables);
};
auto constitutive_law_finite_strain = [this]
(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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto && F = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(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 =
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<Material::stress_measure, Material::strain_measure>
(F, stress);
};
iterable_proxy<NeedTangent::no> 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<typename Material::Strain_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<0>(arglist)))>>::value,
"Type mismatch for strain reference, check iterator "
"value_type");
static_assert(std::is_same<typename Material::Stress_t,
std::remove_reference_t<
decltype(std::get<0>(std::get<1>(arglist)))>>::value,
"Type mismatch for stress reference, check iterator"
"value_type");
switch (Form) {
case Formulation::small_strain: {
apply(constitutive_law_small_strain, std::move(arglist));
break;
}
case Formulation::finite_strain: {
apply(constitutive_law_finite_strain, std::move(arglist));
break;
}
}
}
}
/* ---------------------------------------------------------------------- */
//! this iterator class is a default for simple laws that just take a strain
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
class MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy {
public:
//! Default constructor
iterable_proxy() = delete;
using NeedTangent =
typename MaterialMuSpectre<Material, DimS, DimM>::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<bool DoNeedTgt=(NeedTgt == NeedTangent::yes)>
iterable_proxy(const MaterialMuSpectre & mat,
const StrainField_t & F,
StressField_t & P,
std::enable_if_t<DoNeedTgt, TangentField_t> & K)
:material(mat), strain_field(F), stress_tup{P,K},
internals(material.internal_variables){};
template<bool DontNeedTgt=(NeedTgt == NeedTangent::no)>
iterable_proxy(const MaterialMuSpectre & mat,
const StrainField_t & F,
std::enable_if_t<DontNeedTgt, StressField_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<StressField_t&, TangentField_t&>,
std::tuple<StressField_t&>>;
using StressMapTup = std::conditional_t
<(NeedTgt == NeedTangent::yes),
std::tuple<StressMap_t, TangentMap_t>,
std::tuple<StressMap_t>>;
using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes),
std::tuple<Stress_t, Tangent_t>,
std::tuple<Stress_t>>;
//! Copy constructor
iterable_proxy(const iterable_proxy &other) = default;
//! Move constructor
iterable_proxy(iterable_proxy &&other) = default;
//! Destructor
virtual ~iterable_proxy() = default;
//! Copy assignment operator
iterable_proxy& operator=(const iterable_proxy &other) = default;
//! Move assignment operator
iterable_proxy& operator=(iterable_proxy &&other) = default;
class iterator
{
public:
using InternalReferences = MatTB::ReferenceTuple_t<InternalVariables>;
using value_type =
std::tuple<std::tuple<Strain_t>, 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) = default;
//! Destructor
virtual ~iterator() = default;
//! Copy assignment operator
iterator& operator=(const iterator &other) = default;
//! Move assignment operator
iterator& operator=(iterator &&other) = 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 <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
bool
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
operator!=(const iterator & other) const {
return (this->index != other.index);
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
typename MaterialMuSpectre<Material, DimS, DimM>::
template iterable_proxy<NeedTgt>::
iterator &
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
operator++() {
this->index++;
return *this;
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgT>
typename MaterialMuSpectre<Material, DimS, DimM>::
template iterable_proxy<NeedTgT>::iterator::
value_type
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgT>::iterator::
operator*() {
const Ccoord_t<DimS> pixel{
this->it.material.internal_fields.get_ccoord(this->index)};
auto && strain = std::make_tuple(this->strain_map[pixel]);
auto && stresses =
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 =
apply([id] (auto && ... internals) {
return std::make_tuple(internals[id]...);},
internal);
return std::make_tuple(std::move(strain),
std::move(stresses),
std::move(internals));
}
} // muSpectre
#endif /* MATERIAL_MUSPECTRE_BASE_H */
c4science · Help