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 <till.junge@epfl.ch>
*
* @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 <unsupported/Eigen/CXX11/Tensor>
+#include <Eigen/Dense>
#include <utility>
#include <type_traits>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
namespace internal {
//! Creates a Eigen::Sizes type for a Tensor defined by an order and dim
template <Dim_t order, Dim_t dim, Dim_t... dims>
struct SizesByOrderHelper {
//! type to use
using Sizes = typename SizesByOrderHelper<order-1, dim, dim, dims...>::Sizes;
};
//! Creates a Eigen::Sizes type for a Tensor defined by an order and dim
template <Dim_t dim, Dim_t... dims>
struct SizesByOrderHelper<0, dim, dims...> {
//! type to use
using Sizes = Eigen::Sizes<dims...>;
};
} // internal
//! Creates a Eigen::Sizes type for a Tensor defined by an order and dim
template <Dim_t order, Dim_t dim>
struct SizesByOrder {
static_assert(order > 0, "works only for order greater than zero");
//! `Eigen::Sizes`
using Sizes = typename
internal::SizesByOrderHelper<order-1, dim, dim>::Sizes;
};
/* ---------------------------------------------------------------------- */
namespace internal {
/* ---------------------------------------------------------------------- */
//! Call a passed lambda with the unpacked sizes as arguments
template<Dim_t order, typename Fun_t, Dim_t dim, Dim_t ... args>
struct CallSizesHelper {
//! applies the call
static decltype(auto) call(Fun_t && fun) {
static_assert(order > 0, "can't handle empty sizes b)");
return CallSizesHelper<order-1, Fun_t, dim, dim, args...>::call
(fun);
}
};
/* ---------------------------------------------------------------------- */
template<typename Fun_t, Dim_t dim, Dim_t ... args>
//! 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<Dim_t order, Dim_t dim, typename Fun_t>
inline decltype(auto) call_sizes(Fun_t && fun) {
static_assert(order > 1, "can't handle empty sizes");
return internal::CallSizesHelper<order-1, Fun_t, dim, dim>::
call(std::forward<Fun_t>(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 <class Derived>
struct is_matrix {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
//! evaluated test
constexpr static bool value{std::is_same<typename Eigen::internal::traits<T>::XprKind,
Eigen::MatrixXpr>::value};
};
/**
* Helper class to check whether an `Eigen::Array` or
* `Eigen::Matrix` is statically sized
*/
template <class Derived>
struct is_fixed {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
//! 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 <class Derived>
struct is_square {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
//! true if the object is square and statically sized
constexpr static bool value{
(T::RowsAtCompileTime == T::ColsAtCompileTime) &&
is_fixed<T>::value};
};
/**
* computes the dimension from a second order tensor represented
* square matrix or array
*/
template <class Derived>
struct tensor_dim {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
static_assert(is_matrix<T>::value, "The type of t is not understood as an Eigen::Matrix");
static_assert(is_square<T>::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 <class Derived>
struct tensor_4_dim {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
static_assert(is_matrix<T>::value, "The type of t is not understood as an Eigen::Matrix");
static_assert(is_square<T>::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 <Dim_t dim>
+ using Mat_t = Eigen::Matrix<Real, dim, dim>;
+ template <Dim_t dim>
+ using Vec_t = Eigen::Matrix<Real, dim, 1>;
+
+ //! 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 <Dim_t dim, Dim_t i, Dim_t j = dim-1>
+ struct Proj {
+ static inline decltype(auto) compute(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ static_assert(dim > 0, "only works for positive dimensions");
+ return 1./(eigs(i) -eigs(j))*(T-eigs(j)*Mat_t<dim>::Identity()) * Proj<dim, i, j-1>::compute(eigs, T);
+ }
+ };
+
+ // catch the case when there's nothing to do
+ template <Dim_t dim, Dim_t other>
+ struct Proj<dim,other,other> {
+ static inline decltype(auto) compute(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ static_assert(dim > 0, "only works for positive dimensions");
+ return Proj<dim, other, other-1>::compute(eigs, T);
+ }
+ };
+
+ // catch the normal tail case
+ template <Dim_t dim, Dim_t i>
+ struct Proj<dim,i,0> {
+ static constexpr Dim_t j{0};
+ static inline decltype(auto) compute(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ static_assert(dim > 0, "only works for positive dimensions");
+ return 1./(eigs(i) -eigs(j))*(T-eigs(j)*Mat_t<dim>::Identity());
+ }
+ };
+
+ // catch the tail case when the last dimension is i
+ template <Dim_t dim>
+ struct Proj<dim,0,1> {
+ static constexpr Dim_t i{0};
+ static constexpr Dim_t j{1};
+ static inline decltype(auto) compute(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ static_assert(dim > 0, "only works for positive dimensions");
+ return 1./(eigs(i) -eigs(j))*(T-eigs(j)*Mat_t<dim>::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<dim> & /*eigs*/, const Mat_t<dim> & /*T*/) {
+ return Mat_t<dim>::Identity();
+ }
+ };
+
+
+ template <Dim_t dim, Dim_t i>
+ inline decltype(auto) P(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ return Proj<dim, i>::compute(eigs, T);
+ }
+
+ /* ---------------------------------------------------------------------- */
+ template <Dim_t dim, Dim_t i = dim-1>
+ struct Summand {
+ static inline decltype(auto) compute(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ return std::log(eigs(i))*P<dim, i>(eigs, T) + Summand<dim, i-1>::compute(eigs, T);
+ }
+ };
+
+ template <Dim_t dim>
+ struct Summand <dim, 0>{
+ static constexpr Dim_t i{0};
+ static inline decltype(auto) compute(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ return std::log(eigs(i))*P<dim, i>(eigs, T);
+ }
+ };
+
+ template <Dim_t dim>
+ inline decltype(auto) Sum(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
+ return Summand<dim>::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 <Dim_t dim>
+ inline decltype(auto) logm(const log_comp::Mat_t<dim>& mat) {
+ using Mat = log_comp::Mat_t<dim>;
+ Eigen::SelfAdjointEigenSolver<Mat> 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 <till.junge@epfl.ch>
*
* @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 <Eigen/Dense>
+#include <unsupported/Eigen/MatrixFunctions>
#include <exception>
#include <sstream>
#include <iostream>
#include <tuple>
#include <type_traits>
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 <class... T>
struct ReferenceTuple {
//! use this type
using type = std::tuple<typename T::reference ...>;
};
/**
* specialisation for tuples
*/
//template <>
template <class... T>
struct ReferenceTuple<std::tuple<T...>> {
//! use this type
using type = typename ReferenceTuple<T...>::type;
};
/**
* helper type for ReferenceTuple
*/
template <class... T>
using ReferenceTuple_t = typename ReferenceTuple<T...>::type;
/* ---------------------------------------------------------------------- */
namespace internal {
/** Structure for functions returning one strain measure as a
* function of another
**/
template <StrainMeasure In, StrainMeasure Out = In>
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 <class Strain_t>
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<Strain_t>(input);
}
};
/* ---------------------------------------------------------------------- */
/** Specialisation for getting Green-Lagrange strain from the
transformation gradient
+ E = ¹/₂ (C - I) = ¹/₂ (Fᵀ·F - I)
**/
template <>
struct ConvertStrain<StrainMeasure::Gradient, StrainMeasure::GreenLagrange> {
//! returns the converted strain
template <class Strain_t>
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<StrainMeasure::Gradient, StrainMeasure::LCauchyGreen> {
+
+ //! returns the converted strain
+ template <class Strain_t>
+ 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<StrainMeasure::Gradient, StrainMeasure::RCauchyGreen> {
+
+ //! returns the converted strain
+ template <class Strain_t>
+ 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<StrainMeasure::Gradient, StrainMeasure::Log> {
+
+ //! returns the converted strain
+ template <class Strain_t>
+ inline static decltype(auto)
+ compute(Strain_t&& F) {
+ constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
+ return (.5*logm(Eigen::Matrix<Real, dim, dim>{F.transpose()*F})).eval();
+ }
+ };
+
} // internal
/* ---------------------------------------------------------------------- */
//! set of functions returning one strain measure as a function of
//! another
template <StrainMeasure In, StrainMeasure Out,
class Strain_t>
decltype(auto) convert_strain(Strain_t && strain) {
return internal::ConvertStrain<In, Out>::compute(std::move(strain));
};
/* ---------------------------------------------------------------------- */
namespace internal {
/** Structure for functions returning PK1 stress from other stress measures
**/
template <Dim_t Dim,
StressMeasure StressM,
StrainMeasure StrainM>
struct PK1_stress {
//! returns the converted stress
template <class Strain_t, class Stress_t>
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<PK1,T1, T2> for an example.");
}
//! returns the converted stress and stiffness
template <class Strain_t, class Stress_t, class Tangent_t>
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<PK1,T1, T2> for an example.");
}
};
/* ---------------------------------------------------------------------- */
/** Specialisation for the transparent case, where we already
have PK1 stress
**/
template <Dim_t Dim, StrainMeasure StrainM>
struct PK1_stress<Dim, StressMeasure::PK1, StrainM>:
public PK1_stress<Dim, StressMeasure::no_stress_,
StrainMeasure::no_strain_> {
//! returns the converted stress
template <class Strain_t, class Stress_t>
inline static decltype(auto)
compute(Strain_t && /*dummy*/, Stress_t && P) {
return std::forward<Stress_t>(P);
}
};
/* ---------------------------------------------------------------------- */
/** Specialisation for the transparent case, where we already have PK1
stress *and* stiffness is given with respect to the transformation
gradient
**/
template <Dim_t Dim>
struct PK1_stress<Dim, StressMeasure::PK1, StrainMeasure::Gradient>:
public PK1_stress<Dim, StressMeasure::PK1,
StrainMeasure::no_strain_> {
//! base class
using Parent = PK1_stress<Dim, StressMeasure::PK1,
StrainMeasure::no_strain_>;
using Parent::compute;
//! returns the converted stress and stiffness
template <class Strain_t, class Stress_t, class Tangent_t>
inline static decltype(auto)
compute(Strain_t && /*dummy*/, Stress_t && P, Tangent_t && K) {
return std::make_tuple(std::forward<Stress_t>(P),
std::forward<Tangent_t>(K));
}
};
/* ---------------------------------------------------------------------- */
/**
* Specialisation for the case where we get material stress (PK2)
*/
template <Dim_t Dim, StrainMeasure StrainM>
struct PK1_stress<Dim, StressMeasure::PK2, StrainM>:
public PK1_stress<Dim, StressMeasure::no_stress_,
StrainMeasure::no_strain_> {
//! returns the converted stress
template <class Strain_t, class Stress_t>
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 <Dim_t Dim>
struct PK1_stress<Dim, StressMeasure::PK2, StrainMeasure::GreenLagrange>:
public PK1_stress<Dim, StressMeasure::PK2,
StrainMeasure::no_strain_> {
//! base class
using Parent = PK1_stress<Dim, StressMeasure::PK2,
StrainMeasure::no_strain_>;
using Parent::compute;
//! returns the converted stress and stiffness
template <class Strain_t, class Stress_t, class Tangent_t>
inline static decltype(auto)
compute(Strain_t && F, Stress_t && S, Tangent_t && C) {
using T4 = typename std::remove_reference_t<Tangent_t>::PlainObject;
using Tmap = T4MatMap<Real, Dim>;
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<Strain_t>(F),
std::forward<Stress_t>(S));
return std::make_tuple(std::move(P), std::move(K));
}
};
} // internal
/* ---------------------------------------------------------------------- */
//! set of functions returning an expression for PK2 stress based on
template <StressMeasure StressM, StrainMeasure StrainM,
class Stress_t, class Strain_t>
decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) {
constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
static_assert((dim == EigenCheck::tensor_dim<Stress_t>::value),
"Stress and strain tensors have differing dimensions");
return internal::PK1_stress<dim, StressM, StrainM>::compute
(std::forward<Strain_t>(strain),
std::forward<Stress_t>(stress));
};
/* ---------------------------------------------------------------------- */
//! set of functions returning an expression for PK2 stress based on
template <StressMeasure StressM, StrainMeasure StrainM,
class Stress_t, class Strain_t, class Tangent_t>
decltype(auto) PK1_stress(Strain_t && strain,
Stress_t && stress,
Tangent_t && tangent) {
constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
static_assert((dim == EigenCheck::tensor_dim<Stress_t>::value),
"Stress and strain tensors have differing dimensions");
static_assert((dim == EigenCheck::tensor_4_dim<Tangent_t>::value),
"Stress and tangent tensors have differing dimensions");
return internal::PK1_stress<dim, StressM, StrainM>::compute
(std::forward<Strain_t>(strain),
std::forward<Stress_t>(stress),
std::forward<Tangent_t>(tangent));
};
//! static inline implementation of Hooke's law
template <Dim_t Dim, class Strain_t, class Tangent_t>
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<Real, Eigen::Sizes<Dim, Dim, Dim, Dim>>
compute_C(const Real & lambda, const Real & mu) {
return lambda*Tensors::outer<Dim>(Tensors::I2<Dim>(),Tensors::I2<Dim>()) +
2*mu*Tensors::I4S<Dim>();
}
/**
* 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 <class s_t>
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 <class s_t>
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 <till.junge@altermail.ch>
*
* @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 <Eigen/Dense>
#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<Real, dim, dim>;
- 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<StrainMeasure::Gradient,
+ // checking Green-Lagrange
+ T2 Eref = .5*(F.transpose()*F-T2::Identity());
+
+ T2 E_tb = MatTB::convert_strain<StrainMeasure::Gradient,
StrainMeasure::GreenLagrange>
(Eigen::Map<Eigen::Matrix<Real, dim, dim>>(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<StrainMeasure::Gradient,
+ StrainMeasure::LCauchyGreen>(F);
+
+ error = (Eref-E_tb).norm();
+ BOOST_CHECK_LT(error, tol);
+
+ // checking Right Cauchy-Green
+ Eref = F.transpose()*F;
+ E_tb = MatTB::convert_strain<StrainMeasure::Gradient,
+ StrainMeasure::RCauchyGreen>(F);
+
+ error = (Eref-E_tb).norm();
+ BOOST_CHECK_LT(error, tol);
+
+ // checking Hencky (logarithmic) strain
+ Eref = F.transpose()*F;
+ Eigen::SelfAdjointEigenSolver<T2> 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<StrainMeasure::Gradient,
+ StrainMeasure::Log>(F);
+
+ error = (Eref-E_tb).norm();
BOOST_CHECK_LT(error, tol);
auto F_tb = MatTB::convert_strain<StrainMeasure::Gradient, StrainMeasure::Gradient>(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<Real, dim>;
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<Real, dim, dim>;
using T4 = T4Mat<Real, dim>;
testGoodies::RandRange<Real> rng;
T2 F=T2::Identity()*2 ;
//F.setRandom();
T2 E_tb = MatTB::convert_strain<StrainMeasure::Gradient, StrainMeasure::GreenLagrange>
(Eigen::Map<Eigen::Matrix<Real, dim, dim>>(F.data()));
Real lambda = 3;//rng.randval(1, 2);
Real mu = 4;//rng.randval(1,2);
T4 J = Itrac<dim>();
T2 I = I2<dim>();
T4 I4 = Isymm<dim>();
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<StressMeasure::PK2, StrainMeasure::GreenLagrange>(F, S);
T2 Pref = F*S;
error = (P-Pref).norm();
BOOST_CHECK_LT(error, tol);
auto && stress_tgt = MatTB::PK1_stress<StressMeasure::PK2, StrainMeasure::GreenLagrange>(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<StressMeasure::PK1, StrainMeasure::Gradient>(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