Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Orthotropy in MFront solid material models #2774

Merged
merged 10 commits into from
Feb 1, 2020
Prev Previous commit
Next Next commit
[MatL] MFront; Add back and forth rotations.
The tensors are rotated before being passed to MFront.
The results of MFront including the stiffness tensor.
  • Loading branch information
endJunction committed Feb 1, 2020
commit ca9780d4f7e5ea428ed6cbb57fdd7fba4c85ed52
44 changes: 39 additions & 5 deletions MaterialLib/SolidModels/MFront/MFront.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -258,6 +258,8 @@ MFront<DisplacementDim>::integrateStress(
material_state_variables,
double const T) const
{
using namespace MathLib::KelvinVector;

assert(
dynamic_cast<MaterialStateVariables const*>(&material_state_variables));
// New state, copy of current one, packed in unique_ptr for return.
Expand Down Expand Up @@ -290,19 +292,44 @@ MFront<DisplacementDim>::integrateStress(

auto v = mgis::behaviour::make_view(behaviour_data);

auto const eps_prev_MFront = OGSToMFront(eps_prev);
// rotation tensor
auto const Q = [this, &x]() -> KelvinMatrixType<DisplacementDim> {
if (!_local_coordinate_system)
{
return KelvinMatrixType<DisplacementDim>::Identity();
}
return fourthOrderRotationMatrix(
_local_coordinate_system->transformation<DisplacementDim>(x));
}();

auto const eps_prev_MFront =
OGSToMFront(Q.transpose()
.template topLeftCorner<
KelvinVectorDimensions<DisplacementDim>::value,
KelvinVectorDimensions<DisplacementDim>::value>() *
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is the topLeftCorner stuff needed now that we have a 2D (4x4) variant of Q?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be brief, yes.

eps_prev);
for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
{
v.s0.gradients[i] = eps_prev_MFront[i];
}

auto const eps_MFront = OGSToMFront(eps);
auto const eps_MFront =
OGSToMFront(Q.transpose()
.template topLeftCorner<
KelvinVectorDimensions<DisplacementDim>::value,
KelvinVectorDimensions<DisplacementDim>::value>() *
eps);
for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
{
v.s1.gradients[i] = eps_MFront[i];
}

auto const sigma_prev_MFront = OGSToMFront(sigma_prev);
auto const sigma_prev_MFront =
OGSToMFront(Q.transpose()
.template topLeftCorner<
KelvinVectorDimensions<DisplacementDim>::value,
KelvinVectorDimensions<DisplacementDim>::value>() *
sigma_prev);
for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
{
v.s0.thermodynamic_forces[i] = sigma_prev_MFront[i];
Expand All @@ -321,7 +348,10 @@ MFront<DisplacementDim>::integrateStress(
{
sigma[i] = behaviour_data.s1.thermodynamic_forces[i];
}
sigma = MFrontToOGS(sigma);
sigma = Q.template topLeftCorner<
KelvinVectorDimensions<DisplacementDim>::value,
KelvinVectorDimensions<DisplacementDim>::value>() *
MFrontToOGS(sigma);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We rotate sigma back into OGS coordinate system, so epsilon and sigma will be output in xyz. All MFront internal variables (ElasticStrain for example) will be output in the local coordinate system. in the _behaviour.symmetry != mgis::behaviour::Behaviour::Symmetry::ISOTROPIC case we could output a warning that his is the case and users have to be careful. Having the local coordinate system as an output in ParaView for use in transformations would be convenient (but not necessary).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A todo at the place where the rotation should happen added.


// TODO row- vs. column-major storage order. This should only matter for
// anisotropic materials.
Expand All @@ -330,7 +360,11 @@ MFront<DisplacementDim>::integrateStress(
OGS_FATAL("Stiffness matrix has wrong size.");

KelvinMatrix C =
MFrontToOGS(Eigen::Map<KelvinMatrix>(behaviour_data.K.data()));
Q * MFrontToOGS(Eigen::Map<KelvinMatrix>(behaviour_data.K.data())) *
Q.transpose()
.template topLeftCorner<
KelvinVectorDimensions<DisplacementDim>::value,
KelvinVectorDimensions<DisplacementDim>::value>();

return std::make_optional(
std::make_tuple<typename MFront<DisplacementDim>::KelvinVector,
Expand Down