diff --git a/src/common/eigen_tools.hh b/src/common/eigen_tools.hh index d7c03e8..6a01d09 100644 --- a/src/common/eigen_tools.hh +++ b/src/common/eigen_tools.hh @@ -1,205 +1,307 @@ /** * @file eigen_tools.hh * * @author Till Junge * * @date 20 Sep 2017 * * @brief small tools to be used with Eigen * * 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 EIGEN_TOOLS_H #define EIGEN_TOOLS_H #include "common/common.hh" #include +#include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ namespace internal { //! Creates a Eigen::Sizes type for a Tensor defined by an order and dim template struct SizesByOrderHelper { //! type to use using Sizes = typename SizesByOrderHelper::Sizes; }; //! Creates a Eigen::Sizes type for a Tensor defined by an order and dim template struct SizesByOrderHelper<0, dim, dims...> { //! type to use using Sizes = Eigen::Sizes; }; } // internal //! Creates a Eigen::Sizes type for a Tensor defined by an order and dim template struct SizesByOrder { static_assert(order > 0, "works only for order greater than zero"); //! `Eigen::Sizes` using Sizes = typename internal::SizesByOrderHelper::Sizes; }; /* ---------------------------------------------------------------------- */ namespace internal { /* ---------------------------------------------------------------------- */ //! Call a passed lambda with the unpacked sizes as arguments template struct CallSizesHelper { //! applies the call static decltype(auto) call(Fun_t && fun) { static_assert(order > 0, "can't handle empty sizes b)"); return CallSizesHelper::call (fun); } }; /* ---------------------------------------------------------------------- */ template //! Call a passed lambda with the unpacked sizes as arguments struct CallSizesHelper<0, Fun_t, dim, args...> { //! applies the call static decltype(auto) call(Fun_t && fun) { return fun(args...); } }; } // internal /** * takes a lambda and calls it with the proper `Eigen::Sizes` * unpacked as arguments. Is used to call constructors of a * `Eigen::Tensor` or map thereof in a context where the spatial * dimension is templated */ template inline decltype(auto) call_sizes(Fun_t && fun) { static_assert(order > 1, "can't handle empty sizes"); return internal::CallSizesHelper:: call(std::forward(fun)); } //compile-time square root static constexpr Dim_t ct_sqrt(Dim_t res, Dim_t l, Dim_t r){ if(l == r){ return r; } else { const auto mid = (r + l) / 2; if(mid * mid >= res){ return ct_sqrt(res, l, mid); } else { return ct_sqrt(res, mid + 1, r); } } } static constexpr Dim_t ct_sqrt(Dim_t res){ return ct_sqrt(res, 1, res); } namespace EigenCheck { /** * Structure to determine whether an expression can be evaluated * into a `Eigen::Matrix`, `Eigen::Array`, etc. and which helps * determine compile-time size */ template struct is_matrix { //! raw type for testing using T = std::remove_reference_t; //! evaluated test constexpr static bool value{std::is_same::XprKind, Eigen::MatrixXpr>::value}; }; /** * Helper class to check whether an `Eigen::Array` or * `Eigen::Matrix` is statically sized */ template struct is_fixed { //! raw type for testing using T = std::remove_reference_t; //! evaluated test constexpr static bool value{T::SizeAtCompileTime != Eigen::Dynamic}; }; /** * Helper class to check whether an `Eigen::Array` or `Eigen::Matrix` is a * static-size and square. */ template struct is_square { //! raw type for testing using T = std::remove_reference_t; //! true if the object is square and statically sized constexpr static bool value{ (T::RowsAtCompileTime == T::ColsAtCompileTime) && is_fixed::value}; }; /** * computes the dimension from a second order tensor represented * square matrix or array */ template struct tensor_dim { //! raw type for testing using T = std::remove_reference_t; static_assert(is_matrix::value, "The type of t is not understood as an Eigen::Matrix"); static_assert(is_square::value, "t's matrix isn't square"); //! evaluated dimension constexpr static Dim_t value{T::RowsAtCompileTime}; }; //! computes the dimension from a fourth order tensor represented //! by a square matrix template struct tensor_4_dim { //! raw type for testing using T = std::remove_reference_t; static_assert(is_matrix::value, "The type of t is not understood as an Eigen::Matrix"); static_assert(is_square::value, "t's matrix isn't square"); //! evaluated dimension constexpr static Dim_t value{ct_sqrt(T::RowsAtCompileTime)}; static_assert(value*value == T::RowsAtCompileTime, "This is not a fourth-order tensor mapped on a square " "matrix"); }; }; + namespace log_comp { + template + using Mat_t = Eigen::Matrix; + template + using Vec_t = Eigen::Matrix; + + //! This is a static implementation of the explicit determination + //! of log(Tensor) following Jog, C.S. J Elasticity (2008) 93: + //! 141. https://doi.org/10.1007/s10659-008-9169-x + /* ---------------------------------------------------------------------- */ + template + struct Proj { + static inline decltype(auto) compute(const Vec_t & eigs, const Mat_t & T) { + static_assert(dim > 0, "only works for positive dimensions"); + return 1./(eigs(i) -eigs(j))*(T-eigs(j)*Mat_t::Identity()) * Proj::compute(eigs, T); + } + }; + + // catch the case when there's nothing to do + template + struct Proj { + static inline decltype(auto) compute(const Vec_t & eigs, const Mat_t & T) { + static_assert(dim > 0, "only works for positive dimensions"); + return Proj::compute(eigs, T); + } + }; + + // catch the normal tail case + template + struct Proj { + static constexpr Dim_t j{0}; + static inline decltype(auto) compute(const Vec_t & eigs, const Mat_t & T) { + static_assert(dim > 0, "only works for positive dimensions"); + return 1./(eigs(i) -eigs(j))*(T-eigs(j)*Mat_t::Identity()); + } + }; + + // catch the tail case when the last dimension is i + template + struct Proj { + static constexpr Dim_t i{0}; + static constexpr Dim_t j{1}; + static inline decltype(auto) compute(const Vec_t & eigs, const Mat_t & T) { + static_assert(dim > 0, "only works for positive dimensions"); + return 1./(eigs(i) -eigs(j))*(T-eigs(j)*Mat_t::Identity()); + } + }; + + template <> + struct Proj<1, 0, 0> { + static constexpr Dim_t dim{1}; + static constexpr Dim_t i{0}; + static constexpr Dim_t j{0}; + + static inline decltype(auto) compute(const Vec_t & /*eigs*/, const Mat_t & /*T*/) { + return Mat_t::Identity(); + } + }; + + + template + inline decltype(auto) P(const Vec_t & eigs, const Mat_t & T) { + return Proj::compute(eigs, T); + } + + /* ---------------------------------------------------------------------- */ + template + struct Summand { + static inline decltype(auto) compute(const Vec_t & eigs, const Mat_t & T) { + return std::log(eigs(i))*P(eigs, T) + Summand::compute(eigs, T); + } + }; + + template + struct Summand { + static constexpr Dim_t i{0}; + static inline decltype(auto) compute(const Vec_t & eigs, const Mat_t & T) { + return std::log(eigs(i))*P(eigs, T); + } + }; + + template + inline decltype(auto) Sum(const Vec_t & eigs, const Mat_t & T) { + return Summand::compute(eigs, T); + } + + + } // log_comp + + /** + * computes the matrix logarithm efficiently for dim=1, 2, or 3 for + * a diagonizable tensor. For larger tensors, better use the direct + * eigenvalue/vector computation + */ + template + inline decltype(auto) logm(const log_comp::Mat_t& mat) { + using Mat = log_comp::Mat_t; + Eigen::SelfAdjointEigenSolver Solver{}; + Solver.computeDirect(mat, Eigen::EigenvaluesOnly); + return Mat{log_comp::Sum(Solver.eigenvalues(), mat)}; + } } // muSpectre #endif /* EIGEN_TOOLS_H */ diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh index 9fde2e1..403b681 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,393 +1,444 @@ /** * @file materials_toolbox.hh * * @author Till Junge * * @date 02 Nov 2017 * * @brief collection of common continuum mechanics tools * * 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 MATERIALS_TOOLBOX_H #define MATERIALS_TOOLBOX_H #include "common/common.hh" #include "common/tensor_algebra.hh" #include "common/eigen_tools.hh" #include "common/T4_map_proxy.hh" #include +#include #include #include #include #include #include namespace muSpectre { namespace MatTB { /** * thrown when generic materials-related runtime errors occur * (mostly continuum mechanics problems) */ class MaterialsToolboxError:public std::runtime_error{ public: //! constructor explicit MaterialsToolboxError(const std::string& what) :std::runtime_error(what){} //! constructor explicit MaterialsToolboxError(const char * what) :std::runtime_error(what){} }; /* ---------------------------------------------------------------------- */ /** * Flag used to designate whether the material should compute both stress * and tangent moduli or only stress */ enum class NeedTangent { yes, //!< compute both stress and tangent moduli no //!< compute only stress }; /** * struct used to determine the exact type of a tuple of references obtained * when a bunch of iterators over fiel_maps are dereferenced and their * results are concatenated into a tuple */ template struct ReferenceTuple { //! use this type using type = std::tuple; }; /** * specialisation for tuples */ //template <> template struct ReferenceTuple> { //! use this type using type = typename ReferenceTuple::type; }; /** * helper type for ReferenceTuple */ template using ReferenceTuple_t = typename ReferenceTuple::type; /* ---------------------------------------------------------------------- */ namespace internal { /** Structure for functions returning one strain measure as a * function of another **/ template struct ConvertStrain { static_assert((In == StrainMeasure::Gradient) || (In == StrainMeasure::Infinitesimal), "This situation makes me suspect that you are not using " "MatTb as intended. Disable this assert only if you are " "sure about what you are doing."); //! returns the converted strain 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 + E = ¹/₂ (C - I) = ¹/₂ (Fᵀ·F - I) **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t&& F) { return .5*(F.transpose()*F - Strain_t::PlainObject::Identity()); } }; + /* ---------------------------------------------------------------------- */ + /** Specialisation for getting Left Cauchy-Green strain from the + transformation gradient + B = F·Fᵀ = V² + **/ + template <> + struct ConvertStrain { + + //! returns the converted strain + template + inline static decltype(auto) + compute(Strain_t&& F) { + return F*F.transpose(); + } + }; + + /* ---------------------------------------------------------------------- */ + /** Specialisation for getting Right Cauchy-Green strain from the + transformation gradient + C = Fᵀ·F = U² + **/ + template <> + struct ConvertStrain { + + //! returns the converted strain + template + inline static decltype(auto) + compute(Strain_t&& F) { + return F.transpose()*F; + } + }; + + /* ---------------------------------------------------------------------- */ + /** Specialisation for getting logarithmic (Hencky) strain from the + transformation gradient + E₀ = ¹/₂ ln C = ¹/₂ ln (Fᵀ·F) + **/ + template <> + struct ConvertStrain { + + //! returns the converted strain + template + inline static decltype(auto) + compute(Strain_t&& F) { + constexpr Dim_t dim{EigenCheck::tensor_dim::value}; + return (.5*logm(Eigen::Matrix{F.transpose()*F})).eval(); + } + }; + } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning one strain measure as a function of //! another template decltype(auto) convert_strain(Strain_t && strain) { return internal::ConvertStrain::compute(std::move(strain)); }; /* ---------------------------------------------------------------------- */ namespace internal { /** Structure for functions returning PK1 stress from other stress measures **/ template struct PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && /*strain*/, Stress_t && /*stress*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && /*strain*/, Stress_t && /*stress*/, Tangent_t && /*stiffness*/) { // the following test always fails to generate a compile-time error static_assert((StressM == StressMeasure::Cauchy) && (StressM == StressMeasure::PK1), "The requested Stress conversion is not implemented. " "You either made a programming mistake or need to " "implement it as a specialisation of this function. " "See PK2stress for an example."); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress **/ template struct PK1_stress: public PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && /*dummy*/, Stress_t && P) { return std::forward(P); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for the transparent case, where we already have PK1 stress *and* stiffness is given with respect to the transformation gradient **/ template struct PK1_stress: public PK1_stress { //! base class using Parent = PK1_stress; using Parent::compute; //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && /*dummy*/, Stress_t && P, Tangent_t && K) { return std::make_tuple(std::forward(P), std::forward(K)); } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) */ template struct PK1_stress: public PK1_stress { //! returns the converted stress template inline static decltype(auto) compute(Strain_t && F, Stress_t && S) { return F*S; } }; /* ---------------------------------------------------------------------- */ /** * Specialisation for the case where we get material stress (PK2) derived * with respect to Green-Lagrange strain */ template struct PK1_stress: public PK1_stress { //! base class using Parent = PK1_stress; using Parent::compute; //! returns the converted stress and stiffness template inline static decltype(auto) compute(Strain_t && F, Stress_t && S, Tangent_t && C) { using T4 = typename std::remove_reference_t::PlainObject; using Tmap = T4MatMap; T4 K; Tmap Kmap{K.data()}; K.setZero(); for (int i = 0; i < Dim; ++i) { for (int m = 0; m < Dim; ++m) { for (int n = 0; n < Dim; ++n) { get(Kmap,i,m,i,n) += S(m,n); for (int j = 0; j < Dim; ++j) { for (int r = 0; r < Dim; ++r) { for (int s = 0; s < Dim; ++s) { get(Kmap,i,m,j,n) += F(i,r)*get(C,r,m,n,s)*(F(j,s)); } } } } } } auto && P = compute(std::forward(F), std::forward(S)); return std::make_tuple(std::move(P), std::move(K)); } }; } // internal /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) { constexpr Dim_t dim{EigenCheck::tensor_dim::value}; static_assert((dim == EigenCheck::tensor_dim::value), "Stress and strain tensors have differing dimensions"); return internal::PK1_stress::compute (std::forward(strain), std::forward(stress)); }; /* ---------------------------------------------------------------------- */ //! set of functions returning an expression for PK2 stress based on template decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress, Tangent_t && tangent) { constexpr Dim_t dim{EigenCheck::tensor_dim::value}; static_assert((dim == EigenCheck::tensor_dim::value), "Stress and strain tensors have differing dimensions"); static_assert((dim == EigenCheck::tensor_4_dim::value), "Stress and tangent tensors have differing dimensions"); return internal::PK1_stress::compute (std::forward(strain), std::forward(stress), std::forward(tangent)); }; //! static inline implementation of Hooke's law template struct Hooke { /** * compute Lamé's first constant * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_lambda(const Real & young, const Real & poisson) { return young*poisson/((1+poisson)*(1-2*poisson)); } /** * compute Lamé's second constant (i.e., shear modulus) * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_mu(const Real & young, const Real & poisson) { return young/(2*(1+poisson)); } /** * compute the stiffness tensor * @param lambda: Lamé's first constant * @param mu: Lamé's second constant (i.e., shear modulus) */ inline static Eigen::TensorFixedSize> compute_C(const Real & lambda, const Real & mu) { return lambda*Tensors::outer(Tensors::I2(),Tensors::I2()) + 2*mu*Tensors::I4S(); } /** * return stress * @param lambda: First Lamé's constant * @param mu: Second Lamé's constant (i.e. shear modulus) * @param E: Green-Lagrange or small strain tensor */ template inline static decltype(auto) evaluate_stress(const Real & lambda, const Real & mu, s_t && E) { return E.trace()*lambda * Strain_t::Identity() + 2*mu*E; } /** * return stress and tangent stiffness * @param lambda: First Lamé's constant * @param mu: Second Lamé's constant (i.e. shear modulus) * @param E: Green-Lagrange or small strain tensor */ template inline static decltype(auto) evaluate_stress(const Real & lambda, const Real & mu, Tangent_t && C, s_t && E) { return std::make_tuple (std::move(evaluate_stress(lambda, mu, std::move(E))), std::move(C)); } }; } // MatTB } // muSpectre #endif /* MATERIALS_TOOLBOX_H */ diff --git a/tests/test_materials_toolbox.cc b/tests/test_materials_toolbox.cc index 1a88e33..72e3bd2 100644 --- a/tests/test_materials_toolbox.cc +++ b/tests/test_materials_toolbox.cc @@ -1,173 +1,206 @@ /** * @file test_materials_toolbox.cc * * @author Till Junge * * @date 05 Nov 2017 * * @brief Tests for the materials toolbox * * 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. */ #include #include "tests.hh" #include "materials/materials_toolbox.hh" #include "common/T4_map_proxy.hh" #include "common/tensor_algebra.hh" #include "tests/test_goodies.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(materials_toolbox) BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_strain_conversion, Fix, testGoodies::dimlist, Fix){ constexpr Dim_t dim{Fix::dim}; using T2 = Eigen::Matrix; - T2 F; - F.setRandom(); - auto Eref = .5*(F.transpose()*F-T2::Identity()); + T2 F{(T2::Random() -.5*T2::Ones())/20 + T2::Identity()}; - auto E_tb = MatTB::convert_strain (Eigen::Map>(F.data())); - auto error = (Eref-E_tb).norm(); + Real error = (Eref-E_tb).norm(); + BOOST_CHECK_LT(error, tol); + + // checking Left Cauchy-Green + Eref = F*F.transpose(); + E_tb = MatTB::convert_strain(F); + + error = (Eref-E_tb).norm(); + BOOST_CHECK_LT(error, tol); + + // checking Right Cauchy-Green + Eref = F.transpose()*F; + E_tb = MatTB::convert_strain(F); + + error = (Eref-E_tb).norm(); + BOOST_CHECK_LT(error, tol); + + // checking Hencky (logarithmic) strain + Eref = F.transpose()*F; + Eigen::SelfAdjointEigenSolver EigSolv(Eref); + Eref.setZero(); + for (size_t i{0}; i < dim; ++i) { + auto && vec = EigSolv.eigenvectors().col(i); + auto && val = EigSolv.eigenvalues()(i); + Eref += .5*std::log(val) * vec*vec.transpose(); + } + + E_tb = MatTB::convert_strain(F); + + error = (Eref-E_tb).norm(); BOOST_CHECK_LT(error, tol); auto F_tb = MatTB::convert_strain(F); error = (F-F_tb).norm(); BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_tensor_mult_test, Fix, testGoodies::dimlist, Fix) { constexpr Dim_t dim{Fix::dim}; using T4 = T4Mat; T4 A,B, R1, R2; A.setRandom(); B.setRandom(); R1 = A*B; R2.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t a = 0; a < dim; ++a) { for (Dim_t b = 0; b < dim; ++b) { for (Dim_t k = 0; k < dim; ++k) { for (Dim_t l = 0; l < dim; ++l) { get(R2,i,j,k,l) += get(A, i,j,a,b)*get(B, a,b, k,l); } } } } } } auto error{(R1-R2).norm()}; BOOST_CHECK_LT(error, tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_PK1_stress, Fix, testGoodies::dimlist, Fix) { using namespace Matrices; constexpr Dim_t dim{Fix::dim}; using T2 = Eigen::Matrix; using T4 = T4Mat; testGoodies::RandRange rng; T2 F=T2::Identity()*2 ; //F.setRandom(); T2 E_tb = MatTB::convert_strain (Eigen::Map>(F.data())); Real lambda = 3;//rng.randval(1, 2); Real mu = 4;//rng.randval(1,2); T4 J = Itrac(); T2 I = I2(); T4 I4 = Isymm(); T4 C = lambda*J + 2*mu*I4; T2 S = tensmult(C, E_tb); T2 Sref = lambda*E_tb.trace()*I + 2*mu*E_tb; auto error{(Sref-S).norm()}; BOOST_CHECK_LT(error, tol); T4 K = outer_under(I,S) + outer_under(F,I)*C*outer_under(F.transpose(),I); // See Curnier, 2000, "Méthodes numériques en mécanique des solides", p 252 T4 Kref; Real Fkrkr = (F.array()*F.array()).sum(); T2 Fkmkn = F.transpose()*F; T2 Fisjs = F*F.transpose(); Kref.setZero(); for (Dim_t i = 0; i < dim; ++i) { for (Dim_t j = 0; j < dim; ++j) { for (Dim_t m = 0; m < dim; ++m) { for (Dim_t n = 0; n < dim; ++n) { get(Kref, i, m, j, n) = (lambda*((Fkrkr-dim)/2 * I(i,j)*I(m,n) + F(i,m)*F(j,n)) + mu * (I(i,j)*Fkmkn(m,n) + Fisjs(i,j)*I(m,n) - I(i,j) *I(m,n) + F(i,n)*F(j,m))); } } } } error = (Kref-K).norm(); BOOST_CHECK_LT(error, tol); T2 P = MatTB::PK1_stress(F, S); T2 Pref = F*S; error = (P-Pref).norm(); BOOST_CHECK_LT(error, tol); auto && stress_tgt = MatTB::PK1_stress(F, S, C); T2 P_t = std::move(std::get<0>(stress_tgt)); T4 K_t = std::move(std::get<1>(stress_tgt)); error = (P_t-Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_t-Kref).norm(); BOOST_CHECK_LT(error, tol); auto && stress_tgt_trivial = MatTB::PK1_stress(F, P, K); T2 P_u = std::move(std::get<0>(stress_tgt_trivial)); T4 K_u = std::move(std::get<1>(stress_tgt_trivial)); error = (P_u-Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_u-Kref).norm(); BOOST_CHECK_LT(error, tol); T2 P_g; T4 K_g; std::tie(P_g, K_g) = testGoodies::objective_hooke_explicit(lambda, mu, F); error = (P_g-Pref).norm(); BOOST_CHECK_LT(error, tol); error = (K_g-Kref).norm(); BOOST_CHECK_LT(error, tol); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre