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.
>
>
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import getfem as gf
gf.util('trace level', 1)
def SolveBasic3DMecaProblem(version):
N1 = 2; N2 = 4; h = 20.; DX = 1./N1; DY = (1.*h)/N2;
m = gf.Mesh('cartesian', np.arange(-0.5, 0.5+DX,DX), np.arange(0., h+DY,DY),
np.arange(-1.5, 1.5+3*DX,3*DX))
mfu = gf.MeshFem(m, 3)
mfu.set_fem(gf.Fem('FEM_QK(3,2)'))
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS_PARALLELEPIPED(3,4)'))
ftop = m.outer_faces_with_direction([0., 1., 0.], 0.5)
fbot = m.outer_faces_with_direction([0., -1., 0.], 0.5)
m.set_region(1, ftop)
m.set_region(2, fbot)
m.set_region(3, np.append(ftop,fbot,axis=1));
model = gf.Model('real')
model.add_fem_variable('u', mfu)
model.add_initialized_data('paramsIMR', [1,1])
if version=="Predefined":
lawname = 'Incompressible Mooney Rivlin'
model.add_finite_strain_elasticity_brick(mim, lawname, 'u', 'paramsIMR')
elif version=="AssemblyPredefined":
model.add_nonlinear_generic_assembly_brick(mim, "Incompressible_Mooney_Rivlin_potential(Grad_u, [paramsIMR(1);paramsIMR(2)])")
elif version=="GWFL":
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 )")
else:
raise Exception("Version not available")
#Add finite strain incompressibility condition
mfp = gf.MeshFem(m,1)
mfp.set_classical_discontinuous_fem(1)
model.add_fem_variable('p', mfp)
model.add_finite_strain_incompressibility_brick(mim, 'u', 'p')
model.add_initialized_data('f', [0,-0.1,0])
model.add_source_term_brick(mim, 'u', 'f')
model.add_initialized_data('u_d', [0,0,0])
model.add_Dirichlet_condition_with_multipliers(mim, 'u', 1, 3, "u_d")
model.variable_list()
maxIter=50
solverIt=model.solve('very noisy','max_iter',maxIter, 'max_res',1e-7, 'lsearch','simplest')
if not solverIt[0]<maxIter:
raise Exception("Failed to converge")
displacements = model.variable('u')
return displacements
def CheckConsistency_GWFL_PredefinedBrick():
print("Checking consitency between Predefined brick and AssemblyPredefined brick")
solutionPredefined=SolveBasic3DMecaProblem(version="Predefined")
solutionAssemblyPredefined=SolveBasic3DMecaProblem(version="AssemblyPredefined")
np.testing.assert_almost_equal(solutionPredefined,solutionAssemblyPredefined)
print("Checking consitency between Predefined brick and GWFL")
solutionGWFL=SolveBasic3DMecaProblem(version="GWFL")
np.testing.assert_almost_equal(solutionPredefined,solutionGWFL)
if __name__ == '__main__':
CheckConsistency_GWFL_PredefinedBrick()
