Some testing...

Analyse/PlotVoltage.m

@@ -0,0 +1,21 @@
+%close all;
+clear all;
+clc
+
+tmp = load('../tmp/u1');
+
+t = tmp(:,1);
+u = tmp(:,2);
+
+L=numel(t);
+
+subplot(2,1,1);
+plot(t,u);
+
+dt=t(2)-t(1);
+
+f = (1:L)/L/dt;
+fu = fft(u)/L;
+subplot(2,1,2);
+plot(f(1:L/2),abs(fu(1:L/2)));
+

FDTD/engine.cpp

@@ -36,11 +36,11 @@ bool Engine::IterateTS(unsigned int iterTS)
 	for (unsigned int iter=0;iter<iterTS;++iter)
 	{
 		//voltage updates
-		for (pos[2]=1;pos[2]<Op->numLines[2]-1;++pos[2])
+		for (pos[0]=1;pos[0]<Op->numLines[0]-1;++pos[0])
 		{
 			for (pos[1]=1;pos[1]<Op->numLines[1]-1;++pos[1])
 			{
-				for (pos[0]=1;pos[0]<Op->numLines[0]-1;++pos[0])
+				for (pos[2]=1;pos[2]<Op->numLines[2]-1;++pos[2])
 				{
 					//do the updates here
 					//for x
@@ -78,15 +78,15 @@ bool Engine::IterateTS(unsigned int iterTS)
 					//do the updates here
 					//for x
 					curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]];
-					curr[0][pos[0]][pos[1]][pos[2]] -= Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]+1][pos[2]] - volt[2][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]+1] + volt[1][pos[0]][pos[1]][pos[2]]);
+					curr[0][pos[0]][pos[1]][pos[2]] += Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]][pos[1]+1][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] + volt[1][pos[0]][pos[1]][pos[2]+1]);
 
 					//for y
 					curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]];
-					curr[1][pos[0]][pos[1]][pos[2]] -= Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]+1] - volt[0][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]+1][pos[1]][pos[2]] + volt[2][pos[0]][pos[1]][pos[2]]);
+					curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]);
 
 					//for x
 					curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]];
-					curr[2][pos[0]][pos[1]][pos[2]] -= Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]+1][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]+1][pos[2]] + volt[0][pos[0]][pos[1]][pos[2]]);
+					curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]);
 				}
 			}
 		}

FDTD/processvoltage.cpp

@@ -38,5 +38,5 @@ void ProcessVoltage::Process()
 	}
 //	cerr << voltage << endl;
 	voltages.push_back(voltage);
-	file << Eng->numTS << "\t" << voltage << endl;
+	file << (double)Eng->numTS*Op->GetTimestep() << "\t" << voltage << endl;
 }

main.cpp

@@ -12,6 +12,7 @@
 using namespace std;
 
 void BuildMSL(ContinuousStructure &CSX);
+void BuildDipol(ContinuousStructure &CSX);
 
 int main(int argc, char *argv[])
 {
@@ -19,8 +20,7 @@ int main(int argc, char *argv[])
 
 	//*************** setup/read geometry ************//
 	ContinuousStructure CSX;
-//	CSX.ReadFromXML("csx-files/MSL.xml");
+	BuildDipol(CSX);
-	BuildMSL(CSX);
 
 	//*************** setup operator ************//
 	CartOperator cop;
@@ -38,7 +38,7 @@ int main(int argc, char *argv[])
 //	cop.DumpOperator2File("tmp/Operator");
 
 	cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl;
-	unsigned int NrIter = cop.GetNyquistNum(fmax);
+	unsigned int NrIter = cop.GetNyquistNum(fmax)/3;
 
 	cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl;
 
@@ -50,11 +50,13 @@ int main(int argc, char *argv[])
 	//*************** setup processing ************//
 	ProcessVoltage PV(&cop,&eng);
 	PV.OpenFile("tmp/u1");
-	double start[]={0,0,0};
+	double start[]={-500,-150,0};
-	double stop[]={0,50,0};
+	double stop[]={-500,150,0};
 	PV.DefineStartStopCoord(start,stop);
 	unsigned int maxIter = 5000;
 
+	PV.Process();
+//	NrIter=200;
 	//*************** simulate ************//
 	for (unsigned int i=0;i<maxIter;i+=NrIter)
 	{
@@ -73,6 +75,44 @@ int main(int argc, char *argv[])
     return 0;
 }
 
+void BuildDipol(ContinuousStructure &CSX)
+{
+	CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
+	mat->SetKappa(0.001);
+	CSX.AddProperty(mat);
+
+	CSPrimBox* matbox = new CSPrimBox(CSX.GetParameterSet(),mat);
+	matbox->SetCoord(0,-1000.0);matbox->SetCoord(1,1000.0);
+	matbox->SetCoord(2,-1000.0);matbox->SetCoord(3,1000.0);
+	matbox->SetCoord(4,-1000.0);matbox->SetCoord(5,1000.0);
+	CSX.AddPrimitive(matbox);
+
+	CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet());
+	elec->SetExcitation(1,1);
+	elec->SetExcitType(0);
+//	elec->SetDelay(2.0e-9);
+	CSX.AddProperty(elec);
+
+	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec);
+	box->SetCoord(0,-1.0);box->SetCoord(1,1.0);
+	box->SetCoord(2,-75.0);box->SetCoord(3,75.0);
+	box->SetCoord(4,-1.0);box->SetCoord(5,1.0);
+	CSX.AddPrimitive(box);
+
+	CSRectGrid* grid = CSX.GetGrid();
+
+	for (int n=-1000;n<=1000;n+=20)
+		grid->AddDiscLine(2,(double)n);
+	for (int n=-1000;n<=1000;n+=20)
+		grid->AddDiscLine(0,(double)n);
+	for (int n=-1000;n<=1000;n+=20)
+		grid->AddDiscLine(1,(double)n);
+
+	grid->SetDeltaUnit(1e-3);
+
+	CSX.Write2XML("tmp/Dipol.xml");
+}
+
 void BuildMSL(ContinuousStructure &CSX)
 {
 //	CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
@@ -109,14 +149,14 @@ void BuildMSL(ContinuousStructure &CSX)
 
 	CSRectGrid* grid = CSX.GetGrid();
 
-	for (int n=-100;n<=100;n+=10)
+	for (int n=-300;n<=300;n+=20)
-		grid->AddDiscLine(2,(double)n);
-	for (int n=-100;n<=100;n+=10)
 		grid->AddDiscLine(0,(double)n);
-	for (int n=0;n<=100;n+=10)
+	for (int n=0;n<=300;n+=20)
 		grid->AddDiscLine(1,(double)n);
+	for (int n=-1000;n<=1000;n+=20)
+		grid->AddDiscLine(2,(double)n);
 
 	grid->SetDeltaUnit(1e-3);
 
-	CSX.Write2XML("csx-files/MSL.xml");
+	CSX.Write2XML("tmp/MSL.xml");
 }