Closed thelfer closed 4 years ago
I first want to have comments about my first implementation of the CalculateMaterialResponsePK2
.
First I created an Integrate
method which calls the constitutive law integration. This method will be discussed later. After the call to this method, I assume that I have computed the stress and the consistent tangent operator.
It is worh discussing that for finite strain behaviour, my current choice is to always compute the Cauchy stress and use the derivative of the second Piola-Kirchoff stress with respect the Green-Lagrange strain as the consistent tangent operator (the generic
interface allows various stress measure on input/output and various tangent operator on ouput, but this must be choosen when loading the behaviour, so I made this choice by default).
Once the behaviour integration succeeded, I need to provide the stress and the consistent tangent operator to Kratos.
I would like to have feed-backs on the implemntation. I also would like to have some documentation about the functions already implemented in Kratos
for the conversion between the stress measures and tangent operators.
Here is the implementation:
void MGISConstitutiveLaw::CalculateMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) {
KRATOS_TRY;
this->Integrate(rValues);
if (opts.Is(ConstitutiveLaw::COMPUTE_STRESS)) {
// exporting the stress
if (this->strain_measure == StrainMeasure_Deformation_Gradient) {
auto& pk2 = rValues.GetStressVector();
auto s = pk2; // Cauchy stress
Kratos::MGIS::ConvertStressToKratos(s, this->data.s1.thermodynamic_forces.data());
// HERE, use Kratos to convert s to pk2...
} else if (this->strain_measure == StrainMeasure_Infinitesimal) {
auto& s = rValues.GetStressVector();
Kratos::MGIS::ConvertStressToKratos(s, this->data.s1.thermodynamic_forces.data());
} else {
KRATOS_ERROR << "unsupported behaviour type";
}
}
if (opts.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
// the output of the behaviour is the derivative of the second Piola-Kirchhoff stres with
// respect to the Green-Lagrange strain (see the `MGISConstitutiveLawFactory` for the finite
// strain behaviour options).
Matrix K;
Kratos::MGIS::ConvertTangentOperatorToKratos(K, this->data.K.data(),
this->behaviour->hypothesis);
rValues.SetConstitutiveMatrix(K);
}
KRATOS_CATCH("");
}
I first want to have comments about my first implementation of the
CalculateMaterialResponsePK2
.First I created an
Integrate
method which calls the constitutive law integration. This method will be discussed later. After the call to this method, I assume that I have computed the stress and the consistent tangent operator.It is worh discussing that for finite strain behaviour, my current choice is to always compute the Cauchy stress and use the derivative of the second Piola-Kirchoff stress with respect the Green-Lagrange strain as the consistent tangent operator (the
generic
interface allows various stress measure on input/output and various tangent operator on ouput, but this must be choosen when loading the behaviour, so I made this choice by default).Once the behaviour integration succeeded, I need to provide the stress and the consistent tangent operator to Kratos.
I would like to have feed-backs on the implemntation. I also would like to have some documentation about the functions already implemented in
Kratos
for the conversion between the stress measures and tangent operators.Here is the implementation:
void MGISConstitutiveLaw::CalculateMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) { KRATOS_TRY; this->Integrate(rValues); if (opts.Is(ConstitutiveLaw::COMPUTE_STRESS)) { // exporting the stress if (this->strain_measure == StrainMeasure_Deformation_Gradient) { auto& pk2 = rValues.GetStressVector(); auto s = pk2; // Cauchy stress Kratos::MGIS::ConvertStressToKratos(s, this->data.s1.thermodynamic_forces.data()); // HERE, use Kratos to convert s to pk2... } else if (this->strain_measure == StrainMeasure_Infinitesimal) { auto& s = rValues.GetStressVector(); Kratos::MGIS::ConvertStressToKratos(s, this->data.s1.thermodynamic_forces.data()); } else { KRATOS_ERROR << "unsupported behaviour type"; } } if (opts.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) { // the output of the behaviour is the derivative of the second Piola-Kirchhoff stres with // respect to the Green-Lagrange strain (see the `MGISConstitutiveLawFactory` for the finite // strain behaviour options). Matrix K; Kratos::MGIS::ConvertTangentOperatorToKratos(K, this->data.K.data(), this->behaviour->hypothesis); rValues.SetConstitutiveMatrix(K); } KRATOS_CATCH(""); }
In general looks right, but from my own experience be careful when using the methods SetConstitutiveMatrix and SetStressVector. These methods work with raw pointers, see https://github.com/KratosMultiphysics/Kratos/blob/029cac18c3eb3cefcf7481e8db702d143ee650d5/kratos/includes/constitutive_law.h#L229, and therefore if you assign something temporal, after you exit the CalculateMaterialResponsePK2 the value will disappear from memory. These methods are designed to be used in the element when the matrix is used, and later the recommended is to work directly with the references.
See for example here like we create backups: https://github.com/KratosMultiphysics/Kratos/blob/029cac18c3eb3cefcf7481e8db702d143ee650d5/applications/StructuralMechanicsApplication/custom_advanced_constitutive/generic_finite_strain_isotropic_plasticity.cpp#L131
I will extend my comments later, right now I cannot extend more.
@loumalouomega Thanks for your comments. Indeed, I have done this inspired by what is done in the ViscousGeneralizedKelvin
class:
Matrix constitutive_matrix, inverse_constitutive_matrix;
this->CalculateValue(rValues, CONSTITUTIVE_MATRIX, constitutive_matrix);
...
if (r_flags.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
rValues.SetConstitutiveMatrix(constitutive_matrix);
}
I must have missed something. Would something like the following be better ?
auto& K = rValues.GetConstitutiveMatrix(K);
Kratos::MGIS::ConvertTangentOperatorToKratos(K, this->data.K.data(),
this->behaviour->hypothesis);
@loumalouomega Thanks for your comments. Indeed, I have done this inspired by what is done in the
ViscousGeneralizedKelvin
class:Matrix constitutive_matrix, inverse_constitutive_matrix; this->CalculateValue(rValues, CONSTITUTIVE_MATRIX, constitutive_matrix); ... if (r_flags.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) { rValues.SetConstitutiveMatrix(constitutive_matrix); }
I must have missed something. Would something like the following be better ?
auto& K = rValues.GetConstitutiveMatrix(K); Kratos::MGIS::ConvertTangentOperatorToKratos(K, this->data.K.data(), this->behaviour->hypothesis);
Indeed IMO. @AlejandroCornejo can you check the ViscousGeneralizedKelvin doesn't suffer from memory problems?
BTW, @thelfer, @AlejandroCornejo is our constitutive models expert
Nice to meet you @AlejandroCornejo !
Hello @thelfer :) I'll think about this ASAP
Regarding PushForward/PullBack operations, take a look at the base class "ConstitutiveLaw" and look for "Transform" you will find some useful implementations, although admittedly not all the cases are coded.
also take a look at
PullBackConstitutiveMatrix
I still miss an important information that I did not find in the ConstitutiveLaw
class documentation nor on the wiki ( https://github.com/KratosMultiphysics/Kratos/wiki/How-to-use-the-Constitutive-Law-class): what is the tangent operator expected in each case. For example, I assumed that I shall compute the derivative of the second Piola-Kirchoff stress with respect the Green-Lagrange strain as the consistent tangent operator when the CalculateMaterialResponsePK2
is called. But I don't know what is expected in the other cases. Any help ?
we define C as Dsigma/Deps
where sigma and eps are the stress and strain measures you are assuming.
other than that it shoukd be possible to transform from one to the other via the functions of the base class.
to be honest however i am not too sure of how tested this is for cases other than small strain and pk2 vs cauchy green
Riccardo
On Thu, Nov 21, 2019, 5:53 PM thelfer notifications@github.com wrote:
I still miss an important information that I did not find in the ConstitutiveLaw class documentation nor on the wiki ( https://github.com/KratosMultiphysics/Kratos/wiki/How-to-use-the-Constitutive-Law-class): what is the tangent operator expected in each case. For example, I assumed that I shall compute the derivative of the second Piola-Kirchoff stress with respect the Green-Lagrange strain as the consistent tangent operator when the CalculateMaterialResponsePK2 is called. But I don't know what is expected in the other cases. Any help ?
— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/KratosMultiphysics/Kratos/issues/5957?email_source=notifications&email_token=AB5PWEP6GO3PZCN5T7IOI6DQU24J3A5CNFSM4JP75Z7KYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEE25DJI#issuecomment-557175205, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB5PWEPGIZGOPXGXUKCYSBLQU24J3ANCNFSM4JP75Z7A .
@RiccardoRossi Just to be sure, what is eps for PK1, Kirchhoff or Cauchy ? For example, the conjugate to PK1 is the deformation gadient.. For the Kirchhoff and Cauchy, one usually uses either the spatial modulus or some modulus associated with an objective derivative. But the choice really depends on the solver.
wow, i am pretty sure you will soon catch me on shacky ground, i am definitely not an expert.
we pass F=du/dX (capital x) so that one can compute what is needed and do the necessary transformations internally to the CL.
to the best of my knowledge, all of the elements in the structural mechanics application use the binomium of PK2 stress and GreenLagrange strain tensor E:=1/2(F^t-F-I)
this is so for both our total lagrangian element and for the update lagrangian one (this is a bit strange but we preferred to do it this way)
in small deformation we use eps=symm_gradient(displacement) and the cauchy stress (well ... they all coincide in this context)
hope this answers you better
[Edit:] removing some leftovers
about PK1 i think we are currently not using it anywhere.
we made an attempt to design the CL sufficiently flexible that one could use anything...maybe we made it too complicated in the process
OK. So I will do something simple for the first tests and only implement the PK2 option. Thanks for your help !
With the previous elements, the implementation of CalculateMaterialResponsePK2
is much simpler.
Let us now have a look at the Integrate
method:
/**
* @brief integrate the behaviour over the time step.
* @param rValues The internal values of the law
* @see Parameters
* @note This method does not update the stress and
* the tangent operator on the Kratos side. This must
* be done in the calling methods
* (`CalculateMaterialResponseCauchy` for example)
* which are responsible for converting the stress and
* the tangent operator returned by `MGIS` (which is
* respectively the second Piola-Kirchhoff stress and
* its derivative with respect to the Green-Lagrange
* strain.
*/
void Integrate(ConstitutiveLaw::Parameters& rValues);
and
void MGISConstitutiveLaw::Integrate(ConstitutiveLaw::Parameters& rValues) {
const auto& opts = rValues.GetOptions();
const auto& pi = rValues.GetProcessInfo();
// getting the gradients
if (this->strain_measure == StrainMeasure_Infinitesimal) {
if (opts.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
KRATOS_ERROR << "MGISConstitutiveLaw::Integrate: the strain tensor must be provided by the "
"element\n";
}
auto& ek = rValues.GetStrainVector();
Kratos::MGIS::ConvertStrainToMGIS(this->data.s1.gradients.data(), ek);
} else {
KRATOS_ERROR << "MGISConstitutiveLaw::Integrate: unimplemented yet\n";
}
const auto k = opts.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR);
if (opts.Is(ConstitutiveLaw::COMPUTE_STRESS)) {
this->data.K[0] = k ? 4 : 0;
} else {
if (!k) {
KRATOS_ERROR << "MGISConstitutiveLaw::Integrate: don't know what to do, neither the stress "
"nor the tangent operator is requested\n";
}
// return the elastic operator
this->data.K[0] = -1;
}
data.dt = pi[DELTA_TIME];
auto data_view = mgis::behaviour::make_view(this->data);
if (mgis::behaviour::integrate(data_view, *(this->behaviour)) != 0) {
KRATOS_ERROR << "MGISConstitutiveLaw::Integrate: behaviour integration failed\n";
}
} // end of MGISConstitutiveLaw::Integrate
My questions are:
* Is the strain correctly retrieved ?
Assuming that ConvertStrainToMGIS respect the Voigt notation considered in Kratos, yes
- Is the time increment correctly retrieved ? Yes, it is just DELTA_TIME, nothing else.
- Is the deformation gradient always stored as a 3x3 matrix ? 3x3 in 3D elements and 2x2 in 2D elements. Asymmetric is 3x3 also. If you asking if we consider the Voigt notation for symmetric cases, no, we don't
the point of USE_ELEMENT_PROVIDED_STRAIN is to allow the element OR the CL to compute the strain. I appreciate this is strange, however consider that we always give F which is needed for all the transforms you might need. this implies that instead of writing
if (this->strain_measure == StrainMeasure_Infinitesimal) {
if (opts.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
KRATOS_ERROR << "MGISConstitutiveLaw::Integrate: the strain tensor must be provided by the "
"element\n";
}
auto& ek = rValues.GetStrainVector();
Kratos::MGIS::ConvertStrainToMGIS(this->data.s1.gradients.data(), ek);
}
you could do something like
if (this->strain_measure == StrainMeasure_Infinitesimal) {
auto& ek = rValues.GetStrainVector();
if (opts.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
Matrix E = 1/2*(prod(trans(F),F) - IdentityMatrix(3));
ek = MathUtils<double>::StrainTensorToVector(E);
}
Kratos::MGIS::ConvertStrainToMGIS(this->data.s1.gradients.data(), ek);
}
note incidentally that this also implies that if instead of E you need internally another strain measure you can compute it out of F
@RiccardoRossi Here you use the Green-Lagrange strain by default. Does that mean that the computations are always done in finite strain (taking into account non-linear geometrical terms) ? My question is motivated by my experience with another solver which decided to use GL strain by default even in the small strain hypothesis and this experience is pretty bad, since It only leads to discrepancies between the results of the code an analytical solutions....
Two options here:
What do you think ?
Also, I find strange that the line
Matrix E = 1/2*(prod(trans(F),F) - IdentityMatrix(3));
ek = MathUtils<double>::StrainTensorToVector(E);
does not take the hypothesis into account...
sorry, you are right about this, i missed the "StrainMeasure_Infinitesimal"
what we normally do in KRatos is that for small strain elements we do provide the strain, while for large deformations we don't.
in any case the code should have been something like (withing it out of my memory so it may be wrong)
if (this->strain_measure == StrainMeasure_Infinitesimal) {
auto& ek = rValues.GetStrainVector();
if (opts.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
ek[0] = F(0,0)-1 //dux/dX0
ek[1] = F(1,1)-1 //duy/dY0
ek[2] = (F(0,1)-1) + (F(1,0)-1) //duy/dX0 + dux/dY0
//this is in 2D of course...
ek = MathUtils<double>::StrainTensorToVector(E);
}
Kratos::MGIS::ConvertStrainToMGIS(this->data.s1.gradients.data(), ek);
}
elseif (this->strain_measure == StrainMeasure_GreenLagrange) {
....
}
sorry, you are right about this, i missed the "StrainMeasure_Infinitesimal"
what we normally do in KRatos is that for small strain elements we do provide the strain, while for large deformations we don't.
in any case the code should have been something like (withing it out of my memory so it may be wrong)
if (this->strain_measure == StrainMeasure_Infinitesimal) { auto& ek = rValues.GetStrainVector(); if (opts.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) { ek[0] = F(0,0)-1 //dux/dX0 ek[1] = F(1,1)-1 //duy/dY0 ek[2] = (F(0,1)-1) + (F(1,0)-1) //duy/dX0 + dux/dY0 //this is in 2D of course... ek = MathUtils<double>::StrainTensorToVector(E); } Kratos::MGIS::ConvertStrainToMGIS(this->data.s1.gradients.data(), ek); } elseif (this->strain_measure == StrainMeasure_GreenLagrange) { .... }
Just to point we have a CL utilities with many of these operations included:
@RiccardoRossi I think that there is a misunderstanding here.
I think that Kratos
' design is meant to let the CL handle the kinematic (i.e. choosing a strain measure, or use a corotationnal approach, or use any objective stress rate). This design is a good thing in general.
However this is redundant with how MFront
works.
I'll try to make things clear in the following. Let me know if you agree with me.
In MGIS
(i.e. on the solver side), any mechanical constitutive law that expects a strain is classified as STANDARDSTRAINBASEDBEHAVIOUR
and any mechanical consitutive law that expects that a deformation gradient is called STANDARDFINITESTRAINBEHAVIOUR
.
The point is that you can write a MFront
behaviour based on a specific strain measure (for example the Hencky strain or the Green-Lagrange strain). In this case, MFront
will
Such a behaviour will appear to the calling solver as a STANDARDFINITESTRAINBEHAVIOUR
.
MFront
behaviours in Kratos
So, putting it all together, it seems to me that the only consistent choice (again in the MFront
case) is the following:
STANDARDSTRAINBASEDBEHAVIOUR
type it must be used in Kratos
under the small strain assumptions and the element must provide the linearised strain (i.e. no computation of any strain measure in the CL).STANDARDFINITESTRAINBEHAVIOUR
.type, it must be used in Kratos
in finite strain analys. But I can't guarantee that the CL's is correctly used because all the elements provides the deformation gradient, unless the element queries the law features and forbids any misuse (see below the current implementation of the GetLawFeatures
).Currently, I assigned this->strain_measure
to Constitutive::StrainMeasure_Infinitesimal
for all behaviours that expects a strain, so the MGISConstitutiveLaw
constructor:
if ((this->behaviour->btype == Behaviour::STANDARDSTRAINBASEDBEHAVIOUR) &&
(this->behaviour->kinematic == Behaviour::SMALLSTRAINKINEMATIC)) {
this->strain_measure = StrainMeasure_Infinitesimal;
this->stress_measure = StressMeasure_Cauchy;
} else if ((this->behaviour->btype == Behaviour::STANDARDFINITESTRAINBEHAVIOUR) &&
(this->behaviour->kinematic == Behaviour::FINITESTRAINKINEMATIC_F_CAUCHY)) {
this->strain_measure = StrainMeasure_Deformation_Gradient;
this->stress_measure = StressMeasure_Cauchy;
} else {
(note that the behaviour kinematic is not meaningfull for the generic
interface)
I also do this in the GetLawFeatures
method:
if (this->strain_measure == StrainMeasure_Infinitesimal) {
rFeatures.mOptions.Set(INFINITESIMAL_STRAINS);
} else if (this->strain_measure == StrainMeasure_Deformation_Gradient){
rFeatures.mOptions.Set(FINITE_STRAINS);
} else {
KRATOS_ERROR << "unsupported behaviour type\n";
}
Oh i see. So our design is not too different!
i think that the solution you propose makes sense!
Oh i see. So our design is not too different!
i think that the solution you propose makes sense!
:)
No, they ain't ! Indeed, in most solvers, the kinematic is handled by the element and not the CL, which is bad IMHO. I got confused, because I am not used to see solvers making things right :)
Looks good now (just pulled a few changes). In can now integrate a simple elastic behaviour and results are ok !
from __future__ import print_function, absolute_import, division
import KratosMultiphysics
import KratosMultiphysics.StructuralMechanicsApplication as StructuralMechanicsApplication
import KratosMultiphysics.MGISApplication
import KratosMultiphysics.KratosUnittest as KratosUnittest
import Kratos
mgis_factory = KratosMultiphysics.MGISApplication.MGISConstitutiveLawFactory()
p = Kratos.Parameters()
l = Kratos.Parameters()
l.SetString('src/libBehaviour.so')
p.AddValue('library', l)
b = Kratos.Parameters()
b.SetString('Elasticity')
p.AddValue('behaviour', b)
h = Kratos.Parameters()
h.SetString('Tridimensional')
p.AddValue('hypothesis', h)
b = mgis_factory.Create(p)
current_model = KratosMultiphysics.Model()
model_part = current_model.CreateModelPart("test")
cl_options = KratosMultiphysics.Flags()
cl_options.Set(KratosMultiphysics.ConstitutiveLaw.COMPUTE_STRESS, True)
cl_options.Set(KratosMultiphysics.ConstitutiveLaw.COMPUTE_CONSTITUTIVE_TENSOR, True)
# Define deformation gradient
F = KratosMultiphysics.Matrix(6,6)
for i in range(6):
for j in range(6):
if(i==j):
F[i,j] = 1.0
else:
F[i,j] = 0.0
detF = 1.0
stress_vector = KratosMultiphysics.Vector(b.GetStrainSize())
strain_vector = KratosMultiphysics.Vector(b.GetStrainSize())
constitutive_matrix = KratosMultiphysics.Matrix(b.GetStrainSize(),b.GetStrainSize())
for i in range(0, b.GetStrainSize()):
stress_vector[i] = 0.0
strain_vector[i] = 0.0
for j in range(0, b.GetStrainSize()):
constitutive_matrix[i,j] = 0.0
strain_vector[0] = 1e-3
cl_params = KratosMultiphysics.ConstitutiveLawParameters()
cl_params.SetOptions(cl_options)
cl_params.SetDeformationGradientF(F)
cl_params.SetDeterminantF(detF)
cl_params.SetStrainVector(strain_vector)
cl_params.SetStressVector(stress_vector)
cl_params.SetConstitutiveMatrix(constitutive_matrix)
cl_params.SetProcessInfo(model_part.ProcessInfo)
b.CalculateMaterialResponseCauchy(cl_params)
print(cl_params.GetStressVector())
print(cl_params.GetConstitutiveMatrix())
BTW, is there a way to make this block more pleasant to read:
mgis_factory = KratosMultiphysics.MGISApplication.MGISConstitutiveLawFactory()
p = Kratos.Parameters()
l = Kratos.Parameters()
l.SetString('src/libBehaviour.so')
p.AddValue('library', l)
b = Kratos.Parameters()
b.SetString('Elasticity')
p.AddValue('behaviour', b)
h = Kratos.Parameters()
h.SetString('Tridimensional')
p.AddValue('hypothesis', h)
b = mgis_factory.Create(p)
@thelfer FYI, i am adding a variable properties with Parameters in https://github.com/KratosMultiphysics/Kratos/tree/core/add-CONFIGURATION_PARAMETERS. I am having problems with pybind, I assume will be fixed soon :P
@loumalouomega don't know what it means, but looks great.
@loumalouomega don't know what it means, but looks great.
You will able to store Parameters in Properties, to have fully access in the CL without further modification
Close this issue for the moment, as things seems working correctly in the integration part now.
This issue is meant to discuss the constitutive law integration.