What causes a length interpole warning?

I have the following weak form for a vector wave equation

-\int_{\Omega} \nabla \cdot \vec{v} \nabla \cdot \vec{E}^m d\Omega + \int_{\Omega} (\nabla \times \vec{v})\cdot (\nabla \times \vec{E}^m) - \frac{1}{dt^2} \int_{\Omega} \mu_0 \epsilon_0 \epsilon_r \vec{E}^m \cdot \vec{v}d\Omega = -\frac{2}{dt^2} \int_{\Omega}mu_0 \epsilon_0 \epsilon_r \vec{E}^{m-1}\cdot\vec{v} d\Omega -\frac{1}{dt^2} \int_{\Omega}\mu0 \epsilon_0 \epsilon_r \vec{E}^{m-2} \cdot\vec{v} d\Omega

I am treating it as an equation of the form

\hat{A} [E_x^m,E_y^m]^T = \hat{B}[E_x^{m-1},E_y^{m-1}]^T + \hat{C}[E_x^{m-2},E_y^{m-2}]^T

When I solve this equation I get the following error:

Warning LengthInterpole: ( i = 235 l = 39.8156 sss 0.5 ) 0.1

And the resulting plot is blank. Any ideas why? All input is much appreciated! Here is my code:

real T = 100000;
real dT = 0.06;
real t = 0;
real k = 1;

real w = 2*pi*0.5;
// Set initial conditions
real ExInitialGradient = 1;
real EyInitialGradient = 1; 
real ExInitial = 0;
real EyInitial = 0;

// For mesh setup
int nn = 30;
int boxheight = 1000;
int boxwidth = 1000;
int Sourcex = 500;
int Sourcey = 500;
int source = 10;
int box = 1;
//constants
real mu0 = 1.25663706e-6;
real epsilon0 = 8.85e-12;



real a = 0.1; // radius of the circle;
// Set up mesh
border Source(t = 0,2*pi){x = a*cos(t)+Sourcex;y=a*sin(t)+Sourcey;label = source;}
border rightwall(t = 0,boxheight){x = boxwidth; y = t;label = box;}
border leftwall(t = 0,boxheight){x = 0;y = boxheight - t;label = box; }
border bottom(t = 0,boxwidth){x = t;y =0;label = box;}
border top(t = 0,boxwidth){x = boxwidth - t;y = boxheight;label = box;}
mesh Th = buildmesh(Source(10)+rightwall(10)+leftwall(10)+top(10)+bottom(10)); // Our mesh

fespace Vh(Th,P1);
Vh Ex,Ey,Exold,Eyold,Exoldold,Eyoldold,intensity,vx,vy; // p = charge density
func real p() { 
	if (Th(x,y).region == Source) {return sin(t);} else {return 0;}
	
}
Exold = 0;
Eyold = 0;
Exoldold = 0;
Eyoldold = 0;

func epsilonr = 1000*exp(-100*y); // relative permitivity as a function of height;


macro curl(S) (dx(S#y) - dy(S#x)) //
macro div(S) (dx(S#x) + dy(S#y)) // 
macro dotproduct(Sx,Sy,Zx,Zy) (Sx*Zx + Sy*Zy) //

//Next we set up are variable forms for A[Ex,Ey] = B[Exold,Eyold] + C[Exoldold,Eyoldold] + bcEold + bcEoldold;
varf aa([Ex,Ey],[vx,vy]) = int2d(Th)(
	-(div(E)) * div(v)
) + int2d(Th)(
	curl(E) * curl(v)
) - int2d(Th)(
	(1/dT^2)*mu0*epsilon0*epsilonr*dotproduct(Ex,Ey,vx,vy)
) + on(Source, Ey = sin(4*x)) + on(box, Ex = 0, Ey = 0);

varf bb([Exold,Eyold],[vx,vy]) = int2d(Th)(
	-2*(1/dT^2)*mu0*epsilon0*epsilonr*dotproduct(Ex,Ey,vx,vy)
) +on(box, Exold = 0, Eyold = 0);

varf cc([Exoldold,Eyoldold],[vx,vy]) = int2d(Th)(
	-(1/dT^2)*mu0*epsilon0*epsilonr*dotproduct(Exoldold,Eyoldold,vx,vy)
) +on(box, Exoldold = 0, Eyoldold = 0);

// Set up A, B, and C
// First need a vector fespace
fespace Xh(Th,[P1,P1,P1]); // [Ex, Ey, p]
fespace XhWithoutP(Th,[P1,P1]);
matrix A = aa(XhWithoutP,XhWithoutP,solver = UMFPACK);
matrix B = bb(XhWithoutP,XhWithoutP);
matrix C = cc(XhWithoutP,XhWithoutP);

// Next need to set up vectors that will hold the state of E, Eold, and Eoldold;
real[int] sol(XhWithoutP.ndof), solOld(XhWithoutP.ndof), solOldOld(XhWithoutP.ndof), aux(XhWithoutP.ndof);
XhWithoutP [ExOldTemp,EyOldTemp] = [Exold,Eyold];
XhWithoutP [ExOldOldTemp,EyOldOldTemp] = [Exoldold,Eyoldold];
XhWithoutP [ExNow, EyNow];

int m, M = T/dT;
solOld = ExOldTemp[];
solOldOld = ExOldOldTemp[];

for (m = 0; m < M; m++){
	aux = B*solOld; aux += C*solOldOld;
	sol = A^-1 * aux;
	ExNow[] = sol;
	Ex = ExNow; Ey = EyNow; 
	Exoldold = Exold;
	Eyoldold = Eyold;
	Exold = Ex;
	Eyold = Ey;
	XhWithoutP [ExTemp,EyTemp] = [Ex,Ey];
	solOld =Ex[];
	XhWithoutP [ExOldTemp,EyoldTemp] = [Exold,Eyold];
	solOldOld = ExOldOldTemp[];
	
	t += dT;

	A = aa(XhWithoutP,XhWithoutP,solver = UMFPACK);

	intensity  = Ex^2 + Ey^2;
	plot([Ex,Ey], intensity, wait = true, fill = true);


}

Just a remark, the warning is not important , the mesh is correct.

but E[x,y]old and E[x,y]oldold are 0 so sol is 0 and Ex, EY are also 0.
Then you see 0 and any time !!!

1 Like

Thanks for the feedback! One thing I am noticing is that despite specifying the value of the E field on Source it is not changing

Would you have any advice on how to update the values of Ex and Ey on source in a loop?

Be careful in the copy order !!!

1 Like