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
model.add_macro("epsilon", "(Grad_u+(Grad_u)')/2")
model.add_macro("Green_Lagrange_E", "epsilon+ (Grad_u)'*Grad_u/2")
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.ren...@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" <daviddanan9...@gmail.com>
> *À: *"Yves Renard" <yves.ren...@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.ren...@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.ren...@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 <
>>> 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.
>>>>>>>
>>>>>>>
>
