getfem-users
[Top][All Lists]

## Re: [Getfem-users] Questions about user-defined strain energy

 From: Konstantinos Poulios Subject: Re: [Getfem-users] Questions about user-defined strain energy Date: Wed, 15 Jun 2022 10:53:36 +0200

Dear David,

You can also try the contract operator instead of Trace if you want to get a vector directly:

`Contract(A, i, j)` stands for the contraction of tensor A with respect to its ith and jth indices. The first index is numbered 1. For instance, `Contract(A, 1, 2)` is equivalent to `Trace(A)` for a matrix `A`.

BR
Kostas

On Wed, Jun 15, 2022 at 10:45 AM David Danan <daviddanan9021@gmail.com> wrote:
Dear Yves,

No, i didn't think about trying that, thanks for the idea.

I replaced the following _expression_ (which failed)
value=gf.asm_generic(mim=mim, order=0, _expression_="Norm_sqr(Div(clambda*Trace(Green_Lagrange_E)*Id(2)+2*cmu*Green_Lagrange_E))", model=model,region=-1)
where Green_Lagrange_E is a macro defined by
by
newVal=gf.asm_generic(mim=mim, order=0, _expression_="Norm_sqr(Trace(Grad(clambda*Trace(Green_Lagrange_E)*Id(2)+2*cmu*Green_Lagrange_E)))", model=model,region=-1)

which also failed (Trace operator is for square matrices only) so i tried that instead
newVal_x=gf.asm_generic(mim=mim, order=0, _expression_="Norm_sqr(Trace(Grad(clambda*Trace(Green_Lagrange_E)*Id(2)+2*cmu*Green_Lagrange_E)(1,:,:)))", model=model,region=-1)
newVal_y=gf.asm_generic(mim=mim, order=0, _expression_="Norm_sqr(Trace(Grad(clambda*Trace(Green_Lagrange_E)*Id(2)+2*cmu*Green_Lagrange_E)(2,:,:)))", model=model,region=-1)
It does run without error but does it seem consistent to you regarding the indices, in particular ?
Here, the first line would be the computation of the integral over the whole domain of the squared norm of the x-component of Div(S), where S is a tensor, likewise for the second line for the y-component of Div(S).
Does it seem correct?

Best regards,
David.

Le mar. 14 juin 2022 à 18:57, <yves.renard@insa-lyon.fr> a écrit :
Dear David,

Did you try "Trace(Grad(Expr))" ? It works but with some limitations.

Best regards,

Yves

De: "David Danan" <daviddanan9021@gmail.com>
À: "Yves Renard" <yves.renard@insa-lyon.fr>
Envoyé: Mardi 14 Juin 2022 18:44:05
Objet: Re: [Getfem-users] Questions about user-defined strain energy

Dear Yves,

this time, it's perfect, the test passed.
Thanks a lot!

The only question i have left now is the possibility to compute the divergence of a tensor field ( arising from the gradient of an existing variable in a model), any suggestions?

Best regards,
David.

Le mar. 14 juin 2022 à 13:59, Renard Yves <yves.renard@insa-lyon.fr> a écrit :

Dear Danan,

I made an additional correction. Normally the tangent term of matrix_j2 is ok now.

Best regards,

Yves

Le 13/06/2022 à 17:37, David Danan a écrit :
Dear Yves,

i gave it a try, recompiled the master version available at https://git.savannah.nongnu.org/cgit/getfem.git/commit/ with your last fix
but it still fails to converge, unfortunately.

I have tested this with the program enclosed in the third email in this thread, it seems there is still an issue as far as i can tell.

Best regards,
David.

Le lun. 13 juin 2022 à 15:52, Renard Yves <yves.renard@insa-lyon.fr> a écrit :

Dear Danan and Kostas,

I had a rapid check to the implementation of Mooney-Rivlin hyperelastic law and matrix_j2 operators in getfem_nonlinear_elasticity.cc

I do not see any difference between the implementation and the _expression_ in the documentation (https://getfem.org/userdoc/model_nonlinear_elasticity.html#some-recalls-on-finite-strain-elasticity) for the Mooney-Rivlin law. However, it seems to me that a (det) is missing in the matrix_j2 gradient _expression_. I committed a fix. Could you try again if there is still a difference between the use of Mooney-Rivlin law and its _expression_ using matrix_j2 ?

Best regards,

Yves

Le 05/06/2022 à 20:01, David Danan a écrit :
Dear Kostas,

once again, thank you for pointing out where the problem is!
As soon as i have some time, it would be a pleasure to have a look and trying to fix this. Thanks a lot!

Regarding the question in my second message in this thread about the computation of the divergence of a tensor field, arising from the gradient of an existing variable in a model, do you have any idea?

BR,
David.

Le sam. 28 mai 2022 à 22:10, Konstantinos Poulios <logari81@googlemail.com> a écrit :
I see now that Yves has left a comment "to be verified":

// Second derivative
void second_derivative(const arg_list &args, size_type, size_type,
base_tensor &result) const { // To be verified

in the second derivative of matrix_j2_operator in getfem_nonlinear_elasticity.cc, so maybe you can help by checking/fixing this.

BR
Kostas

On Sat, May 28, 2022 at 10:05 PM Konstantinos Poulios <logari81@googlemail.com> wrote:
Dear David,

It works for me with
"paramsIMR(1)* ( Matrix_j1(F*F') - 3 )\
+ paramsIMR(2)* ( Matrix_i2(F*F')*pow(Det(F*F'),-2/3) - 3 )")
so you must have hit a bug in Matrix_j2.

BR
Kostas

On Sat, May 28, 2022 at 4:42 PM David Danan <daviddanan9021@gmail.com> wrote:
Dear Yves,

For the first question, at this point, i guess it's much more convenient at this point to share a script containing a minimal (not) working test regarding this difference.
I have kept the incompressibility brick in both cases.

For the second question, it seems clear, thanks.
For the third question, providing the _expression_ of the stress tensor directly instead is fine by me, it's not very difficult using the invariants already available in the GWFL combined with some macro.

While i am at it, i have another question.
It is somehow related to the one in the following thread but for the divergence of a tensor field

I would like to check the satisfaction of the equilibrium equation (i.e. div(P) + f0, where Pi is a stress tensor and f0 is a volumic force). To do so, i gave it a try without volumic forces
value=gf.asm_generic(mim=mim, order=0, _expression_="Norm_sqr(Div(clambda*Trace(Green_Lagrange_E)*Id(2)+2*cmu*Green_Lagrange_E))", model=model,region=-1)
where Green_Lagrange_E is a macro defined by
but it failed (Error in macro expansion. Only variable name are allowed for macro parameter preceded by Grad_ Hess_ Test_ or Test2_ prefixes)

Next, i tried to check whether i could at least apply the Div operator to Grad_u directly instead, like that
value=gf.asm_generic(mim=mim, order=0, _expression_="Norm_sqr(Div(Grad_u))", model=model,region=-1)
but it also failed, Grad_u was not recognized.

Is there a way to do that? I see that it's possible to have access to the hessian of the model variables; is rebuilding div(P) from the hessian the only way to compute this quantity?
Best regards,

David.

Le jeu. 26 mai 2022 à 10:17, <yves.renard@insa-lyon.fr> a écrit :
Dear David,

I do not see any reason why your _expression_ for the Mooney-Rivlin hyper elastic law would give a different result than the brick. At least if of cours you keep the
incompressible brick
(note that using sqr(expr) instead of pow(expr,2) is slightly more efficient).

Concerning the plane strain approximation, it is a priori possible to do this in GWFL, yes. You can define a 3D matrix from your 2D gradient and use the 3D hyperelastic law.

For the computation of the Von Mises stress from the potential, it is a little bit more tricky because the automatic differentiation of Getfem will give you the directional derivative (in the direction of the test function). So may be it is possible to extract the components with a specific choice of test functions, but it could be not so easy ... I never done that. Unfortunately it is preferable to give the _expression_ of the law in term of stress tensor directly.

Best regards,

Yves

De: "David Danan" <daviddanan9021@gmail.com>
À: "getfem-users" <getfem-users@nongnu.org>
Envoyé: Lundi 23 Mai 2022 13:41:14
Objet: [Getfem-users] Questions about user-defined strain energy

Dear Getfem-users,

to check whether i understand correctly the implementation (and because it's actually much clearer in my code that way), i am trying to replace some predefined bricks for strain energy by the equivalent version using GWFL.
For that, i have several questions

First question
For instance
lawname = 'SaintVenant Kirchhoff'
clambda,cmu = params["clambda"],params["cmu"]
idBrick=model.add_finite_strain_elasticity_brick(mim, lawname, 'u', 'paramsSVK')
becomes
clambda,cmu = params["clambda"],params["cmu"]

And the results seems to be the same.
However, for the incompressible Mooney Rivlin strain energy case, it does not behave as i expected

Using the example there as a basis in a 3D case
and after some simplifications

I have tried to replace the last line in
lawname = 'Incompressible Mooney Rivlin'
model.add_finite_strain_elasticity_brick(mim, lawname, 'u', 'paramsIMR')
By this
Which was exactly the same, of course. Next, i have tried to replace it by these lines
model.add_nonlinear_generic_assembly_brick(mim, "paramsIMR(1)* ( Matrix_j1(Right_Cauchy_Green(F)) - 3 )+ paramsIMR(2)* ( Matrix_j2(Right_Cauchy_Green(F)) - 3 )")
Which failed to converge

I thought this _expression_ was consistent with the implementation of Mooney_Rivlin_hyperelastic_law::strain_energy and the compute_invariants here
and the documentation
But, clearly, i am missing something. Could you explain what I am doing wrong?

Second question
In a 2D case, i would like to be able to use either a plane strain approximation based on a given strain energy _expression_ in 3D.
In the implementation, it is nicely done in  scalar_type plane_strain_hyperelastic_law::strain_energy there

Is it possible to do the same using the GWFL?

Third question
I would like to be able to compute the von-mises field for any strain energy function.

If It is an existing brick, the method
md.compute_finite_strain_elasticity_Von_Mises(lawname, 'u', 'params', mfdu)
will do the trick just fine.

If it's not the case, i guess i can use something akin to the actual implementation of compute_finite_strain_elasticity_Von_Mises

std::string expr = "sqrt(3/2)*Norm(Deviator(Cauchy_stress_from_PK2("
+ adapted_lawname + "_PK2(Grad_" + varname + "," + params + "),Grad_"
+ varname + ")))";
ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM, rg);
}

combined with local_projection to get the value at the elements.
The question is: is it possible to compute Piola Kirchhoff 2 from the strain energy within the GWFL _expression_ given to local_projection?
I have the impression it's the only thing i need to be able to do what i want.
Let W be a strain energy function and E be the Green-Lagrange tensor (defined as macros, let's say), is ```Diff(W, E)``` doing exactly what i am expecting for this purpose?
If it's not the case, is there another way to do it?

Thank you in advance,
kind regards,
David.

reply via email to