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_Lagrangian(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 <[email protected]> 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. > >
