# The example in FreeFem Documentation about Neo-Hookean Compressible Material

**URL:** <https://community.freefem.org/t/the-example-in-freefem-documentation-about-neo-hookean-compressible-material/3044>\
**Category:** General Discussion\
**Created:** [March 16, 2024, 8:29am UTC](https://community.freefem.org/t/the-example-in-freefem-documentation-about-neo-hookean-compressible-material/3044 "2024-03-16T08:29:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![NYKong](https://avatars.discourse-cdn.com/v4/letter/n/8491ac/32.png) [@NYKong](https://community.freefem.org/u/NYKong)\
**Post date:** [March 16, 2024, 8:29am UTC](https://community.freefem.org/t/the-example-in-freefem-documentation-about-neo-hookean-compressible-material/3044/1 "2024-03-16T08:29:14Z")

</div>

Hi , the problem is that the program is misbehaved in the official FreeFem-Document (page 636), which is an example about a Neo-Hookean Compressible Material. I meet errors in line 9 and other definitions of macro. The whole code as follows(totally copy from the official document):  
// Macro  
//Macros for the gradient of a vector field (u1, u2)  
macro grad11(u1, u2) (dx(u1)) //  
macro grad21(u1, u2) (dy(u1)) //  
macro grad12(u1, u2) (dx(u2)) //  
macro grad22(u1, u2) (dy(u2)) //

//Macros for the deformation gradient  
macro F11(u1, u2) (1.0 + grad11(u1, u2)) //  
macro F12(u1, u2) (0.0 + grad12(u1, u2)) //  
macro F21(u1, u2) (0.0 + grad21(u1, u2)) //  
macro F22(u1, u2) (1.0 + grad22(u1, u2)) //

//Macros for the incremental deformation gradient  
macro dF11(varu1, varu2) (grad11(varu1, varu2)) //  
macro dF12(varu1, varu2) (grad12(varu1, varu2)) //  
macro dF21(varu1, varu2) (grad21(varu1, varu2)) //  
macro dF22(varu1, varu2) (grad22(varu1, varu2)) //

//Macro for the determinant of the deformation gradient  
macro J(u1, u2) (  
F11(u1, u2)\*F22(u1, u2)  
- F12(u1, u2)\*F21(u1, u2)  
) //

//Macros for the inverse of the deformation gradient  
macro Finv11 (u1, u2) (  
F22(u1, u2) / J(u1, u2)  
) //  
macro Finv22 (u1, u2) (  
F11(u1, u2) / J(u1, u2)  
) //  
macro Finv12 (u1, u2) (  
- F12(u1, u2) / J(u1, u2)  
) //  
macro Finv21 (u1, u2) (  
- F21(u1, u2) / J(u1, u2)  
) //

//Macros for the square of the inverse of the deformation gradient  
macro FFinv11 (u1, u2) (  
Finv11(u1, u2)^2  
+ Finv12(u1, u2)\*Finv21(u1, u2)  
) //

macro FFinv12 (u1, u2) (  
Finv12(u1, u2)\*(Finv11(u1, u2)  
+ Finv22(u1, u2))  
) //

macro FFinv21 (u1, u2) (  
Finv21(u1, u2)\*(Finv11(u1, u2)  
+ Finv22(u1, u2))  
) //

macro FFinv22 (u1, u2) (  
Finv12(u1, u2)\*Finv21(u1, u2)  
+ Finv22(u1, u2)^2  
) //

//Macros for the inverse of the transpose of the deformation gradient  
macro FinvT11(u1, u2) (Finv11(u1, u2)) //  
macro FinvT12(u1, u2) (Finv21(u1, u2)) //  
macro FinvT21(u1, u2) (Finv12(u1, u2)) //  
macro FinvT22(u1, u2) (Finv22(u1, u2)) //

//The left Cauchy-Green strain tensor  
macro B11(u1, u2) (  
F11(u1, u2)^2 + F12(u1, u2)^2  
)//

macro B12(u1, u2) (  
F11(u1, u2)\*F21(u1, u2)  
+ F12(u1, u2)\*F22(u1, u2)  
)//

macro B21(u1, u2) (  
F11(u1, u2)\*F21(u1, u2)  
+ F12(u1, u2)\*F22(u1, u2)  
)//

macro B22(u1, u2)(  
F21(u1, u2)^2 + F22(u1, u2)^2  
)//

//The macros for the auxiliary tensors (D0, D1, D2, …): Begin  
//The tensor quantity D0 = F{n} (delta F{n+1})^T  
macro d0Aux11 (u1, u2, varu1, varu2) (  
dF11(varu1, varu2) \* F11(u1, u2)  
+ dF12(varu1, varu2) \* F12(u1, u2)  
)//

macro d0Aux12 (u1, u2, varu1, varu2) (  
dF21(varu1, varu2) \* F11(u1, u2)  
+ dF22(varu1, varu2) \* F12(u1, u2)  
)//

macro d0Aux21 (u1, u2, varu1, varu2) (  
dF11(varu1, varu2) \* F21(u1, u2)  
+ dF12(varu1, varu2) \* F22(u1, u2)  
)//

macro d0Aux22 (u1, u2, varu1, varu2) (  
dF21(varu1, varu2) \* F21(u1, u2)  
+ dF22(varu1, varu2) \* F22(u1, u2)  
)//

///The tensor quantity D1 = D0 + D0^T  
macro d1Aux11 (u1, u2, varu1, varu2) (  
2.0 \* d0Aux11 (u1, u2, varu1, varu2)  
)//

macro d1Aux12 (u1, u2, varu1, varu2) (  
d0Aux12 (u1, u2, varu1, varu2)  
+ d0Aux21 (u1, u2, varu1, varu2)  
)//

macro d1Aux21 (u1, u2, varu1, varu2) (  
d1Aux12 (u1, u2, varu1, varu2)  
)//

macro d1Aux22 (u1, u2, varu1, varu2)(  
2.0 \* d0Aux22 (u1, u2, varu1, varu2)  
)//

///The tensor quantity D2 = F^{-T}_{n} dF_{n+1}  
macro d2Aux11 (u1, u2, varu1, varu2) (  
dF11(varu1, varu2) \* FinvT11(u1, u2)  
+ dF21(varu1, varu2) \* FinvT12(u1, u2)  
)//

macro d2Aux12 (u1, u2, varu1, varu2) (  
dF12(varu1, varu2) \* FinvT11(u1, u2)  
+ dF22(varu1, varu2) \* FinvT12(u1, u2)  
)//

macro d2Aux21 (u1, u2, varu1, varu2) (  
dF11(varu1, varu2) \* FinvT21(u1, u2)  
+ dF21(varu1, varu2) \* FinvT22(u1, u2)  
)//

macro d2Aux22 (u1, u2, varu1, varu2) (  
dF12(varu1, varu2) \* FinvT21(u1, u2)  
+ dF22(varu1, varu2) \* FinvT22(u1, u2)  
)//

///The tensor quantity D3 = F^{-2}_{n} dF_{n+1}  
macro d3Aux11 (u1, u2, varu1, varu2, w1, w2) (  
dF11(varu1, varu2) \* FFinv11(u1, u2) \* grad11(w1, w2)  
+ dF21(varu1, varu2) \* FFinv12(u1, u2) \* grad11(w1, w2)  
+ dF11(varu1, varu2) \* FFinv21(u1, u2) \* grad12(w1, w2)  
+ dF21(varu1, varu2) \* FFinv22(u1, u2) \* grad12(w1, w2)  
)//

macro d3Aux12 (u1, u2, varu1, varu2, w1, w2) (  
dF12(varu1, varu2) \* FFinv11(u1, u2) \* grad11(w1, w2)  
+ dF22(varu1, varu2) \* FFinv12(u1, u2) \* grad11(w1, w2)  
+ dF12(varu1, varu2) \* FFinv21(u1, u2) \* grad12(w1, w2)  
+ dF22(varu1, varu2) \* FFinv22(u1, u2) \* grad12(w1, w2)  
)//

macro d3Aux21 (u1, u2, varu1, varu2, w1, w2) (  
F11(varu1, varu2) \* FFinv11(u1, u2) \* grad21(w1, w2)  
+ dF21(varu1, varu2) \* FFinv12(u1, u2) \* grad21(w1, w2)  
+ dF11(varu1, varu2) \* FFinv21(u1, u2) \* grad22(w1, w2)  
+ dF21(varu1, varu2) \* FFinv22(u1, u2) \* grad22(w1, w2)  
)//

macro d3Aux22 (u1, u2, varu1, varu2, w1, w2) (  
dF12(varu1, varu2) \* FFinv11(u1, u2) \* grad21(w1, w2)  
+ dF22(varu1, varu2) \* FFinv12(u1, u2) \* grad21(w1, w2)  
+ dF12(varu1, varu2) \* FFinv21(u1, u2) \* grad22(w1, w2)  
+ dF22(varu1, varu2) \* FFinv22(u1, u2) \* grad22(w1, w2)  
)//

///The tensor quantity D4 = (grad w) \* Finv  
macro d4Aux11 (w1, w2, u1, u2) (  
Finv11(u1, u2)\*grad11(w1, w2)  
+ Finv21(u1, u2)\*grad12(w1, w2)  
)//

macro d4Aux12 (w1, w2, u1, u2) (  
Finv12(u1, u2)\*grad11(w1, w2)  
+ Finv22(u1, u2)\*grad12(w1, w2)  
)//

macro d4Aux21 (w1, w2, u1, u2) (  
Finv11(u1, u2)\*grad21(w1, w2)  
+ Finv21(u1, u2)\*grad22(w1, w2)  
)//

macro d4Aux22 (w1, w2, u1, u2) (  
Finv12(u1, u2)\*grad21(w1, w2)  
+ Finv22(u1, u2)\*grad22(w1, w2)  
)//  
//The macros for the auxiliary tensors (D0, D1, D2, …): End

//The macros for the various stress measures: BEGIN  
//The Kirchhoff stress tensor  
macro StressK11(u1, u2) (  
mu \* (B11(u1, u2) - 1.0)  
)//

//The Kirchhoff stress tensor  
macro StressK12(u1, u2) (  
mu \* B12(u1, u2)  
)//

//The Kirchhoff stress tensor  
macro StressK21(u1, u2) (  
mu \* B21(u1, u2)  
)//

//The Kirchhoff stress tensor  
macro StressK22(u1, u2) (  
mu \* (B22(u1, u2) - 1.0)  
)//

//The tangent Kirchhoff stress tensor  
macro TanK11(u1, u2, varu1, varu2) (  
mu \* d1Aux11(u1, u2, varu1, varu2)  
)//

macro TanK12(u1, u2, varu1, varu2) (  
mu \* d1Aux12(u1, u2, varu1, varu2)  
)//

macro TanK21(u1, u2, varu1, varu2) (  
mu \* d1Aux21(u1, u2, varu1, varu2)  
)//

macro TanK22(u1, u2, varu1, varu2) (  
mu \* d1Aux22(u1, u2, varu1, varu2)  
)//  
//The macros for the stress tensor components: END

// Parameters  
real mu = 5.e2; //Elastic coefficients (kg/cm^2)  
real D = 1.e3; //(1 / compressibility)  
real Pa = -3.e2; //Stress loads

real InnerRadius = 1.e0; //The wound radius  
real OuterRadius = 4.e0; //The outer (truncated) radius  
real tol = 1.e-4; //Tolerance (L^2)  
real InnerEllipseExtension = 1.e0; //Extension of the inner ellipse ((major axis)(minor axis))

int m = 40, n = 20;

// Mesh  
border InnerEdge(t=0, 2._pi){x=(1.0 + InnerEllipseExtension)InnerRadiuscos(t);y=InnerRadius_sin(t); label=1;}  
border OuterEdge(t=0, 2._pi){x=(1.0 + 0.0_InnerEllipseExtension)_OuterRadius_cos(t);y=OuterRadius\*sin(t); label=2;}  
mesh Th = buildmesh(InnerEdge(-m) + OuterEdge(n));  
int bottom = 1, right = 2, upper = 3, left = 4;  
plot(Th);

// Fespace  
fespace Wh(Th, P1dc);  
fespace Vh(Th, [P1, P1]);  
Vh [w1, w2], [u1n,u2n], [varu1, varu2];  
Vh [ehat1x, ehat1y], [ehat2x, ehat2y];  
Vh [auxVec1, auxVec2]; //The individual elements of the total 1st Piola-Kirchoff stress  
Vh [ef1, ef2];

fespace Sh(Th, P1);  
Sh p, ppp;  
Sh StrK11, StrK12, StrK21, StrK22;  
Sh u1, u2;

// Problem  
varf vfMass1D(p, q) = int2d(Th)(p\*q);  
matrix Mass1D = vfMass1D(Sh, Sh, solver=CG);

p = 1;  
ppp = Mass1D \* p;

real DomainMass = ppp.sum;  
cout \<\< "DomainMass = " \<\< DomainMass \<\< endl;

varf vmass ([u1, u2], [v1, v2], solver=CG)  
= int2d(Th)( (u1_v1 + u2_v2) / DomainMass );

matrix Mass = vmass(Vh, Vh);

matrix Id = vmass(Vh, Vh);

//Define the standard Euclidean basis functions  
[ehat1x, ehat1y] = [1.0, 0.0];  
[ehat2x, ehat2y] = [0.0, 1.0];

real ContParam, dContParam;

problem neoHookeanInc ([varu1, varu2], [w1, w2], solver=LU)  
= int2d(Th, qforder=1)(  
- (  
StressK11 (u1n, u2n) \* d3Aux11(u1n, u2n, varu1, varu2, w1, w2)  
+ StressK12 (u1n, u2n) \* d3Aux12(u1n, u2n, varu1, varu2, w1, w2)  
+ StressK21 (u1n, u2n) \* d3Aux21(u1n, u2n, varu1, varu2, w1, w2)  
+ StressK22 (u1n, u2n) \* d3Aux22(u1n, u2n, varu1, varu2, w1, w2)  
)  
+ TanK11 (u1n, u2n, varu1, varu2) \* d4Aux11(w1, w2, u1n, u2n)  
+ TanK12 (u1n, u2n, varu1, varu2) \* d4Aux12(w1, w2, u1n, u2n)  
+ TanK21 (u1n, u2n, varu1, varu2) \* d4Aux21(w1, w2, u1n, u2n)  
+ TanK22 (u1n, u2n, varu1, varu2) \* d4Aux22(w1, w2, u1n, u2n)  
)  
+ int2d(Th, qforder=1)(  
StressK11 (u1n, u2n) \* d4Aux11(w1, w2, u1n, u2n)  
+ StressK12 (u1n, u2n) \* d4Aux12(w1, w2, u1n, u2n)  
+ StressK21 (u1n, u2n) \* d4Aux21(w1, w2, u1n, u2n)  
+ StressK22 (u1n, u2n) \* d4Aux22(w1, w2, u1n, u2n)  
)  
//Choose one of the following two boundary conditions involving Pa:  
// Load vectors normal to the boundary:  
- int1d(Th, 1)(  
Pa \* (w1_N.x + w2_N.y)  
)  
//Load vectors tangential to the boundary:  
//- int1d(Th, 1)(  
// Pa \* (w1_N.y - w2_N.x)  
//)  
+ on(2, varu1=0, varu2=0)  
;

//Auxiliary variables  
matrix auxMat;

// Newton’s method  
ContParam = 0.;  
dContParam = 0.01;

//Initialization:  
[varu1, varu2] = [0., 0.];  
[u1n, u2n] = [0., 0.];  
real res = 2. \* tol;  
real eforceres;  
int loopcount = 0;  
int loopmax = 45;

// Iterations  
while (loopcount \<= loopmax && res \>= tol){  
loopcount ++;  
cout \<\< "Loop " \<\< loopcount \<\< endl;

```
// Solve
neoHookeanInc;

// Update
u1 = varu1;
u2 = varu2;

// Residual
w1[] = Mass*varu1[];
res = sqrt(w1[]' * varu1[]); //L^2 norm of [varu1, varu2]
cout << " L^2 residual = " << res << endl;

// Newton
u1n[] += varu1[];

// Plot
plot([u1n,u2n], cmm="displacement");

```

}

// Plot  
plot(Th, wait=true);

// Movemesh  
mesh Th1 = movemesh(Th, [x+u1n, y+u2n]);

// Plot  
plot(Th1, wait=true);  
plot([u1n,u2n]);

What can I do to debug this code? Thanks so much.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [March 17, 2024, 4:34pm UTC](https://community.freefem.org/t/the-example-in-freefem-documentation-about-neo-hookean-compressible-material/3044/2 "2024-03-17T16:34:58Z")

</div>

Hello,  
maybe you got copy-paste errors, and it is true that the code in freefem doc has also some copy-paste errors. Here is a file that seems to work  
[neo-hookean.edp](https://community.freefem.org/uploads/short-url/iG5S33gTGUJGtbCxwdcgH1Ib3rl.edp) (9.1 KB)