Add vanishing strain wrapper (for plane strain) (#317)

Author: henrij22

CHANGELOG.md

@@ -43,7 +43,7 @@ SPDX-License-Identifier: LGPL-3.0-or-later
     - Python:
         - Bindings for Enums can now be done conveniently with the `ENUM_BINDINGS` macro.
         - The finite element functions `calculateMatrix`, `calculateVector`, and `calculateScalar` now directly accept the affordances.
-        - The assembler bindings now also accept affordances and `DBCOption`, and they are also renamed to simply `matrix`, `vector`
+        - The assembler bindings now also accept affordances and `DBCOption`, and they are also renamed  to simply `matrix`, `vector`
           and `scalar`.
         - The assemblers also export the binding functions to bind the assemblers.
 - Added a new class `AssemblerManipulator` that wraps an existing assembler
@@ -55,6 +55,7 @@ SPDX-License-Identifier: LGPL-3.0-or-later
 - Add a truss element ([#302](https://github.com/ikarus-project/ikarus/pull/302))
 - Add an About Ikarus page in the documentation ([#291](https://github.com/ikarus-project/ikarus/pull/291))
 - Add new class `Vtk::Writer`, which implements some convenience methods over the existing `dune-vtk` module ([#309](https://github.com/ikarus-project/ikarus/pull/309))
+- Add `VanishingStrain` material (useful for example for plane strain case), also refactor the constructor of `LinearElastic` to take any linear material law ([#317](https://github.com/ikarus-project/ikarus/pull/317))
 
 ## Release v0.4 (Ganymede)
 

ikarus/finiteelements/mechanics/linearelastic.hh

@@ -26,18 +26,21 @@
 
 namespace Ikarus {
 
-template <typename PreFE, typename FE>
+template <typename PreFE, typename FE, typename PRE>
 class LinearElastic;
 
 /**
  * \brief A PreFE struct for linear elastic elements.
+ * \tparam MAT Type of the material.
  */
+template <Concepts::GeometricallyLinearMaterial MAT>
 struct LinearElasticPre
 {
-  YoungsModulusAndPoissonsRatio material;
+  using Material = MAT;
+  MAT material;
 
   template <typename PreFE, typename FE>
-  using Skill = LinearElastic<PreFE, FE>;
+  using Skill = LinearElastic<PreFE, FE, LinearElasticPre>;
 };
 
 /**
@@ -48,7 +51,7 @@ struct LinearElasticPre
  * \tparam PreFE The type of the  pre finite element.
  * \tparam FE The type of the finite element.
  */
-template <typename PreFE, typename FE>
+template <typename PreFE, typename FE, typename PRE>
 class LinearElastic : public ResultTypeBase<ResultTypes::linearStress>
 {
 public:
@@ -61,7 +64,8 @@ public:
   using Geometry  = typename Traits::Geometry;
   using GridView  = typename Traits::GridView;
   using Element   = typename Traits::Element;
-  using Pre       = LinearElasticPre;
+  using Material  = PRE::Material;
+  using Pre       = PRE;
 
   static constexpr int myDim = Traits::mydim;
   using LocalBasisType       = decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
@@ -115,7 +119,8 @@ public:
     return uFunction;
   }
   /**
-   * \brief Gets the strain function for the given Requirement and optional displacement vector.
+   * \brief Gets the strain function for the given Requirement and optional di
+   splacement vector.
    *
    * \tparam ScalarType The scalar type for the strain vector.
    * \param par The Requirement object.
@@ -135,10 +140,8 @@ public:
    * \return The material tangent matrix.
    */
   auto materialTangent() const {
-    if constexpr (myDim == 2)
-      return planeStressLinearElasticMaterialTangent(mat_.emodul, mat_.nu);
-    else if constexpr (myDim == 3)
-      return linearElasticMaterialTangent3D(mat_.emodul, mat_.nu);
+    // Since that material is independent of the strains, a zero strain is passed here
+    return mat_.template tangentModuli<StrainTags::linear, true>(Eigen::Vector<double, 6>::Zero());
   }
 
   /**
@@ -187,7 +190,7 @@ private:
 
   std::shared_ptr<const Geometry> geo_;
   Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>> localBasis_;
-  YoungsModulusAndPoissonsRatio mat_;
+  Material mat_;
   size_t numberOfNodes_{0};
   int order_{};
 
@@ -257,11 +260,13 @@ protected:
 
 /**
  * \brief A helper function to create a linear elastic pre finite element.
- * \param mat Material parameters for the linear elastic element.
+ * \tparam MAT Type of the material.
+ * \param mat Material parameters for the non-linear elastic element.
  * \return A linear elastic pre finite element.
  */
-auto linearElastic(const YoungsModulusAndPoissonsRatio& mat) {
-  LinearElasticPre pre(mat);
+template <typename MAT>
+auto linearElastic(const MAT& mat) {
+  LinearElasticPre<MAT> pre(mat);
 
   return pre;
 }

ikarus/finiteelements/mechanics/materials.hh

@@ -13,4 +13,5 @@
 #include <ikarus/finiteelements/mechanics/materials/neohooke.hh>
 #include <ikarus/finiteelements/mechanics/materials/svk.hh>
 #include <ikarus/finiteelements/mechanics/materials/tags.hh>
+#include <ikarus/finiteelements/mechanics/materials/vanishingstrain.hh>
 #include <ikarus/finiteelements/mechanics/materials/vanishingstress.hh>

ikarus/finiteelements/mechanics/materials/CMakeLists.txt

@@ -10,5 +10,7 @@ install(
         svk.hh
         tags.hh
         vanishingstress.hh
+        vanishingstrain.hh
+        vanishinghelpers.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics/materials
 )

ikarus/finiteelements/mechanics/materials/interface.hh

@@ -24,6 +24,9 @@ struct Material;
 
 template <auto stressIndexPair, typename MImpl>
 struct VanishingStress;
+
+template <auto strainIndexPair, typename MImpl>
+struct VanishingStrain;
 #endif
 
 /**
@@ -78,9 +81,10 @@ struct Material
   using MaterialImpl = MI; ///< Type of material implementation
 
   /**
-   * \brief Static constant for determining if the material has vanishing stress components (is reduced).
+   * \brief Static constant for determining if the material has vanishing stress or strain components (is reduced).
    */
-  static constexpr bool isReduced = traits::isSpecializationNonTypeAndTypes<VanishingStress, MaterialImpl>::value;
+  static constexpr bool isReduced = traits::isSpecializationNonTypeAndTypes<VanishingStress, MaterialImpl>::value or
+                                    traits::isSpecializationNonTypeAndTypes<VanishingStrain, MaterialImpl>::value;
 
   /**
    * \brief Const accessor to the underlying material (CRTP).
@@ -101,7 +105,21 @@ struct Material
    *
    * \return Name of the material.
    */
-  [[nodiscard]] constexpr std::string name() const { return impl().nameImpl(); }
+  [[nodiscard]] constexpr static std::string name() { return MI::nameImpl(); }
+
+  /**
+   * \brief Returns the material parameters stored in the implemented material.
+   * \return Material parameter.
+   */
+  [[nodiscard]] auto materialParameters() const { return impl().materialParametersImpl(); }
+
+  /**
+   * \brief This factor denotes the differences between the returned stresses and moduli and the passed strain
+   * \details For neoHooke the inserted quantity is $C$ the Green-Lagrangian strain tensor, the function relation
+   * between the energy and the stresses is $S = 1 \dfrac{\partial\psi(E)}{\partial E}$. This factor is the pre factor,
+   * which is the difference to the actual derivative, which is written here
+   */
+  static constexpr double derivativeFactor = MI::derivativeFactorImpl;
 
   /**
    * \brief Return the stored potential energy of the material.

ikarus/finiteelements/mechanics/materials/linearelasticity.hh

@@ -34,38 +34,39 @@ namespace Ikarus {
 template <typename ST>
 struct LinearElasticityT : Material<LinearElasticityT<ST>>
 {
-  [[nodiscard]] constexpr std::string nameImpl() const noexcept { return "LinearElasticity"; }
-
   using ScalarType = ST;
   using Base       = StVenantKirchhoffT<ScalarType>;
 
+  using field_type                    = ScalarType;
+  static constexpr int worldDimension = 3;
+  using StrainMatrix                  = Eigen::Matrix<ScalarType, worldDimension, worldDimension>;
+  using StressMatrix                  = StrainMatrix;
+  using MaterialParameters            = typename Base::MaterialParameters;
+
+  static constexpr auto strainTag              = StrainTags::linear;
+  static constexpr auto stressTag              = StressTags::linear;
+  static constexpr auto tangentModuliTag       = TangentModuliTags::Material;
+  static constexpr bool energyAcceptsVoigt     = Base::energyAcceptsVoigt;
+  static constexpr bool stressToVoigt          = Base::stressToVoigt;
+  static constexpr bool stressAcceptsVoigt     = Base::stressAcceptsVoigt;
+  static constexpr bool moduliToVoigt          = Base::moduliToVoigt;
+  static constexpr bool moduliAcceptsVoigt     = Base::moduliAcceptsVoigt;
+  static constexpr double derivativeFactorImpl = Base::derivativeFactorImpl;
+
+  [[nodiscard]] constexpr static std::string nameImpl() noexcept { return "LinearElasticity"; }
+
   /**
    * \brief Constructor for LinearElasticityT.
    *
    * \param mpt The LamesFirstParameterAndShearModulus object.
    */
-  explicit LinearElasticityT(const LamesFirstParameterAndShearModulus& mpt)
+  explicit LinearElasticityT(const MaterialParameters& mpt)
       : svk_{mpt} {}
 
-  using field_type                    = ScalarType;
-  static constexpr int worldDimension = 3;
-  using StrainMatrix                  = Eigen::Matrix<ScalarType, worldDimension, worldDimension>;
-  using StressMatrix                  = StrainMatrix;
-
-  static constexpr auto strainTag          = StrainTags::linear;
-  static constexpr auto stressTag          = StressTags::linear;
-  static constexpr auto tangentModuliTag   = TangentModuliTags::Material;
-  static constexpr bool energyAcceptsVoigt = Base::energyAcceptsVoigt;
-  static constexpr bool stressToVoigt      = Base::stressToVoigt;
-  static constexpr bool stressAcceptsVoigt = Base::stressAcceptsVoigt;
-  static constexpr bool moduliToVoigt      = Base::moduliToVoigt;
-  static constexpr bool moduliAcceptsVoigt = Base::moduliAcceptsVoigt;
-
-  // This factor denotes the differences between the returned stresses and moduli
-  // and the passed strain. For neoHooke the inserted quantity is C, the right Cauchy-Green tensor.
-  // The function relation between the energy and the stresses is S = 2*\partial \psi(C)/ \partial C.
-  // This factor is a pre-factor, which is the difference to the actual derivative is written here.
-  static constexpr double derivativeFactor = 1;
+  /**
+   * \brief Returns the material parameters stored in the material
+   */
+  MaterialParameters materialParametersImpl() const { return svk_.materialParametersImpl(); }
 
   /**
    * \brief Calculates the stored energy in the material.

ikarus/finiteelements/mechanics/materials/neohooke.hh

@@ -35,29 +35,35 @@ namespace Ikarus {
 template <typename ST>
 struct NeoHookeT : public Material<NeoHookeT<ST>>
 {
-  [[nodiscard]] constexpr std::string nameImpl() const noexcept { return "NeoHooke"; }
+  using ScalarType                    = ST;
+  static constexpr int worldDimension = 3;
+  using StrainMatrix                  = Eigen::Matrix<ScalarType, worldDimension, worldDimension>;
+  using StressMatrix                  = StrainMatrix;
+  using MaterialParameters            = LamesFirstParameterAndShearModulus;
+
+  static constexpr auto strainTag              = StrainTags::rightCauchyGreenTensor;
+  static constexpr auto stressTag              = StressTags::PK2;
+  static constexpr auto tangentModuliTag       = TangentModuliTags::Material;
+  static constexpr bool energyAcceptsVoigt     = false;
+  static constexpr bool stressToVoigt          = false;
+  static constexpr bool stressAcceptsVoigt     = false;
+  static constexpr bool moduliToVoigt          = false;
+  static constexpr bool moduliAcceptsVoigt     = false;
+  static constexpr double derivativeFactorImpl = 2;
+
+  [[nodiscard]] constexpr static std::string nameImpl() noexcept { return "NeoHooke"; }
 
   /**
    * \brief Constructor for NeoHookeT.
    * \param mpt The Lame's parameters (first parameter and shear modulus).
    */
-  explicit NeoHookeT(const LamesFirstParameterAndShearModulus& mpt)
-      : lambdaAndmu_{mpt} {}
+  explicit NeoHookeT(const MaterialParameters& mpt)
+      : materialParameter_{mpt} {}
 
-  using ScalarType                    = ST;
-  static constexpr int worldDimension = 3;
-  using StrainMatrix                  = Eigen::Matrix<ScalarType, worldDimension, worldDimension>;
-  using StressMatrix                  = StrainMatrix;
-
-  static constexpr auto strainTag          = StrainTags::rightCauchyGreenTensor;
-  static constexpr auto stressTag          = StressTags::PK2;
-  static constexpr auto tangentModuliTag   = TangentModuliTags::Material;
-  static constexpr bool energyAcceptsVoigt = false;
-  static constexpr bool stressToVoigt      = false;
-  static constexpr bool stressAcceptsVoigt = false;
-  static constexpr bool moduliToVoigt      = false;
-  static constexpr bool moduliAcceptsVoigt = false;
-  static constexpr double derivativeFactor = 2;
+  /**
+   * \brief Returns the material parameters stored in the material
+   */
+  MaterialParameters materialParametersImpl() const { return materialParameter_; }
 
   /**
    * \brief Computes the stored energy in the Neo-Hookean material model.
@@ -71,7 +77,8 @@ struct NeoHookeT : public Material<NeoHookeT<ST>>
     if constexpr (!Concepts::EigenVector<Derived>) {
       const auto traceC  = C.trace();
       const auto logdetF = log(sqrt(C.determinant()));
-      return lambdaAndmu_.mu / 2.0 * (traceC - 3 - 2 * logdetF) + lambdaAndmu_.lambda / 2.0 * logdetF * logdetF;
+      return materialParameter_.mu / 2.0 * (traceC - 3 - 2 * logdetF) +
+             materialParameter_.lambda / 2.0 * logdetF * logdetF;
     } else
       static_assert(!Concepts::EigenVector<Derived>,
                     "NeoHooke energy can only be called with a matrix and not a vector in Voigt notation");
@@ -91,7 +98,8 @@ struct NeoHookeT : public Material<NeoHookeT<ST>>
       if constexpr (!Concepts::EigenVector<Derived>) {
         const auto logdetF = log(sqrt(C.determinant()));
         const auto invC    = C.inverse().eval();
-        return (lambdaAndmu_.mu * (StrainMatrix::Identity() - invC) + lambdaAndmu_.lambda * logdetF * invC).eval();
+        return (materialParameter_.mu * (StrainMatrix::Identity() - invC) + materialParameter_.lambda * logdetF * invC)
+            .eval();
       } else
         static_assert(!Concepts::EigenVector<Derived>,
                       "NeoHooke can only be called with a matrix and not a vector in Voigt notation");
@@ -115,8 +123,9 @@ struct NeoHookeT : public Material<NeoHookeT<ST>>
       const auto CTinv   = tensorView(invC, std::array<Eigen::Index, 2>({3, 3}));
       static_assert(Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3>>::NumIndices == 2);
       Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>> moduli =
-          (lambdaAndmu_.lambda * dyadic(CTinv, CTinv) +
-           2 * (lambdaAndmu_.mu - lambdaAndmu_.lambda * logdetF) * symTwoSlots(fourthOrderIKJL(invC, invC), {2, 3}))
+          (materialParameter_.lambda * dyadic(CTinv, CTinv) +
+           2 * (materialParameter_.mu - materialParameter_.lambda * logdetF) *
+               symTwoSlots(fourthOrderIKJL(invC, invC), {2, 3}))
               .eval();
       return moduli;
     } else
@@ -130,11 +139,11 @@ struct NeoHookeT : public Material<NeoHookeT<ST>>
    */
   template <typename STO>
   auto rebind() const {
-    return NeoHookeT<STO>(lambdaAndmu_);
+    return NeoHookeT<STO>(materialParameter_);
   }
 
 private:
-  LamesFirstParameterAndShearModulus lambdaAndmu_;
+  MaterialParameters materialParameter_;
 };
 
 /**

ikarus/finiteelements/mechanics/materials/strainconversions.hh

@@ -83,7 +83,7 @@ decltype(auto) createDeformationGradient(const Eigen::MatrixBase<Derived>& eMB)
   else if constexpr (tag == StrainTags::displacementGradient)
     return (e + Derived::Identity()).eval();
   else if constexpr (tag == StrainTags::rightCauchyGreenTensor) {
-    return (e.sqrt()).eval(); // this looses information, since the rotation information from an original F is lost!
+    return (e.sqrt()).eval(); // this looses information, since the rotation information from the original F is lost!
   }
 }
 
@@ -125,14 +125,14 @@ decltype(auto) createRightCauchyGreen(const Eigen::MatrixBase<Derived>& eMB) {
  * \tparam from Type of the source strain tag.
  * \tparam to Type of the target strain tag.
  * \tparam Derived Type of the Eigen matrix.
- * \param eRaw Eigen matrix representing the input strain.
+ * \param eRaw Eigen matrix representing the input strain (can be in Voigt notation).
  * \return The transformed strain matrix.
  */
 template <StrainTags from, StrainTags to, typename Derived>
 decltype(auto) transformStrain(const Eigen::MatrixBase<Derived>& eRaw) {
   static_assert((from == to) or (from != StrainTags::linear and to != StrainTags::linear),
                 "No useful transformation available for linear strains.");
-  decltype(auto) e = Impl::mayTransformFromVoigt(eRaw, true);
+  decltype(auto) e = Impl::mayTransformFromVoigt(eRaw.eval(), true);
   if constexpr (from == to)
     return e;
   else if constexpr (to == StrainTags::greenLagrangian)