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