Mesh dissapears

General Discussion
1

So for some reason my mesh vanishes after the first iteration. Any ideas why?
For some context I am trying to solve Maxwells equations using the Yee scheme for time discretisation:

My aim is to get a wave simulation, however I don’t seem to be getting very far, this is the first time my mesh has vanished after an iteration…

I suspect the issue is in the for loop at the bottom. Any input is much appreciated!





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

real w = 2*pi*10;
// Set initial conditions

real ExInitial = 0;
real EyInitial = 0;
real HxInitial = 0;
real HyInitial = 0;

// For mesh setup
int nn = 30;
int boxheight = 10;
int boxwidth = 10;
int Sourcex = 5;
int Sourcey = 5;
int source = 10;
int box = 1;
//constants
real mu0 = 1.25663706e-6;
real epsilon0 = 8.85e-12;
func epsilonr = 10*exp(-2*y);

macro Curl(Sx,Sy,Sz) ((dy(Sz) - dz(Sy))),((dz(Sx) - dx(Sz))),((dx(Sy) - dy(Sx))) //
macro dotproduct(Ax,Ay,Az,Bx,By,Bz) ((Ax*Bx) +  (Ay*By) + (Az*Bz))  //
macro dotproductcurl(Sx,Sy,Sz,Ax,Ay,Az) (
	((dy(Sz) - dz(Sy))) * Ax +
	((dz(Sx) - dx(Sz))) * Ay +
	((dx(Sy) - dy(Sx))) * Az
) //



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

// fespaces
fespace Vh(Th,P1);

Vh Exnext, Eynext, Eznext;
Vh Exold, Eyold, Ezold;

Vh Hxnext, Hynext, Hznext;
Vh Hxold, Hyold , Hzold;

Vh HIntensity;
Vh EIntensity;
Vh vx,vy, vz;

Vh Jz, vj;
fespace Xh(Th,[P1,P1,P1]); // For vector fields
// Set initial conditions

Exold = 0;
Eyold = 0;
Ezold = 0;



Hxold = 0;
Hyold = 0;
Hzold = 0;


// Want to get it in the form A[Exnext,Eynext,Eznext] = B[Exold,Eyold,Ezold] + C([Hxnext,Hynext,Hznext]) + bc[Exold,Eyold,Ezold] + bc[Hznext,Hynext,Hznext]
// D[Hx,Hy,Hz] = E[Hxold,Hyold,Hzold] + F([Exnext,Eynext,Eznext]) + bc[Hxold,Hyold,Hzold] + bc[Eznext,Eynext,Eznext]

// Efield!

varf aa([Exnext,Eynext,Eznext],[vx,vy,vz]) = int2d(Th)(
	dotproduct(Exnext,Eynext,Eznext,vx,vy,vz)
);

varf bb([Exold,Eyold,Ezold],[vx,vy,vz]) = int2d(Th)(
	dotproduct(Exold,Eyold,Ezold,vx,vy,vz)
);

varf cc([Hxnext,Hynext,Hznext],[vx,vy,vz]) = int2d(Th)(
	(dT/(epsilon0*epsilonr)) * dotproductcurl(Hxnext,Hynext,Hznext,vx,vy,vz)
);

varf bcB([Exold,Eyold,Ezold],[vx,vy,vz]) = on(Source,Eyold = sin(w*(t)));
varf bcC([Hxnext,Hynext,Hznext],[vx,vy,vz]) = on(Source,Hznext = sin(w*(t + 0.5*dT)));


// Hfield!

varf dd([Hxnext,Hynext,Hznext],[vx,vy,vz]) = int2d(Th)(
	dotproduct(Hxnext,Hynext,Hznext,vx,vy,vz)
);

varf ee([Hxold,Hyold,Hzold],[vx,vy,vz]) = int2d(Th)(
	dotproduct(Hxold,Hyold,Hzold,vx,vy,vz)
);

varf ff([Exold,Eyold,Ezold],[vx,vy,vz]) = int2d(Th)(
	-(dT/mu0)* dotproductcurl(Exold,Eyold,Ezold,vx,vy,vz)
);

varf bcE([Hxold,Hyold,Hzold],[vx,vy,vz]) = on(Source,Hyold = sin(w*(t-0.5*dT)));
varf bcF([Exold,Eyold,Ezold],[vx,vy,vz]) = on(Source,Eyold = sin(w*(t)));

// matrices

matrix A = aa(Xh,Xh,solver  = UMFPACK);
matrix B = bb(Xh,Xh);
matrix C = cc(Xh,Xh);
matrix D = dd(Xh,Xh,solver = UMFPACK);
matrix E = ee(Xh,Xh);
matrix F = ff(Xh,Xh);

// boundary condition vector(s)
real[int] bcCVec = bcC(0,Xh);
real[int] bcBVec = bcB(0,Xh);
real[int] bcEVec = bcE(0,Xh);
real[int] bcFVec = bcF(0,Xh);

// State vectors
real[int] EnextSol(Xh.ndof);
real[int] EoldSol(Xh.ndof);
real[int] HnextSol(Xh.ndof);
real[int] HoldSol(Xh.ndof);
real[int] Eaux(Xh.ndof);
real[int] Haux(Xh.ndof);


// Temporary values
Xh [ExOldtemp, EyOldtemp, EzOldtemp] = [Exold,Eyold,Ezold];
Xh [HxOldtemp, HyOldtemp, HzOldtemp] = [Hxold,Hyold,Hzold];
Xh [ExNexttemp, EyNexttemp, EzNexttemp] = [Exnext,Eynext,Eznext];
Xh [Hxnexttemp, Hynexttemp, Hznexttemp] = [Hxnext,Hynext,Hznext];
EoldSol = ExOldtemp[]; // Initialse state vectors
HoldSol = HxOldtemp[];
HnextSol = Hxnexttemp[];


// Set up loop

int m, M = T/dT;
//loop!

for (m = 0; m < M; m++) {
	// Initialise everything for appropriate time.
	bcCVec = bcC(0,Xh);
	bcBVec = bcB(0,Xh);
	bcEVec = bcE(0,Xh);
	bcFVec = bcF(0,Xh);

	// Set up the right hand side of the equation to calculate E_{t+dt}
	Eaux += B*EoldSol; Eaux += bcBVec; Eaux += C*HnextSol; Eaux += bcCVec;
	EnextSol = A^-1 * Eaux;

	// Exctract data from ENextSol
	Xh [ExNexttemp, EyNexttemp, EzNexttemp] = EnextSol;
	Exnext = ExNexttemp; Eynext = EyNexttemp; Eznext = EzNexttemp;

	//Now calculate for H field
	Haux += E*HoldSol; Haux += bcEVec; Haux += F*EoldSol; Haux += bcFVec;
	HnextSol = D^-1 * Haux;
	// Extract data from HnextSol
	Xh [Hxnexttemp, Hynexttemp, Hznexttemp] = HnextSol;
	Hxnext = Hxnexttemp; Hynext = Hznexttemp; Hznext = Hznexttemp;

	// Increment time
	t += dT;

	// Plot data;
	EIntensity = ExNexttemp^2 + EyNexttemp^2;
	plot(EIntensity,[ExNexttemp,EyNexttemp],wait = true, fill = true);
	cout << "t = " << t << "\n";

	

}