Skip to content

Commit

Permalink
Merge branch 'MFrontUpdatesForLargeDeformations' into 'master'
Browse files Browse the repository at this point in the history
MFront update in preparation for LargeDeformation process

See merge request ogs/ogs!4794
  • Loading branch information
endJunction committed Nov 27, 2023
2 parents db89e1a + 27af6d3 commit 74b394f
Show file tree
Hide file tree
Showing 8 changed files with 255 additions and 28 deletions.
54 changes: 54 additions & 0 deletions MaterialLib/MPL/Utils/Tensor.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
/*
* \file
* \copyright
* Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/

#pragma once

#include <Eigen/Core>

#include "MaterialLib/MPL/Property.h"

namespace MaterialPropertyLib
{
/// See Tensor type for details.
constexpr int tensorSize(int dim)
{
if (dim == 1)
{
return 3; // Diagonal entries.
}
if (dim == 2)
{
return 5; // 2x2 matrix and the 3rd diagonal entry.
}
if (dim == 3)
{
return 9; // Full 3x3 matrix.
}
OGS_FATAL("Tensor size for dimension {} is not defined.", dim);
}

/// The tensor's components in 3D case are ordered in the usual way:
/// (1,1), (1,2), (1,3)
/// (2,1), (2,2), (2,3)
/// (3,1), (3,2), (3,3).
///
/// For the 2D case the 2x2 block is as usual and is followed by the (3,3)
/// component:
/// (1,1), (1,2),
/// (2,1), (2,2),
/// (3,3).
///
/// For the 1D case only the diagonal is stored:
/// (1,1),
/// (2,2),
/// (3,3).
template <int Dim>
using Tensor = Eigen::Matrix<double, tensorSize(Dim), 1>;

} // namespace MaterialPropertyLib
14 changes: 13 additions & 1 deletion MaterialLib/MPL/VariableType.h
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ enum class Variable : int
{
capillary_pressure,
concentration,
deformation_gradient,
density,
effective_pore_pressure,
enthalpy,
Expand Down Expand Up @@ -58,6 +59,7 @@ static const std::array<std::string,
static_cast<int>(Variable::number_of_variables)>
variable_enum_to_string{{"capillary_pressure",
"concentration",
"deformation_gradient",
"density",
"effective_pore_pressure",
"enthalpy",
Expand Down Expand Up @@ -86,7 +88,9 @@ static const std::array<std::string,
using VariableType = std::variant<std::monostate,
double,
Eigen::Matrix<double, 4, 1>,
Eigen::Matrix<double, 6, 1>>;
Eigen::Matrix<double, 5, 1>,
Eigen::Matrix<double, 6, 1>,
Eigen::Matrix<double, 9, 1>>;

class VariableArray
{
Expand All @@ -102,6 +106,8 @@ class VariableArray
return capillary_pressure;
case Variable::concentration:
return concentration;
case Variable::deformation_gradient:
return std::visit(identity, deformation_gradient);
case Variable::density:
return density;
case Variable::effective_pore_pressure:
Expand Down Expand Up @@ -156,6 +162,12 @@ class VariableArray

double capillary_pressure = nan_;
double concentration = nan_;
// Compare to GMatrixPolicy::GradientVectorType. The 1d case = Matrix<3, 1>
// is not used so far.
std::variant<std::monostate,
Eigen::Matrix<double, 5, 1>,
Eigen::Matrix<double, 9, 1>>
deformation_gradient;
double density = nan_;
double effective_pore_pressure = nan_;
double enthalpy = nan_;
Expand Down
61 changes: 61 additions & 0 deletions MaterialLib/SolidModels/MFront/MFrontGeneric.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include <boost/mp11.hpp>

#include "MaterialLib/MPL/Utils/GetSymmetricTensor.h"
#include "MaterialLib/MPL/Utils/Tensor.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "NumLib/Exceptions.h"
#include "ParameterLib/Parameter.h"
Expand Down Expand Up @@ -47,6 +48,45 @@ constexpr auto eigenSwap45View(Eigen::MatrixBase<Derived> const& matrix)
});
}

/// Converts between OGS' and MFront's tensors, which are represented as
/// vectors. An OGS tensor
/// 11 12 13 21 22 23 31 32 33
/// 0 1 2 3 4 5 6 7 8
/// is converted to MFront tensor
/// 11 22 33 12 21 13 31 23 32
/// 0 4 8 1 3 2 6 5 7.
template <int DisplacementDim, typename Derived>
constexpr auto ogsTensorToMFrontTensor(Eigen::MatrixBase<Derived> const& matrix)
{
using Matrix =
Eigen::Matrix<typename Derived::Scalar, Derived::RowsAtCompileTime,
Derived::ColsAtCompileTime>;

if constexpr (DisplacementDim == 2)
{
return Matrix::NullaryExpr(
matrix.rows(), matrix.cols(),
[&m = matrix.derived()](Eigen::Index const row,
Eigen::Index const col)
{
constexpr std::ptrdiff_t result[5] = {0, 3, 4, 1, 2};
return m(result[row], result[col]);
});
}
if constexpr (DisplacementDim == 3)
{
return Matrix::NullaryExpr(
matrix.rows(), matrix.cols(),
[&m = matrix.derived()](Eigen::Index const row,
Eigen::Index const col)
{
constexpr std::ptrdiff_t result[9] = {0, 4, 8, 1, 3,
2, 6, 5, 7};
return m(result[row], result[col]);
});
}
}

const char* varTypeToString(int v);

int getEquivalentPlasticStrainOffset(mgis::behaviour::Behaviour const& b);
Expand Down Expand Up @@ -120,6 +160,12 @@ struct MapToMPLType<DisplacementDim, mgis::behaviour::Variable::Type::STENSOR>
using type = MaterialPropertyLib::SymmetricTensor<DisplacementDim>;
};

template <int DisplacementDim>
struct MapToMPLType<DisplacementDim, mgis::behaviour::Variable::Type::TENSOR>
{
using type = MaterialPropertyLib::Tensor<DisplacementDim>;
};

template <int DisplacementDim>
struct MapToMPLType<DisplacementDim, mgis::behaviour::Variable::Type::SCALAR>
{
Expand Down Expand Up @@ -160,6 +206,21 @@ struct SetGradient
: eigenSwap45View(grad_ogs).eval();
std::copy_n(grad_mfront.data(), num_comp, target);
}
else if constexpr (Grad::type ==
mgis::behaviour::Variable::Type::TENSOR)
{
using MPLType = MapToMPLType_t<DisplacementDim, Grad::type>;
auto const& grad_ogs =
std::get<MPLType>(variable_array.*Grad::mpl_var);

if (Q.has_value())
{
OGS_FATAL("Rotations of tensors are not implemented.");
}

Eigen::Map<Eigen::Vector<double, num_comp>>{target} =
ogsTensorToMFrontTensor<DisplacementDim>(grad_ogs);
}
else
{
OGS_FATAL("Unsupported gradient type {}.",
Expand Down
3 changes: 2 additions & 1 deletion MaterialLib/SolidModels/MFront/TangentOperatorBlocksView.h
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include <range/v3/view/zip.hpp>

#include "BaseLib/cpp23.h"
#include "MaterialLib/MPL/Utils/Tensor.h"
#include "MathLib/KelvinVector.h"

namespace MaterialLib::Solids::MFront
Expand Down Expand Up @@ -221,7 +222,7 @@ class OGSMFrontTangentOperatorBlocksView
case VT::VECTOR:
return DisplacementDim;
case VT::TENSOR:
return DisplacementDim * DisplacementDim;
return MaterialPropertyLib::tensorSize(DisplacementDim);
}

OGS_FATAL("Unsupported variable type {}", BaseLib::to_underlying(vt));
Expand Down
55 changes: 53 additions & 2 deletions MaterialLib/SolidModels/MFront/Variable.h
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@

#include <MGIS/Behaviour/Variable.hxx>

#include "MaterialLib/MPL/Utils/Tensor.h"
#include "MaterialLib/MPL/VariableType.h"
#include "MathLib/KelvinVector.h"

Expand Down Expand Up @@ -47,7 +48,7 @@ struct Variable
return MathLib::KelvinVector::kelvin_vector_dimensions(
DisplacementDim);
case T::TENSOR:
return DisplacementDim;
return MaterialPropertyLib::tensorSize(DisplacementDim);
}
}

Expand All @@ -65,7 +66,7 @@ struct Variable
case T::STENSOR:
return 1;
case T::TENSOR:
return DisplacementDim;
return 1;
}
}
};
Expand All @@ -90,6 +91,44 @@ struct Strain : Variable<Strain>
/// Instance that can be used for overload resolution/template type deduction.
static constexpr Strain strain;

/// Meta data for Green-Lagrange strain.
struct GreenLagrangeStrain : Variable<GreenLagrangeStrain>
{
/// The name of the variable in MFront.
constexpr static const char* name = "GreenLagrangeStrain";

/// The type of the variable in MFront.
constexpr static mgis::behaviour::Variable::Type type =
mgis::behaviour::Variable::Type::STENSOR;

/// The VariableArray entry that holds this variable in OGS.
///
/// \note Currently we always pass strain via mechanical_strain.
constexpr static auto mpl_var =
&MaterialPropertyLib::VariableArray::mechanical_strain;
};

/// Instance that can be used for overload resolution/template type deduction.
static constexpr GreenLagrangeStrain green_lagrange_strain;

/// Meta data for deformation gradient.
struct DeformationGradient : Variable<DeformationGradient>
{
/// The name of the variable in MFront.
constexpr static const char* name = "DeformationGradient";

/// The type of the variable in MFront.
constexpr static mgis::behaviour::Variable::Type type =
mgis::behaviour::Variable::Type::TENSOR;

/// The VariableArray entry that holds this variable in OGS.
constexpr static auto mpl_var =
&MaterialPropertyLib::VariableArray::deformation_gradient;
};

/// Instance that can be used for overload resolution/template type deduction.
static constexpr DeformationGradient deformation_gradient;

struct LiquidPressure : Variable<LiquidPressure>
{
constexpr static const char* name = "LiquidPressure";
Expand All @@ -115,6 +154,18 @@ struct Stress : Variable<Stress>

static constexpr Stress stress;

struct SecondPiolaKirchhoffStress : Variable<SecondPiolaKirchhoffStress>
{
constexpr static const char* name = "SecondPiolaKirchhoffStress";

constexpr static mgis::behaviour::Variable::Type type =
mgis::behaviour::Variable::Type::STENSOR;

constexpr static auto mpl_var = &MaterialPropertyLib::VariableArray::stress;
};

static constexpr SecondPiolaKirchhoffStress second_piola_kirchhoff_stress;

struct Saturation : Variable<Saturation>
{
constexpr static const char* name = "Saturation";
Expand Down
40 changes: 40 additions & 0 deletions Tests/MaterialLib/OgsToMFrontConversion.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
/**
* \file
* \copyright
* Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/

#ifdef OGS_USE_MFRONT

#include <gmock/gmock-matchers.h>
#include <gtest/gtest.h>

#include "MaterialLib/SolidModels/MFront/MFrontGeneric.h"
#include "Tests/TestTools.h"

namespace MPL = MaterialPropertyLib;
namespace MSM = MaterialLib::Solids::MFront;

TEST(MaterialLib_OgsToMFrontConversion, Tensor3D)
{
MPL::Tensor<3> ogs_tensor;
ogs_tensor << 11, 12, 13, 21, 22, 23, 31, 32, 33;

EXPECT_THAT(MSM::ogsTensorToMFrontTensor<3>(ogs_tensor).eval(),
testing::Pointwise(testing::DoubleEq(),
{11, 22, 33, 12, 21, 13, 31, 23, 32}));
}

TEST(MaterialLib_OgsToMFrontConversion, Tensor2D)
{
MPL::Tensor<2> ogs_tensor;
ogs_tensor << 11, 12, 21, 22, 33;

EXPECT_THAT(MSM::ogsTensorToMFrontTensor<2>(ogs_tensor).eval(),
testing::Pointwise(testing::DoubleEq(), {11, 22, 33, 12, 21}));
}

#endif
Loading

0 comments on commit 74b394f

Please sign in to comment.