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=InnerRadiussin(t); label=1;}
border OuterEdge(t=0, 2.pi){x=(1.0 + 0.0InnerEllipseExtension)OuterRadiuscos(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)( (u1v1 + u2v2) / 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 * (w1N.x + w2N.y)
)
//Load vectors tangential to the boundary:
//- int1d(Th, 1)(
// Pa * (w1N.y - w2N.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.