# Traction Forces - BilinearOperator NM error

**URL:** <https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857>\
**Category:** General Discussion\
**Created:** [April 11, 2025, 3:31pm UTC](https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857 "2025-04-11T15:31:20Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![prayuktha](https://avatars.discourse-cdn.com/v4/letter/p/2acd7d/32.png) [@prayuktha](https://community.freefem.org/u/prayuktha)\
**Post date:** [April 11, 2025, 3:31pm UTC](https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857/1 "2025-04-11T15:31:21Z")

</div>

Hi all, I am writing a code to compute the displacement field in Contact Mechanics. But Its giving an error I cant figure out. Can anyone help me? Thanks in Advance. The code is the following one:  
1 : load “msh3”  
2 :  
3 : // Parameter arrays  
4 : real[int] Pvals = [5e4, 1e5, 2e5]; // Load (N)  
5 : real[int] avals = [0.005, 0.01, 0.015]; // Contact radius (m)  
6 : real[int] Evals = [1e9, 7e10, 2e11]; // Young’s modulus ¶  
7 : real[int] muvals = [0.25, 0.3, 0.35]; // Poisson’s ratio  
8 : real[int] Rvals = [0.01, 0.03, 0.05]; // Sphere radius (m)  
9 :  
10 : int numR = 50; // Number of radial sampling points  
11 :  
12 : // Print header for output  
13 : cout \<\< “P,E,mu,a,R”;  
14 : for (int i = 0; i \< numR; ++i)  
15 : cout \<\< “,uz” + i;  
16 : cout \<\< “\n”;  
17 :  
18 : // Looping over parameters using classic for loops  
19 : for (int ip = 0; ip \< Pvals.n; ++ip)  
20 : for (int ia = 0; ia \< avals.n; ++ia)  
21 : for (int ie = 0; ie \< Evals.n; ++ie)  
22 : for (int imu = 0; imu \< muvals.n; ++imu)  
23 : for (int ir = 0; ir \< Rvals.n; ++ir)  
24 : {  
25 : real loadValue = Pvals[ip];  
26 : real aparam = avals[ia];  
27 : real E = Evals[ie];  
28 : real mu = muvals[imu];  
29 : real R = Rvals[ir];  
30 :  
31 : real L = 2 \* R; // domain width  
32 : real H = R; // domain height  
33 :  
34 : // Mesh  
35 : int nx = 50, ny = 50;  
36 : mesh Th = square(nx, ny, [x_L, y_H]);  
37 :  
38 : // FEM space  
39 : fespace Vh(Th, P1);  
40 : Vh ur, uz, vr, vz;  
41 :  
42 : // Macros for strain components  
43 : macro epsrr(ur, uz) dx(ur) ) //  
44 : macro epszz(ur, uz) dy(uz) ) //  
45 : macro epsrz(ur, uz) (dy(ur) + dx(uz)) / 2 ) //  
46 : macro divu(ur, uz) (dx(ur) + dy(uz)) ) // total divergence  
47 :  
48 : // Lame parameters  
49 : real lambda = E \* mu / ((1. + mu) \* (1. - 2. \* mu));  
50 : real muL = E / (2. \* (1. + mu));  
51 :  
52 : // Variational form (bilinear form)  
53 : varf avar([ur, uz], [vr, vz]) =  
54 : int2d(Th)(  
55 : lambda \* divu(ur, uz) (dx(ur) + dy( uz)) \* divu(vr, vz) (dx(vr) + dy( vz))  
56 : + 2. \* muL \* (  
57 : epsrr(ur, uz) dx(ur) \* epsrr(vr, vz) dx(vr)  
58 : + epszz(ur, uz) dy( uz) \* epszz(vr, vz) dy( vz)  
59 : + 2. \* epsrz(ur, uz) (dy(ur) + dx( uz)) / 2 \* epsrz(vr, vz) (dy(vr) + dx( vz)) / 2  
60 : )  
61 : )  
62 : + on(1, ur = 0, uz = 0); // Apply zero displacement on boundary 1  
63 :  
64 : // Load vector  
65 : varf Lrhs([vr, vz], [vrtest, vztest]) =  
66 : int1d(Th, 4)(  
67 : (x \<= aparam) \* (-loadValue / (pi \* aparam^2) \* vztest) // If x \<= aparam, apply load; otherwise 0  
68 : );  
69 :  
70 : // Assemble the system of equations  
71 : matrix A = avar(Vh, Vh, tgv = -1);  
72 : real[int] b = Lrhs(0, Vh);  
73 : set(A, solver = sparsesolver);  
74 :  
75 :  
76 : // Solving equations  
77 : real[int] sol = A^-1 \* b;  
78 : [ur, uz] = sol;  
79 :  
80 : // Write results to console  
81 : cout \<\< loadValue \<\< “,” \<\< E \<\< “,” \<\< mu \<\< “,” \<\< aparam \<\< “,” \<\< R; // Output parameters  
82 : for (int i = 0; i \< numR; ++i) {  
83 : cout \<\< “,” \<\< uz[i];  
84 : }  
85 : cout \<\< “\n”;  
86 : }

where as i am getting the Error as following 87 : sizestack + 1024 =3356 ( 2332 )

P,E,mu,a,R,uz0,uz1,uz2,uz3,uz4,uz5,uz6,uz7,uz8,uz9,uz10,uz11,uz12,uz13,uz14,uz15,uz16,uz17,uz18,uz19,uz20,uz21,uz22,uz23,uz24,uz25,uz26,uz27,uz28,uz29,uz30,uz31,uz32,uz33,uz34,uz35,uz36,uz37,uz38,uz39,uz40,uz41,uz42,uz43,uz44,uz45,uz46,uz47,uz48,uz49  
– Square mesh : nb vertices =2601 , nb triangles = 5000 , nb boundary edges 200 rmdup= 0

- 0 1\*0 1  
check ii.first =0 \< M 1  
check jj.first =0 \< N 1
- 0 1\*1 2  
check ii.first =1 \< M 1  
check jj.first =0 \< N 1
- 1 2\*0 1  
check ii.first =0 \< M 1  
check jj.first =1 \< N 1
- 1 2\*1 2  
check ii.first =1 \< M 1  
check jj.first =1 \< N 1
- 0 2\*0 2  
check ii.first =0 \< M 1  
check jj.first =0 \< N 1
- 0 2\*1 1  
check ii.first =1 \< M 1  
check jj.first =0 \< N 1
- 1 1\*0 2  
check ii.first =0 \< M 1  
check jj.first =1 \< N 1
- 1 1\*1 1  
check ii.first =1 \< M 1  
check jj.first =1 \< N 1  
value of N=1  
value of M=1  
current line = 71  
Exec error : Check BilinearOperator N M  
– number :1  
Exec error : Check BilinearOperator N M  
– number :1  
err code 8 , mpirank 0

---

<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:** [April 11, 2025, 4:10pm UTC](https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857/2 "2025-04-11T16:10:25Z")

</div>

You need to post the true source, not the output of the preprocessor, so that we can try it!

---

<div class="post-metadata">

**Author:** ![prayuktha](https://avatars.discourse-cdn.com/v4/letter/p/2acd7d/32.png) [@prayuktha](https://community.freefem.org/u/prayuktha)\
**Post date:** [April 12, 2025, 11:31am UTC](https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857/3 "2025-04-12T11:31:32Z")

</div>

I apologise for the confusion before, Please see the attached code for your reference. Thanks in advance.

```
 load "msh3" 

  real[int] Pvals = [5e4, 1e5, 2e5];     
  real[int] avals = [0.005, 0.01, 0.015];  
  real[int] Evals = [1e9, 7e10, 2e11];     
  real[int] muvals = [0.25, 0.3, 0.35];    
  real[int] Rvals = [0.01, 0.03, 0.05];    

int numR = 50; 
cout << "P,E,mu,a,R";
for (int i = 0; i < numR; ++i)
cout << ",uz" + i;
cout << "\n";

for (int ip = 0; ip < Pvals.n; ++ip)
for (int ia = 0; ia < avals.n; ++ia)
for (int ie = 0; ie < Evals.n; ++ie)
for (int imu = 0; imu < muvals.n; ++imu)
for (int ir = 0; ir < Rvals.n; ++ir)
{
real loadValue = Pvals[ip];     
real aparam = avals[ia];    
real E = Evals[ie];
real mu = muvals[imu];
real R = Rvals[ir];

real L = 2 * R; // domain width
real H = R; // domain height

// Mesh
int nx = 50, ny = 50;
mesh Th = square(nx, ny, [x*L, y*H]);

// FEM space
fespace Vh(Th, P1);
Vh ur, uz, vr, vz;

// Macros for strain components
macro epsrr(ur, uz) dx(ur) //
macro epszz(ur, uz) dy(uz) //
macro epsrz(ur, uz) (dy(ur) + dx(uz)) / 2 //
macro divu(ur, uz) (dx(ur) + dy(uz)) // total divergence

// Lame parameters
real lambda = E * mu / ((1. + mu) * (1. - 2. * mu));
real muL = E / (2. * (1. + mu));

// Variational form (bilinear form)
varf avar([ur, uz], [vr, vz]) =
    int2d(Th)(
        lambda * divu(ur, uz) * divu(vr, vz)
      + 2. * muL * (
            epsrr(ur, uz) * epsrr(vr, vz)
          + epszz(ur, uz) * epszz(vr, vz)
          + 2. * epsrz(ur, uz) * epsrz(vr, vz)
        )
    )
    + on(1, ur = 0, uz = 0); // Applying zero displacement on boundary 1

varf Lrhs([vr, vz], [vrtest, vztest]) =
    int1d(Th, 4)(
        (x <= aparam) * (-loadValue / (pi * aparam^2) * vztest) // If x <= aparam, apply load; otherwise 0
    );

matrix A = avar(Vh, Vh, tgv = -1); // System matrix
real[int] b = Lrhs(0, Vh); // Right-hand side vector
set(A, solver = sparsesolver);       
 real[int] sol = A^-1 * b;               
  [ur[], uz[]] = sol;                  
 cout << loadValue << "," << E << "," << mu << "," << aparam << "," << R;  
for (int i = 0; i < numR; ++i) {
  cout << "," << uz[][i];      
}
   cout << "\n";
}

```

---

<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:** [April 12, 2025, 2:16pm UTC](https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857/4 "2025-04-12T14:16:55Z")

</div>

The problem is that the space `Vh` you give in `matrix A = avar(Vh, Vh, tgv = -1)`  
has not the right size, since `varf avar([ur, uz], [vr, vz])` has two unknowns `ur,uz`.

You have to set

```auto
fespace Vh2(Th,[P1,P1]);
Vh2 [u1,u2];

```

The space `Vh2` has the right size. Then you correct your system resolution by

```auto
matrix A = avar(Vh2, Vh2, tgv = -1); // System matrix
real[int] b = Lrhs(0, Vh2); // Right-hand side vector
set(A, solver = sparsesolver);       
 real[int] sol = A^-1 * b;  
 u1[]=sol;// distribute the solution to [u1,u2], that belongs to Vh2
 ur=u1;
 uz=u2;

```

---

<div class="post-metadata">

**Author:** ![prayuktha](https://avatars.discourse-cdn.com/v4/letter/p/2acd7d/32.png) [@prayuktha](https://community.freefem.org/u/prayuktha)\
**Post date:** [April 12, 2025, 2:40pm UTC](https://community.freefem.org/t/traction-forces-bilinearoperator-nm-error/3857/5 "2025-04-12T14:40:42Z")

</div>

Thank you for your feedback! That the issue arises from using the incorrect finite element space `Vh` for the matrix assembly. I did the necessary changes. My code worked. I appreciate your guidance on this