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