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 <daviddanan9...@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
> 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.
>>>>>>>>
>>>>>>>>
>>
