Traction Forces - BilinearOperator NM error

prayuktha · April 11, 2025, 3:31pm

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, [xL, yH]);
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

fb77 · April 11, 2025, 4:10pm

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

prayuktha · April 12, 2025, 11:31am

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";
}

fb77 · April 12, 2025, 2:16pm

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

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

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

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;

prayuktha · April 12, 2025, 2:40pm

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

Topic		Replies	Views
Operator Error in problem formulation General Discussion	3	35	April 14, 2025
Elasticity neo-Hookean material General Discussion	0	252	February 17, 2022
Numerical results error General Discussion	8	170	September 7, 2023
Elasticity problem General Discussion	11	270	September 8, 2023
Solve U=K^-1.F using stiffness matrix for vector fields General Discussion	12	1350	September 19, 2021

Traction Forces - BilinearOperator NM error

Related topics