# Error in the solve routine for dynamic loading of an axisymmetric disc

**URL:** <https://community.freefem.org/t/error-in-the-solve-routine-for-dynamic-loading-of-an-axisymmetric-disc/4341>\
**Category:** General Discussion\
**Created:** [August 15, 2026, 12:14pm UTC](https://community.freefem.org/t/error-in-the-solve-routine-for-dynamic-loading-of-an-axisymmetric-disc/4341 "2026-08-15T12:14:17Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![EChase](https://avatars.discourse-cdn.com/v4/letter/e/b9e5f3/32.png) [@EChase](https://community.freefem.org/u/EChase)\
**Post date:** [August 15, 2026, 12:14pm UTC](https://community.freefem.org/t/error-in-the-solve-routine-for-dynamic-loading-of-an-axisymmetric-disc/4341/1 "2026-08-15T12:14:17Z")

</div>

Good morning

I’m trying to build a simple model of a disc, constrained at the center, with a dynamic axial load applied at its radius. Ultimately, I’d like to use the displacement at the outer diameter as a chekc on resonant frequency.

However, I am stuck with an error I do not understand with the solve function.

Any advice would be greatly appreciated.

THANK YOU!

here is the code:

/\* Forced vibration of a disc to look at dynamic transmissibility.  
Axisymmetric modeling to try to identify resonance to match with predictions  
\*/

/\*

* * *

------ Frequency Definitions ------

* * *

\*/

real[int] ff=[10.0,280.,550.,650.,800.]; // the table of selected frequencies  
// include “freq\_disc.edp”;

/\*

* * *

------ Geometry -------

* * *

\*/

real a=4.75_2.54/100; // disc radius  
real l=0.1_2.54/100; // disc thickness

/\*

* * *

------ Output Filename -------

* * *

\*/

string ftitle = “Point\_transmissibility.dat”;

/\*

* * *

------ Material Definitions and Loading-------

* * *

\*/

complex V0=0.0, V1=1.0; // electrode potentials  
real rho = 7850; // material density

complex E = 210E9; // Young’s modulus steel  
complex nu = 1/3; // Poisson ratio steel

// Stiff - Stiffness Matrix  
func Stiff = (E/((1+nu)_(1-2_nu)))\* [[1-nu, nu, nu, 0],  
[nu, 1-nu, nu, 0],  
[nu, nu, 1-nu, 0],  
[0 , 0 , 0 , 1-2\*nu]];

// force applied at the end of the disc  
complex force = 1.0;  
/\*

* * *

------ Mesh Definition and Creation -------

* * *

\*/

int MM=10, NN=5; // grid resolution  
mesh Sp = square(MM, NN, [a_x, l_y]); // grid generation  
plot(Sp, ps=“ff\_mesh.eps”); // plot and save mesh

//\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*  
// ------ Identify Gridpoint Coordinates --------

// Vh xx=x, yy=y;  
real[int] xx(Sp.nv);  
real[int] yy(Sp.nv);  
for(int i = 0; i \<Sp.nv; i++){  
xx[i] = Sp(i).x;  
yy[i] = Sp(i).y;  
}

/\*

* * *

------ Macros --------

* * *

_/  
real sqrt2 = sqrt(2.);  
//macro Strain(ur,uz) [dx(ur), ur/x, dy(uz), 0.5_(dy(ur)+dx(uz))] //EOM  
macro Strain(ur,uz) [dx(ur), ur/x, dy(uz), (dy(ur)+dx(uz))/sqrt2] //EOM

//\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*

// ------ Problem Definition -------  
fespace Uh(Sp,[P1,P1]); // piecewise linear FE  
fespace Vh(Sp,[P1]);  
Uh [ur,uz]; // variational variables  
Uh [vr,vz]; // variational variables

Vh xcoord=x;

// Iterate through all the frequencies

for(int jj=0; jj\<ff.n; jj++) { // for all frequencies  
real f0 = ff[jj], w0 = 2_pi_f0; cout \<\< f0/1e3 \<\< “kHz” \<\< endl;

solve AxiDisc([ur,uz],[vr,vz,]) = // variational equation!  
int2d(Sp)(0.5_x_rho_xcoord_w0^2\*[vr,vz]'_[ur,uz])  
-int2d(Sp)(x_Strain(vr,vz)'_Stiff_Strain(ur,uz)  
)  
-int1D(Uh, 2)(  
force\*uz  
)

- on(4, ur=0) // BC: axis of symmetry
- on(4, uz=0); // BC: axis of symmetry

Vh urreal=real(ur);  
Vh uzreal=real(uz);  
mesh Sp2 = movemesh (Sp,[x+1000.\*urreal, y+1000.\*uzreal]);  
plot(Sp,Sp2,urreal, wait=true, cmm=“f = " + f0/1000 + " kHz”, fill = true, value=true);  
// plot(Sp,Sp2,uzreal, wait=true, cmm=“f = " + f0/1000 + " kHz”, fill = true, value=true);

}

---

<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:** [August 16, 2026, 2:57pm UTC](https://community.freefem.org/t/error-in-the-solve-routine-for-dynamic-loading-of-an-axisymmetric-disc/4341/2 "2026-08-16T14:57:01Z")

</div>

There is an extra comma in `solve AxiDisc([ur,uz],[vr,vz,])` after “vz”.  
Moreover:  
– you cannot put complex numbers if the `fespace` is real,  
– the 1d integral should be `-int1d(Sp, 2)(force*vz)`

---

<div class="post-metadata">

**Author:** ![EChase](https://avatars.discourse-cdn.com/v4/letter/e/b9e5f3/32.png) [@EChase](https://community.freefem.org/u/EChase)\
**Post date:** [August 24, 2026, 3:47pm UTC](https://community.freefem.org/t/error-in-the-solve-routine-for-dynamic-loading-of-an-axisymmetric-disc/4341/3 "2026-08-24T15:47:10Z")

</div>

> [@fb77](#):
>
> – you cannot put complex numbers if the `fespace` is real,  
> – the 1d integral should be `-int1d(Sp, 2)(force*vz)`

Thank you immensely for the help!