diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt index 261e40f102f0a0b43a3b151a739812ca2c445baf..a3a7536d90948486b2e83c2f656446e00cb2c2a9 100644 --- a/MaterialLib/CMakeLists.txt +++ b/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 diff --git a/MaterialLib/SolidModels/KelvinVector-impl.h b/MaterialLib/SolidModels/KelvinVector-impl.h new file mode 100644 index 0000000000000000000000000000000000000000..0474792a9da02cad37b039716d198912b9489ec1 --- /dev/null +++ b/MaterialLib/SolidModels/KelvinVector-impl.h @@ -0,0 +1,106 @@ +/** + * \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 diff --git a/MaterialLib/SolidModels/KelvinVector.cpp b/MaterialLib/SolidModels/KelvinVector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ac872fd2aa8b7e0174ba8637e8bda7bb37e7ea4b --- /dev/null +++ b/MaterialLib/SolidModels/KelvinVector.cpp @@ -0,0 +1,61 @@ +/** + * \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 diff --git a/MaterialLib/SolidModels/KelvinVector.h b/MaterialLib/SolidModels/KelvinVector.h new file mode 100644 index 0000000000000000000000000000000000000000..fc41442d620e6583002bf0f1671d6a27dcb94e30 --- /dev/null +++ b/MaterialLib/SolidModels/KelvinVector.h @@ -0,0 +1,84 @@ +/** + * \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