[MatL] Kelvin vector invariants and utilities.

MaterialLib/CMakeLists.txt

@@ -3,6 +3,8 @@ GET_SOURCE_FILES(SOURCES)
 GET_SOURCE_FILES(SOURCES_ADSORPTION Adsorption)
 set(SOURCES ${SOURCES} ${SOURCES_ADSORPTION})
 
+APPEND_SOURCE_FILES(SOURCES SolidModels)
+
 add_library(MaterialLib ${SOURCES} )
 target_link_libraries(MaterialLib
     BaseLib

MaterialLib/SolidModels/KelvinVector-impl.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+namespace MaterialLib
+{
+namespace SolidModels
+{
+template <int KelvinVectorSize>
+double Invariants<KelvinVectorSize>::equivalentStress(
+    Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
+{
+    assert(std::abs(trace(deviatoric_v)) <
+           std::numeric_limits<double>::epsilon());
+    return std::sqrt(3 * J2(deviatoric_v));
+}
+
+template <int KelvinVectorSize>
+double Invariants<KelvinVectorSize>::J2(
+    Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
+{
+    assert(std::abs(trace(deviatoric_v)) <
+           std::numeric_limits<double>::epsilon());
+    return 0.5 * deviatoric_v.transpose() * deviatoric_v;
+}
+
+/// Third invariant, equal to determinant of a deviatoric tensor.
+/// \note The input vector must have trace equal zero.
+template <int KelvinVectorSize>
+double Invariants<KelvinVectorSize>::J3(
+    Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v)
+{
+    assert(std::abs(trace(deviatoric_v)) <
+           std::numeric_limits<double>::epsilon());
+    return determinant(deviatoric_v);
+}
+
+/// Trace of the corresponding tensor.
+template <int KelvinVectorSize>
+double Invariants<KelvinVectorSize>::trace(
+    Eigen::Matrix<double, KelvinVectorSize, 1> const& v)
+{
+    return v.template topLeftCorner<3, 1>().sum();
+}
+
+//
+// Initialization of static Invariant variables.
+//
+
+namespace detail
+{
+template <int KelvinVectorSize>
+Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize>
+initDeviatoricProjection()
+{
+    Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> P_dev =
+        Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize>::Identity();
+
+    P_dev.template topLeftCorner<3, 3>() -=
+        1. / 3. * Eigen::Matrix<double, 3, 3>::Ones();
+    return P_dev;
+}
+
+template <int KelvinVectorSize>
+Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize>
+initSphericalProjection()
+{
+    Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> P_sph =
+        Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize>::Zero();
+
+    P_sph.template topLeftCorner<3, 3>().setConstant(1. / 3.);
+    return P_sph;
+}
+
+template <int KelvinVectorSize>
+Eigen::Matrix<double, KelvinVectorSize, 1> initIdentity2()
+{
+    Eigen::Matrix<double, KelvinVectorSize, 1> ivec =
+        Eigen::Matrix<double, KelvinVectorSize, 1>::Zero();
+
+    ivec.template topLeftCorner<3, 1>().setConstant(1.);
+    return ivec;
+}
+}  // namespace detail
+
+template <int KelvinVectorSize>
+const Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize>
+    Invariants<KelvinVectorSize>::deviatoric_projection =
+        detail::initDeviatoricProjection<KelvinVectorSize>();
+
+template <int KelvinVectorSize>
+Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> const
+    Invariants<KelvinVectorSize>::spherical_projection =
+        detail::initSphericalProjection<KelvinVectorSize>();
+
+template <int KelvinVectorSize>
+const Eigen::Matrix<double, KelvinVectorSize, 1> Invariants<
+    KelvinVectorSize>::identity2 = detail::initIdentity2<KelvinVectorSize>();
+
+}  // namespace SolidModels
+}  // namespace MaterialLib

MaterialLib/SolidModels/KelvinVector.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "KelvinVector.h"
+
+namespace MaterialLib
+{
+namespace SolidModels
+{
+template <>
+double Invariants<6>::determinant(Eigen::Matrix<double, 6, 1> const& v)
+{
+    return v(0) * v(1) * v(2) + v(3) * v(4) * v(5) / std::sqrt(2.) -
+           v(3) * v(3) * v(2) / 2. - v(4) * v(4) * v(0) / 2. -
+           v(5) * v(5) * v(1) / 2.;
+}
+
+template <>
+double Invariants<4>::determinant(Eigen::Matrix<double, 4, 1> const& v)
+{
+    return v(2) * (v(0) * v(1) - v(3) * v(3) / 2.);
+}
+
+template <>
+Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> inverse(
+    Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> const& v)
+{
+    assert(Invariants<4>::determinant(v) != 0);
+
+    Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> inv;
+    inv(0) = v(1) * v(2);
+    inv(1) = v(0) * v(2);
+    inv(2) = v(0) * v(1) - v(3) * v(3) / 2.;
+    inv(3) = -v(3) * v(2);
+    return inv / Invariants<4>::determinant(v);
+}
+
+template <>
+Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> inverse(
+    Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> const& v)
+{
+    assert(Invariants<6>::determinant(v) != 0);
+
+    Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> inv;
+    inv(0) = v(1) * v(2) - v(4) * v(4) / 2.;
+    inv(1) = v(0) * v(2) - v(5) * v(5) / 2.;
+    inv(2) = v(0) * v(1) - v(3) * v(3) / 2.;
+    inv(3) = v(4) * v(5) / std::sqrt(2.) - v(3) * v(2);
+    inv(4) = v(3) * v(5) / std::sqrt(2.) - v(4) * v(0);
+    inv(5) = v(4) * v(3) / std::sqrt(2.) - v(1) * v(5);
+    return inv / Invariants<6>::determinant(v);
+}
+
+}  // namespace SolidModels
+}  // namespace MaterialLib

MaterialLib/SolidModels/KelvinVector.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+#ifndef MATERIALLIB_SOLIDMODELS_KELVINVECTOR_H
+#define MATERIALLIB_SOLIDMODELS_KELVINVECTOR_H
+
+#include <Eigen/Dense>
+
+namespace MaterialLib
+{
+/// The invariants and the Kelving mapping are explained in detail in the
+/// article "On Advantages of the Kelvin Mapping in Finite Element
+/// Implementations of Deformation Processes" \cite Nagel2016.
+namespace SolidModels
+{
+/// Invariants used in mechanics, based on Kelvin representation of the vectors
+/// and matrices.
+/// The invariants are computed at process creation time.
+template <int KelvinVectorSize>
+struct Invariants final
+{
+    static_assert(KelvinVectorSize == 4 || KelvinVectorSize == 6,
+                  "KelvinVector invariants for vectors of size different than "
+                  "4 or 6 is not allowed.");
+    /// Kelvin mapping of deviatoric projection tensor. \f$A_{\rm dev} = P_{\rm
+    /// dev}:A\f$ for \f$A\f$ being a second order tensor.
+    static Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> const
+        deviatoric_projection;
+    /// Kelvin mapping of spherical projection tensor. \f$A_{\rm sph} = P_{\rm
+    /// sph}:A\f$ for \f$A\f$ being a second order tensor.
+    static Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> const
+        spherical_projection;
+    /// Kelvin mapping of 2nd order identity tensor.
+    static Eigen::Matrix<double, KelvinVectorSize, 1> const identity2;
+
+    /// Determinant of a matrix in Kelvin vector representation.
+    static double determinant(
+        Eigen::Matrix<double, KelvinVectorSize, 1> const& v);
+
+    /// The von Mises equivalent stress.
+    /// \note The input vector must have trace equal zero.
+    static double equivalentStress(
+        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);
+
+    /// Second invariant of deviatoric tensor.
+    /// \note The input vector must have trace equal zero.
+    static double J2(
+        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);
+
+    /// Third invariant, equal to determinant of a deviatoric tensor.
+    /// \note The input vector must have trace equal zero.
+    static double J3(
+        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);
+
+    /// Trace of the corresponding tensor.
+    static double trace(Eigen::Matrix<double, KelvinVectorSize, 1> const& v);
+};
+
+//
+// Inverses of a Kelvin vector.
+//
+
+/// Inverse of a matrix in Kelvin vector representation.
+/// There are only implementations for the Kelvin vector size 4 and 6.
+template <int KelvinVectorSize>
+Eigen::Matrix<double, KelvinVectorSize, 1, Eigen::ColMajor, KelvinVectorSize, 1>
+inverse(Eigen::Matrix<double,
+                      KelvinVectorSize,
+                      1,
+                      Eigen::ColMajor,
+                      KelvinVectorSize,
+                      1> const& v);
+
+}  // namespace SolidModels
+}  // namespace MaterialLib
+
+#include "KelvinVector-impl.h"
+
+#endif  // MATERIALLIB_SOLIDMODELS_KELVINVECTOR_H