diff --git a/src/common/T4_map_proxy.hh b/src/common/T4_map_proxy.hh new file mode 100644 index 0000000..c01c71c --- /dev/null +++ b/src/common/T4_map_proxy.hh @@ -0,0 +1,112 @@ +/** + * file T4_map_proxy.hh + * + * @author Till Junge + * + * @date 19 Nov 2017 + * + * @brief Map type to allow fourth-order tensor-like maps on 2D matrices + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#include + + +#ifndef T4_MAP_PROXY_H +#define T4_MAP_PROXY_H + +namespace muSpectre { + + /* ---------------------------------------------------------------------- */ + /** Proxy class mapping a fourth-order tensor onto a 2D matrix (in + order to avoid the use of Eigen::Tensor + **/ + template + class T4Map: + public Eigen::MapBase, MapOptions> { + { + public: + typedef MapBase Base; + EIGEN_DENSE_PUBLIC_INTERFACE(T4Map); + + typedef typename Base::PointerType PointerType; + typedef PointerType PointerArgType; + EIGEN_DEVICE_FUNC + inline PointerType cast_to_pointer_type(PointerArgType ptr) { return ptr; } + + EIGEN_DEVICE_FUNC + inline Index innerStride() const + { + return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1; + } + + EIGEN_DEVICE_FUNC + inline Index outerStride() const + { + return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer() + : IsVectorAtCompileTime ? this->size() + : int(Flags)&RowMajorBit ? this->cols() + : this->rows(); + } + + /** Constructor in the fixed-size case. + * + * \param dataPtr pointer to the array to map + * \param stride optional Stride object, passing the strides. + */ + EIGEN_DEVICE_FUNC + explicit inline T4Map(PointerArgType dataPtr, const StrideType& stride = StrideType()) + : Base(cast_to_pointer_type(dataPtr)), m_stride(stride) + { + PlainObjectType::Base::_check_template_params(); + } + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(T4Map); + + /** My accessor to mimick tensorial access + **/ + inline const Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) const { + constexpr auto myColStride{ + (colStride() == 1) ? colStride(): colStride()/Dim}; + constexpr auto myRowStride{ + (rowStride() == 1) ? rowStride(): rowStride()/Dim}; + return this->operator()(i * myRowStride + j * myColStride, + k * myRowStride + l * myColStride); + } + + inline Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) { + constexpr auto myColStride{ + (colStride() == 1) ? colStride(): colStride()/Dim}; + constexpr auto myRowStride{ + (rowStride() == 1) ? rowStride(): rowStride()/Dim}; + return this->operator()(i * myRowStride + j * myColStride, + k * myRowStride + l * myColStride); + } + + protected: + StrideType m_stride; + }; + + +} // muSpectre + +#endif /* T4_MAP_PROXY_H */ diff --git a/src/common/common.cc b/src/common/common.cc index 199f333..b03485d 100644 --- a/src/common/common.cc +++ b/src/common/common.cc @@ -1,51 +1,142 @@ /** * file common.cc * * @author Till Junge * * @date 15 Nov 2017 * * @brief Implementation for common functions * * @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 "common/common.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ std::ostream & operator<<(std::ostream & os, Formulation f) { switch (f) { case Formulation::small_strain: {os << "small_strain"; break;} case Formulation::finite_strain: {os << "finite_strain"; break;} default: throw std::runtime_error ("unknown formulation."); break; } return os; } + /* ---------------------------------------------------------------------- */ + std::ostream & operator<<(std::ostream & os, StressMeasure s) { + switch (s) { + case StressMeasure::Cauchy: {os << "Cauchy"; break;} + case StressMeasure::PK1: {os << "PK1"; break;} + case StressMeasure::PK2: {os << "PK2"; break;} + case StressMeasure::Kirchhoff: {os << "Kirchhoff"; break;} + case StressMeasure::Biot: {os << "Biot"; break;} + case StressMeasure::Mandel: {os << "Mandel"; break;} + default: + throw std::runtime_error + ("a stress measure must be missing"); + break; + } + return os; + } + + /* ---------------------------------------------------------------------- */ + std::ostream & operator<<(std::ostream & os, StrainMeasure s) { + switch (s) { + case StrainMeasure::Gradient: {os << "Gradient"; break;} + case StrainMeasure::Infinitesimal: {os << "Infinitesimal"; break;} + case StrainMeasure::GreenLagrange: {os << "Green-Lagrange"; break;} + case StrainMeasure::Biot: {os << "Biot"; break;} + case StrainMeasure::Log: {os << "Logarithmic"; break;} + case StrainMeasure::Almansi: {os << "Almansi"; break;} + case StrainMeasure::RCauchyGreen: {os << "Right Cauchy-Green"; break;} + case StrainMeasure::LCauchyGreen: {os << "Left Cauchy-Green"; break;} + default: + throw std::runtime_error + ("a strain measure must be missing"); + + } + return os; + } + + /* ---------------------------------------------------------------------- */ + /** Compile-time functions to set the stress and strain measures + stored by mu_spectre depending on the formulation + **/ + constexpr StrainMeasure get_stored_strain_type(Formulation form) { + switch (form) { + case Formulation::finite_strain: { + return StrainMeasure::Gradient; + break; + } + case Formulation::small_strain: { + return StrainMeasure::Infinitesimal; + break; + } + default: + return StrainMeasure::__nostrain__; + break; + } + } + + constexpr StressMeasure get_stored_stress_type(Formulation form) { + switch (form) { + case Formulation::finite_strain: { + return StressMeasure::PK1; + break; + } + case Formulation::small_strain: { + return StressMeasure::Cauchy; + break; + } + default: + return StressMeasure::__nostress__; + break; + } + } + + /* ---------------------------------------------------------------------- */ + constexpr StrainMeasure get_formulation_strain_type(Formulation form, + StrainMeasure expected) { + switch (form) { + case Formulation::finite_strain: { + return expected; + break; + } + case Formulation::small_strain: { + return get_stored_strain_type(form); + break; + } + default: + return StrainMeasure::__nostrain__; + break; + } + } + + } // muSpectre diff --git a/src/common/common.hh b/src/common/common.hh index 24aa1ca..cac9d3f 100644 --- a/src/common/common.hh +++ b/src/common/common.hh @@ -1,90 +1,122 @@ /** * file common.hh * * @author Till Junge * * @date 01 May 2017 * * @brief Small definitions of commonly used types througout µSpectre * * @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 #ifndef COMMON_H #define COMMON_H namespace muSpectre { //! Eigen uses signed integers for dimensions. For consistency, µSpectre uses //! them througout the code using Dim_t = int;// needs to represent -1 for eigen const Dim_t oneD{1}; const Dim_t twoD{2}; const Dim_t threeD{3}; //! Ccoord_t are cell coordinates, i.e. integer coordinates template using Ccoord_t = std::array; template std::ostream & operator << (std::ostream & os, const Ccoord_t & index) { os << "("; for (size_t i = 0; i < dim-1; ++i) { os << index[i] << ", "; } os << index.back() << ")"; return os; } //! Scalar types used for mathematical calculations using Uint = unsigned int; using Int = int; using Real = double; using Complex = std::complex; //! compile-time potentiation required for field-size computations template constexpr I ipow(I base, I exponent) { static_assert(std::is_integral::value, "Type must be integer"); I retval{1}; for (I i = 0; i < exponent; ++i) { retval *= base; } return retval; } //! continuum mechanics flags enum class Formulation{finite_strain, small_strain}; std::ostream & operator<<(std::ostream & os, Formulation f); + /* ---------------------------------------------------------------------- */ + //! Material laws can declare which type of stress measure they provide, + //! and µSpectre will handle conversions + enum class StressMeasure { + Cauchy, PK1, PK2, Kirchhoff, Biot, Mandel, __nostress__}; + std::ostream & operator<<(std::ostream & os, StressMeasure s); + + /* ---------------------------------------------------------------------- */ + //! Material laws can declare which type of strain measure they require and + //! µSpectre will provide it + enum class StrainMeasure { + Gradient, Infinitesimal, GreenLagrange, Biot, Log, Almansi, + RCauchyGreen, LCauchyGreen, __nostrain__}; + std::ostream & operator<<(std::ostream & os, StrainMeasure s); + + /* ---------------------------------------------------------------------- */ + /** Compile-time functions to set the stress and strain measures + stored by mu_spectre depending on the formulation + **/ + constexpr StrainMeasure get_stored_strain_type(Formulation form); + constexpr StressMeasure get_stored_stress_type(Formulation form); + + /* ---------------------------------------------------------------------- */ + /** Compile-time functions to get the stress and strain measures + after they may have been modified by choosing a formulation. + + For instance, a law that expecs a Green-Lagrange strain as input + will get the infinitesimal strain tensor instead in a small + strain computation + **/ + constexpr StrainMeasure get_formulation_strain_type(Formulation form, + StrainMeasure expected); } // muSpectre #ifndef EXPLICITLY_TURNED_ON_CXX17 #include "common/utilities.hh" #endif #endif /* COMMON_H */ diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh index 686081a..2bd0a71 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,194 +1,195 @@ /** * 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" namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ - enum class Map_t{Matrix, Array}; + enum class Map_t{Matrix, Array, T4Matrix}; template struct NameWrapper { const static std::string field_info_root; }; template<> const std::string NameWrapper::field_info_root{"Array"}; template<> const std::string NameWrapper::field_info_root{"Matrix"}; + template<> + const std::string NameWrapper::field_info_root{"T4Matrix"}; /* ---------------------------------------------------------------------- */ template class MatrixLikeFieldMap: public FieldMap { public: using Parent = FieldMap; using Ccoord = Ccoord_t; using value_type = EigenArray; using reference = value_type; // since it's a resource handle using const_reference = const value_type; using size_type = typename Parent::size_type; using pointer = std::unique_ptr; using TypedField = typename Parent::TypedField; using Field = typename TypedField::Parent; using iterator = typename Parent::template iterator; using const_iterator= typename Parent::template iterator; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor MatrixLikeFieldMap() = delete; MatrixLikeFieldMap(Field & field); //! Copy constructor MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = delete; //! Move constructor MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = delete; //! Destructor virtual ~MatrixLikeFieldMap() noexcept = default; //! Copy assignment operator MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete; //! Move assignment operator MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) noexcept = delete; //! give human-readable field map type inline std::string info_string() const override final; //! member access template inline ref operator[](size_type index); inline reference operator[](const Ccoord& ccoord); //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} inline const_iterator cbegin(){return const_iterator(*this);} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);}; inline const_iterator cend(){return const_iterator(*this, false);} protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); const static std::string field_info_root; private: }; template const std::string MatrixLikeFieldMap:: field_info_root{NameWrapper::field_info_root}; /* ---------------------------------------------------------------------- */ template MatrixLikeFieldMap:: MatrixLikeFieldMap(Field & field) :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string MatrixLikeFieldMap:: info_string() const { std::stringstream info; info << field_info_root << "(" << typeid(typename EigenArray::value_type).name() << ", " << EigenArray::RowsAtCompileTime << "x" << EigenArray::ColsAtCompileTime << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template template ref MatrixLikeFieldMap:: operator[](size_type index) { return ref(this->get_ptr_to_entry(index)); } template typename MatrixLikeFieldMap::reference MatrixLikeFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); return reference(this->get_ptr_to_entry(std::move(index))); } /* ---------------------------------------------------------------------- */ template typename MatrixLikeFieldMap::pointer MatrixLikeFieldMap:: ptr_to_val_t(size_type index) { return std::make_unique (this->get_ptr_to_entry(std::move(index))); } } // internal - /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using MatrixFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Matrix>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen array map as iterate template using ArrayFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Array>; } // muSpectre #endif /* FIELD_MAP_MATRIXLIKE_H */ diff --git a/src/materials/CMakeLists.txt b/src/materials/CMakeLists.txt index d22f63a..f80ac63 100644 --- a/src/materials/CMakeLists.txt +++ b/src/materials/CMakeLists.txt @@ -1,9 +1,9 @@ set (materials_SRC material_base.cc - #material_hyper_elastic.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_base.hh b/src/materials/material_base.hh index 519fc86..6ed7e2f 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,116 +1,116 @@ /** * file material_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * @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 "common/common.hh" -#include "common/field_map_tensor.hh" +#include "common/field_map_matrixlike.hh" #include "common/field_collection.hh" #ifndef MATERIAL_BASE_H #define MATERIAL_BASE_H namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for system-wide fields, like stress, strain, etc using GFieldCollection_t = FieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = FieldCollection; - using StressMap_t = TensorFieldMap; + using StressMap_t = MatrixFieldMap; using StrainMap_t = StressMap_t; - using StiffnessMap_t = TensorFieldMap; + using TangentMap_t = MatrixFieldMap; using Ccoord = Ccoord_t; //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) noexcept = delete; //! Destructor virtual ~MaterialBase() noexcept = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) noexcept = delete; //! take responsibility for a pixel identified by its cell coordinates //! TODO: this won't work. for materials with additional info per pixel (as, e.g. for eigenstrain), we need to pass more parameters, so this needs to be a variadic function. void add_pixel(const Ccoord & ccord); //! allocate memory, etc //! TODO: this won't work. for materials with additional info per pixel (see above TODO), we neet to allocate before we know for sure how many pixels the material is responsible for. virtual void initialize(bool stiffness = false) = 0; //! return the materil's name const std::string & get_name() const; //! for static inheritance stuff constexpr static Dim_t sdim() {return DimS;} constexpr static Dim_t mdim() {return DimM;} //! computes stress virtual void compute_stress(const StrainMap_t & F, StressMap_t & P, Formulation form) = 0; //! computes stress and tangent stiffness virtual void compute_stress_stiffness(const StrainMap_t & F, StressMap_t & P, - StiffnessMap_t & K, + TangentMap_t & K, Formulation form) = 0; protected: //! members const std::string name; MFieldCollection_t internal_fields; private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/materials/material_hyperelastic1.cc b/src/materials/material_hyper_elastic1.cc similarity index 67% rename from src/materials/material_hyperelastic1.cc rename to src/materials/material_hyper_elastic1.cc index 486b555..b280702 100644 --- a/src/materials/material_hyperelastic1.cc +++ b/src/materials/material_hyper_elastic1.cc @@ -1,54 +1,66 @@ /** - * file material_hyperelastic1.cc + * 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_hyperelastic1.hh" +#include "materials/material_hyper_elastic1.hh" #include "common/tensor_algebra.hh" +#include namespace muSpectre { /* ---------------------------------------------------------------------- */ template - Materialhyperelastic1::MaterialHyperElastic1(std::string name, + MaterialHyperElastic1::MaterialHyperElastic1(std::string name, Real young, Real poisson) - :young{young}, poisson{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()} + C{lambda*Tensors::outer(Tensors::I2(),Tensors::I2())}// + + //2*mu*Tensors::I4S()} {} /* ---------------------------------------------------------------------- */ template - MaterialHyperElastic1::evaluate_stress(const Strain_t & E, - Stress_t & S) { - S = E.trace()*lambda * Strain_t::Identity() + 2*mu*E; + decltype(auto) + MaterialHyperElastic1::evaluate_stress(const Strain_t & E) { + return E.trace()*lambda * Strain_t::Identity() + 2*mu*E; } + /* ---------------------------------------------------------------------- */ + template + decltype(auto) + MaterialHyperElastic1::evaluate_stress_tangent(const Strain_t & E) { + return std::forward_as_tuple(this->evaluate_stress(E), this->C); + } + + template class MaterialHyperElastic1<2, 2>; + template class MaterialHyperElastic1<2, 3>; + template class MaterialHyperElastic1<3, 3>; + } // muSpectre diff --git a/src/materials/material_hyperelastic1.hh b/src/materials/material_hyper_elastic1.hh similarity index 78% rename from src/materials/material_hyperelastic1.hh rename to src/materials/material_hyper_elastic1.hh index 88e8734..4e11676 100644 --- a/src/materials/material_hyperelastic1.hh +++ b/src/materials/material_hyper_elastic1.hh @@ -1,94 +1,94 @@ /** - * file material_hyperelastic1.hh + * 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_HYPERELASTIC1_H -#define MATERIAL_HYPERELASTIC1_H +#ifndef MATERIAL_HYPER_ELASTIC1_H +#define MATERIAL_HYPER_ELASTIC1_H namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template - class MaterialHyperElastic1: public MaterialMuSpectre + class MaterialHyperElastic1: public MaterialMuSpectre, DimS, DimM> { public: - using Parent = public MaterialMuSpectre; + using Parent = MaterialMuSpectre; // declare what type of strain measure your law takes as input - constexpr static auto strain_measure{MatTB::StrainMeasure::GreenLagrange}; + constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; // declare what type of stress measure your law yields as output - constexpr static auto strain_measure{MatTB::StressMeasure::PK2}; + 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 Strain_t = Eigen::Matrix; using Stress_t = Strain_t; - using Stiffness_t = Eigen::TensorFixedSize - , Eigen::RowMajor>; + using Tangent_t = Eigen::TensorFixedSize + , Eigen::RowMajor>; //! Default constructor MaterialHyperElastic1() = delete; //! Copy constructor MaterialHyperElastic1(const MaterialHyperElastic1 &other) = delete; //! Construct by name, Young's modulus and Poisson's ratio - MaterialMuSpectre(std::string name, Real young, real Poisson); + 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; decltype(auto) evaluate_stress(const Strain_t & E); + decltype(auto) evaluate_stress_tangent(const Strain_t & E); - const Stiffness_t & get_stiffness() const; - + const Tangent_t & get_stiffness() const; protected: const Real young, poisson, lambda, mu; - const Stiffness_t C; + const Tangent_t C; private: }; } // muSpectre -#endif /* MATERIAL_HYPERELASTIC1_H */ +#endif /* MATERIAL_HYPER_ELASTIC1_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 06cfcfb..1a0af29 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,275 +1,358 @@ /** * 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 Stiffness. This class uses the + * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include "stdexcept" -#include "materials/material.hh" +#include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #ifndef MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H namespace muSpectre { //! 'Material' is a CRTP - template - class MaterialMuSpectre: public MaterialBase + template + class MaterialMuSpectre: public MaterialBase { public: - using Parent = public MaterialBase; - using MFieldCollection_t = Parent::MFieldCollection_t; + enum class NeedTangent {yes, no}; + using Parent = MaterialBase; + using MFieldCollection_t = typename Parent::MFieldCollection_t; + using StressMap_t = typename Parent::StressMap_t; + using StrainMap_t = typename Parent::StrainMap_t; + using TangentMap_t = typename Parent::TangentMap_t; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) noexcept = delete; //! Destructor virtual ~MaterialMuSpectre() noexcept = default; //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) noexcept = delete; //! computes stress virtual void compute_stresses(const StrainMap_t & F, - StressMap_t & P, - Formulation form); - //! computes stress and tangent stiffness - virtual void compute_stresses_stiffness(const StrainMap_t & F, + StressMap_t & P, + Formulation form); + //! computes stress and tangent modulus + virtual void compute_stresses_tangent(const StrainMap_t & F, StressMap_t & P, - StiffnessMap_t & K, + TangentMap_t & K, Formulation form); protected: //! computes stress with the formulation available at compile time template inline void compute_stresses_worker(const StrainMap_t & F, StressMap_t & P); //! computes stress with the formulation available at compile time template inline void compute_stresses_worker(const StrainMap_t & F, StressMap_t & P, - StiffnessMap_t & K); + TangentMap_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 stiffness - template + //! that does or does not include tangent moduli + template class iterable_proxy; private: }; /* ---------------------------------------------------------------------- */ - template - void MaterialMuSpectre:: + template + void MaterialMuSpectre:: compute_stresses(const StrainMap_t &F, StressMap_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_stiffness(const StrainMap_t & F, StressMap_t & P, - StiffnessMap_t & K + template + void MaterialMuSpectre:: + compute_stresses_tangent(const StrainMap_t & F, StressMap_t & P, + TangentMap_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 template - void MaterialMuSpectre:: + void MaterialMuSpectre:: compute_stresses_worker(const StrainMap_t & F, - StressMap_t & P, - StiffnessMap_t & K){ + StressMap_t & P, + TangentMap_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), - Material::iterable_proxy, - Material::iterable_proxy>; + IterableSmallStrain_t, + IterableFiniteStrain_t>; - /* This lambda is executed for every intergration point. + constexpr auto stored_strain_m = get_stored_strain_type(Form); + constexpr auto expected_strain_m = get_formulation_strain_type(Form, Material::strain_measure); - The input_tuple contains strain (in the appropriate strain - measure), and stress and potentially tangent moduli in the - appropriate stress measure. + /* This lambda is executed for every integration point. - The internals tuple is potentially empty (for laws without - internal variables, such as e.g., simple hyperelastic - materials), else it contains internal variables in the - appropriate form, as specified in the iterable + 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 fun = [this](auto && input_args...) { - - if (Material::uniform_stiffness) { - + auto constitutive_law = [this](auto && F, auto && ... internal_variables) { + auto && strain = MatTB::convert_strain(F); + // return value contains a tuple of rvalue_refs to both stress and tangent moduli + auto && stress_tgt = this->evaluate_stress_tangent(std::move(strain), internal_variables...); + if (Form == Formulation::small_strain) { + return stress_tgt; + } else { + auto && stress = std::get<0>(stress_tgt); + auto && tangent = std::get<1>(stress_tgt); + return MatTB::PK1_stress + (F, stress, tangent); } }; + + iterable it{this}; + for (auto && tuples: it) { + // the iterator yields a pair of tuples. this first tuple contains + // references to stress and stiffness in the global arrays, the second + // contains references to the deformation gradient and internal variables + // (some of them are const). + auto && stress_tgt = std::get<0>(tuples); + auto && inputs = std::get<1>(tuples); + stress_tgt = std::apply(constitutive_law, std::move(inputs)); + } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain - template - template - class MaterialMuSpectre::iterable_proxy { + template + template::NeedTangent NeedTgt> + class MaterialMuSpectre::iterable_proxy { public: - using Strain_t = typename Material::Strain_t::const_reference; - using Stress_t = typename Material::Stress_t::reference; - using Stiffness_t = typename Material::Stiffness_t::reference; - using value_type = - std::conditional_t, - std::tuple>; - using iterator_category = std::forward_iterator_tag; - //! Default constructor - iterator() = delete; + iterable_proxy() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ - iterator(const MaterialMuSpectre & mat, bool begin = true); + iterable_proxy(const MaterialMuSpectre & mat); + + static constexpr auto strain_measure{StrainMeasure}; + using Strain_t = typename Material::Strain_t::const_reference; + using Stress_t = typename Material::Stress_t::reference; + using Tangent_t = typename Material::Tangent_t::reference; //! Copy constructor - iterator(const iterator &other) = default; + iterable_proxy(const iterable_proxy &other) = default; //! Move constructor - iterator(iterator &&other) noexcept = default; + iterable_proxy(iterable_proxy &&other) noexcept = default; //! Destructor - virtual ~iterator() noexcept = default; + virtual ~iterable_proxy() noexcept = default; //! Copy assignment operator - iterator& operator=(const iterator &other) = default; + iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator - iterator& operator=(iterator &&other) noexcept = default; + iterable_proxy& operator=(iterable_proxy &&other) noexcept = default; + + class iterator + { + public: + using value_type = + std::conditional_t<(NeedTgt == NeedTangent::Yes), + std::tuple, + std::tuple>; + using iterator_category = std::forward_iterator_tag; + + //! Default constructor + iterator() = delete; + + /** Iterator uses the material's internal variables field + collection to iterate selectively over the global fields + (such as the transformation gradient F and first + Piola-Kirchhoff stress P. + **/ + iterator(const MaterialMuSpectre & mat, bool begin = true); + + + //! Copy constructor + iterator(const iterator &other) = default; - //! pre-increment - inline iterator & operator++(); - //! dereference - inline value_type operator*(); - //! inequality - inline bool operator!=(const iterator & other) const; + //! Move constructor + iterator(iterator &&other) noexcept = default; + //! Destructor + virtual ~iterator() noexcept = default; + + //! Copy assignment operator + iterator& operator=(const iterator &other) = default; + + //! Move assignment operator + iterator& operator=(iterator &&other) noexcept = default; + + //! pre-increment + inline iterator & operator++(); + //! dereference + inline value_type operator*(); + //! inequality + inline bool operator!=(const iterator & other) const; + + + protected: + const MaterialMuSpectre & material; + size_t index; + private: + }; protected: const MaterialMuSpectre & material; - size_t index; private: }; /* ---------------------------------------------------------------------- */ - template - template - MaterialMuSpectre::iterator:: - iterator(const MaterialMuSpectre & mat, bool begin) + template + template::NeedTangent Tangent> + MaterialMuSpectre::template iterable_proxy:: + iterable_proxy(const MaterialMuSpectre & mat) + : material{mat} + {} + + /* ---------------------------------------------------------------------- */ + template + template::NeedTangent Tangent> + MaterialMuSpectre::iterable_proxy:: + iterator::iterator(const MaterialMuSpectre & mat, bool begin) : material{mat}, index{begin ? 0:mat.internal_fields.size()} {} /* ---------------------------------------------------------------------- */ //! pre-increment - template - template - typename MaterialMuSpectre::iterator & - MaterialMuSpectre::iterator:: + template + template::NeedTangent Tangent> + + typename MaterialMuSpectre:: + template iterable_proxy::iterator & + + MaterialMuSpectre:: + template iterable_proxy::iterator:: operator++() { ++this->index; + return *this; } /* ---------------------------------------------------------------------- */ //! dereference - template - template - typename MaterialMuSpectre::iterator::value_type - MaterialMuSpectre::iterator:: + template + template::NeedTangent Tangent> + + typename MaterialMuSpectre:: + template iterable_proxy::iterator::value_type + + MaterialMuSpectre:: + template iterable_proxy::iterator:: operator*() { - static_assert(Stiffness, "Template specialisation failed"); + static_assert(Tangent, "Template specialisation failed"); return std::make_tuple - (this->material.internal_fields) + (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 6eba03c..da2a6a6 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,313 +1,306 @@ /** * 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 "common/common.hh" #include "common/tensor_algebra.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){} }; + /* ---------------------------------------------------------------------- */ - //! Material laws can declare which type of stress measure they provide, - //! and µSpectre will handle conversions - enum class StressMeasure { - Cauchy, PK1, PK2, Kirchhoff, Biot, Mandel, __nostress__}; - std::ostream & operator<<(std::ostream & os, StressMeasure s) { - switch (s) { - case StressMeasure::Cauchy: {os << "Cauchy"; break;} - case StressMeasure::PK1: {os << "PK1"; break;} - case StressMeasure::PK2: {os << "PK2"; break;} - case StressMeasure::Kirchhoff: {os << "Kirchhoff"; break;} - case StressMeasure::Biot: {os << "Biot"; break;} - case StressMeasure::Mandel: {os << "Mandel"; break;} - default: - throw MaterialsToolboxError - ("a stress measure must be missing"); - break; + /** 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()); } - return os; - } + + } // internal /* ---------------------------------------------------------------------- */ - //! Material laws can declare which type of strain measure they require and - //! µSpectre will provide it - enum class StrainMeasure { - Gradient, Infinitesimal, GreenLagrange, Biot, Log, Almansi, - RCauchyGreen, LCauchyGreen, __nostrain__}; - std::ostream & operator<<(std::ostream & os, StrainMeasure s) { - switch (s) { - case StrainMeasure::Gradient: {os << "Gradient"; break;} - case StrainMeasure::Infinitesimal: {os << "Infinitesimal"; break;} - case StrainMeasure::GreenLagrange: {os << "Green-Lagrange"; break;} - case StrainMeasure::Biot: {os << "Biot"; break;} - case StrainMeasure::Log: {os << "Logarithmic"; break;} - case StrainMeasure::Almansi: {os << "Almansi"; break;} - case StrainMeasure::RCauchyGreen: {os << "Right Cauchy-Green"; break;} - case StrainMeasure::LCauchyGreen: {os << "Left Cauchy-Green"; break;} - default: - throw MaterialsToolboxError - ("a strain measure must be missing"); - break; - } - return os; - } + //! 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(Stress_t && /*stress*/, Strain_t && /*strain*/) { // 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 + template inline static decltype(auto) - compute(Stress_t && /*stress*/, Stiffness_t && /*stiffness*/, + compute(Stress_t && /*stress*/, Tangent_t && /*stiffness*/, Strain_t && /*strain*/) { // 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(Stress_t && P, Strain_t && /*dummy*/) { 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 { - template + template decltype(auto) - compute(Stress_t && P, Stiffness_t && K, Strain_t && /*dummy*/) { + compute(Stress_t && P, Tangent_t && K, Strain_t && /*dummy*/) { 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(Stress_t && S, Strain_t && F) { 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 { - template + template inline static decltype(auto) - compute(Stress_t && S, Stiffness_t && C, Strain_t && F) { - using T4 = typename Stiffness_t::PlainObject; + compute(Stress_t && S, Tangent_t && C, Strain_t && F) { + using T4 = typename Tangent_t::PlainObject; T4 K; K.setZero(); constexpr int dim{Strain_t::ColsAtCompileTime}; for (int i = 0; i < dim; ++i) { for (int m = 0; m < dim; ++m) { for (int n = 0; n < dim; ++n) { K(i,m,i,n) += S(m,n); for (int j = 0; j < dim; ++j) { for (int r = 0; r < dim; ++r) { for (int s = 0; s < dim; ++s) { K(i,m,j,n) += F(i,r)*C(r,m,n,s)*(F(j,s)); } } } } } } return K; } }; } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template - decltype(auto) PK1_stress(Stress_t && stress, - Strain_t && strain) { + decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) { 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(Stress_t && stress, - Strain_t && strain, - Stiffness_t && stiffness) { + class Stress_t, class Strain_t, class Tangent_t> + decltype(auto) PK1_stress(Strain_t && strain, + Stress_t && stress, + Tangent_t && stiffness) { return internal::PK1_stress::compute(std::move(stress), std::move(strain), std::move(stiffness)); }; /* ---------------------------------------------------------------------- */ /** 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()); }; - /* ---------------------------------------------------------------------- */ - //! function returning expressions for PK2 stress and stiffness - template - decltype(auto) PK2_stress_stiffness(Tens1 && stress, Tens2 && F) { - } - } // MatTB } // muSpectre #endif /* MATERIALS_TOOLBOX_H */