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