Re: [Getfem-users] Hyperelastic law Benchmark

[franz chouly]

Dear David,

Maybe using a specific variant of Newton can also help (see, e.g., 
https://scholar.harvard.edu/files/vasios/files/ArcLength.pdf) ? 

Best regards,

  Franz


----- Mail original -----
> De: "Yves Renard" <address@hidden>
> À: "David Danan" <address@hidden>
> Cc: address@hidden
> Envoyé: Mercredi 28 Juin 2017 20:03:42
> Objet: Re: [Getfem-users] Hyperelastic law Benchmark

> Dear David,
> 
> No, no idea. Often, if the load is too large, the Newton fails to find the
> solution. It is quite normal. You can decompose into small load steps.
> Alternatively, you can also use the (beautiful) continuations tools of Getfem
> developed by Tomas Ligursky which automatically adapt the step to try to find 
> a
> solution. If the continuation procedure do not find a solution, this may be of
> course due to the fact that there could be no solution to the problem for a
> large load ... This could happen for very large loads.
> 
> Best regards,
> 
> Yves.
> 
> ----- Mail original -----
> De: "David Danan" <address@hidden>
> À: "logari81" <address@hidden>
> Cc: address@hidden
> Envoyé: Jeudi 22 Juin 2017 12:27:23
> Objet: Re: [Getfem-users] Hyperelastic law Benchmark
> 
> Dear getfem users,
> 
> i encounter stability issues with the second problem of the same benchmark.
> http://rspa.royalsocietypublishing.org/content/471/2184/20150641#sec-15
> 
> By using [39], i have tried to improve the stability but the computation
> still failed at the first iteration anyway.
> 
> So far i can use either:
> -The classical version of the hyperelastic law with an incompressibility
> condition (lagrange multiplier)
> -The classical version of the hyperelastic law with a
> quasi-incompressibility condition (penalization)
> -A stabilized version of the law (incompressibility with Lagrange
> multiplier and penalization)
> -An isochoric version that use the isochoric component of the deformation
> gradient with either an incompressibility condition, a penalization or both
> of them.
> 
> Note that the convergence is reached if i decrease the value of the
> pressure.
> 
> Please, find enclosed:
> -the program Ellipsoid.cc and the parameter file Ellipsoid.param
> -Two meshes provided by the benchmark. I have also tried with several other
> meshes, by using the methods described above, with the same result.
> 
> Do you have any idea/advices regarding these issues?
> Thanks a lot.
> 
> Yours sincerely,
> David.
> 
> 
> 
> 
> 2017-06-15 20:53 GMT+02:00 David Danan <address@hidden>:
> 
>> Dear Kostas,
>>
>> first, thank you for your very fast answer.
>>
>> So far, i didn't even think about using the high level generic assembly
>> language for this (and it's a shame given how easy it is to use in this
>> configuration).
>> I followed your suggestion and tried your expression (which seems
>> perfectly correct to me) and the deformation was, as far as i could tell,
>> identical to the one in the Results_nosym.png.
>>
>> However, your second remark was spot on because the paper actually said it
>> explicitly:
>> "Please note that in all problems the direction of the pressure boundary
>> condition changes with the deformed surface orientation, and its magnitude
>> scales with the deformed area."
>>
>> Therefore, thanks a lot for providing your notes here, it is very
>> helpful...And it seems to actually work just fine! There is no visible
>> difference between the deformed configuration given by Getfem and the
>> reference solution with this modification.
>>
>> Once again, thank you!
>>
>> David.
>>
>>
>> 2017-06-14 17:52 GMT+02:00 Konstantinos Poulios <address@hidden>:
>>
>>> It doesn't seem to be explicitly stated in the paper but since it is
>>> about cardiovascular simulation I guess that applying the pressure as a
>>> follower load is the standard thing to do, so I am providing you here with
>>> my notes on follower loads:
>>>
>>> [image: Inline image 1]
>>>
>>>
>>> So if your p is the actual blood pressure you need the upper right case
>>> of the table with q=-p. If at some point you also need shear stresses from
>>> the fluid you can also use the second row of the table.
>>>
>>> BR
>>> Kostas
>>>
>>>
>>> On Wed, Jun 14, 2017 at 4:57 PM, Konstantinos Poulios <
>>> address@hidden> wrote:
>>>
>>>> Dear David
>>>>
>>>> Have you tried the high level generic assembly language for this?
>>>>
>>>> In principle you should be able to provide GetFEM with your energy
>>>> density function and let GetFEM do the necessary derivations.
>>>>
>>>> Instead of
>>>>
>>>> getfem::add_finite_strain_elasticity_brick(...)
>>>>
>>>> you have to call
>>>>
>>>>  getfem::add_nonlinear_generic_assembly_brick
>>>> (md, mim, "0.5*C*(exp([[bf,bfs,bfs],[bfs,bt,bt],[bfs,bt,bt]]:(Green_La
>>>> grangian(Grad_u+Id(3)).*Green_Lagrangian(Grad_u+Id(3))))-1)");
>>>>
>>>> with C,bf,bfs and bt scalar parameters defined with
>>>> md.add_initialized_scalar_data(...).
>>>>
>>>> I hope I got the expression from the paper right. Can you give it a try?
>>>>
>>>> Then the other question is how the applied surface pressure p is
>>>> distributed, if it is a follower load you need a more complex expression
>>>> than
>>>>
>>>> "-Pressure*Normal.Test_u"
>>>>
>>>> Because "Normal" is in the undeformed configuration.
>>>>
>>>> BR
>>>> Kostas
>>>>
>>>>
>>>> On Wed, Jun 14, 2017 at 3:58 PM, David Danan <address@hidden>
>>>> wrote:
>>>>
>>>>> Dear Getfem users,
>>>>>
>>>>> i am trying to implement a new hyperelastic law and, in order to
>>>>> validate my results, i am using the following Benchmark
>>>>> http://rspa.royalsocietypublishing.org/content/471/2184/20150641#sec-15
>>>>>
>>>>> There are 3 problems, for now i am working on the first one that is to
>>>>> say the deformation of a 3D rectangular beam clamped on one side and with 
>>>>> a
>>>>> pressure applied to the bottom face.
>>>>>
>>>>> While the deformed configuration given by Getfem is relatively close to
>>>>> the reference(s) solution(s) provided by the benchmark, a visible
>>>>> difference between them still remains and i don't understand where it 
>>>>> comes
>>>>> from.
>>>>>
>>>>> The material is governed by a transversely isotropic constitutive law
>>>>> with an incompressibility constraint, often used in cardiac modelling,
>>>>> where the strain energy function is a function of the components of the
>>>>> Green–Lagrange strain tensor *E.*
>>>>>
>>>>> I tried 2 differents implementations of this law:
>>>>> -the first use the symmetry of the Green-Lagrange strain tensor to
>>>>> simplify the strain energy function
>>>>> -The second does not (ergo it is necessary to write the 9 components of
>>>>> the second piola Kirchhoff stress tensor and the 81 components of the
>>>>> fourth order tensor)
>>>>>
>>>>> Please find enclosed
>>>>> -the comparison in the first case: Results.png
>>>>> -the comparison in the second case: Results_nosym.png (slightly better
>>>>> results but 15 times as slow as the first version)
>>>>> -the python program used to compute the derivative and second
>>>>> derivative of the strain energy function in the first case.
>>>>> -the implementation of the laws in getfem_nonlinear_elasticity.cc and
>>>>> getfem_nonlinear_elasticity.h
>>>>> -The program Guccione.cc and Guccione.param used to produce these very
>>>>> pictures
>>>>>
>>>>> in both pictures, the reference solution is in grey.
>>>>> The computation uses Q2/Q1 elements (displacement/lagrange multiplier),
>>>>> since there is no restrictions regarding these aspects.
>>>>> I have tried with a quasi-incompressibility condition instead of the
>>>>> Lagrangian multiplier: same result (which was to be expected).
>>>>> I have also tried with other meshes (more or less refined) used by
>>>>> other teams but in vain.
>>>>>
>>>>> Could someone have a look and provide some advices regarding this
>>>>> case/tell me what i am doing wrong?
>>>>>
>>>>> Thanks a lot.
>>>>>
>>>>> Yours sincerely,
>>>>> David.
>>>>>
>>>>>
>>>>
>>>