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
<logar...@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
<logar...@googlemail.com> wrote:
Dear David,
It works for me with
model.add_nonlinear_generic_assembly_brick(mim,
"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
<daviddanan9...@gmail.com> wrote:
Dear Yves,
first, thanks for your detailed answers!
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
https://lists.nongnu.org/archive/html/getfem-users/2015-03/msg00009.html
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
model.add_macro("epsilon", "(Grad_u+(Grad_u)')/2")
model.add_macro("Green_Lagrange_E", "epsilon+
(Grad_u)'*Grad_u/2")
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.ren...@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
md.add_finite_strain_incompressibility_brick(mim,
'u', 'p')
(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" <daviddanan9...@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"]
model.add_initialized_data('paramsSVK', [clambda,
cmu])
idBrick=model.add_finite_strain_elasticity_brick(mim,
lawname, 'u', 'paramsSVK')
becomes
clambda,cmu= params["clambda"],params["cmu"]
model.add_initialized_data('paramsSVK', [clambda,
cmu])
idBrick=model.add_nonlinear_generic_assembly_brick(mim,"(paramsSVK(1)/2)*pow(Trace(Green_Lagrange_E),2)+paramsSVK(2)*Trace(Green_Lagrange_E'*Green_Lagrange_E)")
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
https://github.com/getfem-doc/getfem/blob/master/interface/tests/python/demo_nonlinear_elasticity.py
and after some simplifications
I have tried to replace the last line in
lawname= 'Incompressible Mooney Rivlin'
model.add_initialized_data('paramsIMR', [1,1])
model.add_finite_strain_elasticity_brick(mim,
lawname, 'u', 'paramsIMR')
By this
model.add_nonlinear_generic_assembly_brick(mim,
"Incompressible_Mooney_Rivlin_potential(Grad_u,
[paramsIMR(1);paramsIMR(2)])")
Which was exactly the same, of course. Next, i
have tried to replace it by these lines
model.add_macro("F", "Grad_u+Id(3)")
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
http://download-mirror.savannah.gnu.org/releases/getfem/doc/getfem_reference/getfem__nonlinear__elasticity_8cc_source.html
and the documentation
https://getfem.org/userdoc/model_nonlinear_elasticity.html?highlight=finite%20strain
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
http://download-mirror.savannah.gnu.org/releases/getfem/doc/getfem_reference/getfem__nonlinear__elasticity_8cc_source.html
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
<http://download-mirror.savannah.gnu.org/releases/getfem/doc/getfem_reference/namespacegetfem.html#a9166bc2065a4f2b5633e2e933ec3e693>
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.
