Engine fix!, processing class, lots of fixes

.gitignore

@@ -2,6 +2,7 @@ release/
 debug/
 doc/
 obj/
+tmp*
 Makefile*
 *~
 *.so*

FDTD/cartoperator.cpp

@@ -1,7 +1,7 @@
 #include "cartoperator.h"
 #include "tools/array_ops.h"
 
-CartOperator::CartOperator()
+CartOperator::CartOperator() : Operator()
 {
     Init();
 }
@@ -325,6 +325,7 @@ bool CartOperator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
 	inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
 	inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
 
+	return true;
 }
 
 

FDTD/engine.cpp

@@ -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->vi[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]+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]]);
 
 					//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->vi[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]+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]]);
 
 					//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->vi[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]+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]]);
 				}
 			}
 		}

FDTD/engine.h

@@ -5,6 +5,7 @@
 
 class Engine
 {
+	friend class ProcessVoltage;
 public:
 	Engine(Operator* op);
 	virtual ~Engine();
@@ -15,6 +16,8 @@ public:
 	//!Iterate a number of timesteps
 	virtual bool IterateTS(unsigned int iterTS);
 
+	unsigned int GetNumberOfTimesteps() {return numTS;};
+
 protected:
 	Operator* Op;
 

FDTD/operator.cpp

@@ -1,3 +1,4 @@
+#include <fstream>
 #include "operator.h"
 #include "tools/array_ops.h"
 
@@ -46,6 +47,37 @@ void Operator::Reset()
 	Operator::Init();
 }
 
+unsigned int Operator::GetNyquistNum(double fmax)
+{
+	if (dT==0) return 1;
+	double T0 = 1/fmax;
+	return floor(T0/2/dT);
+}
+
+bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord)
+{
+	bool ok=true;
+	for (int n=0;n<3;++n)
+	{
+		if (dcoord[n]<discLines[n][0]) {ok=false;uicoord[n]=0;};
+		if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1;};
+		for (unsigned int i=0;i<numLines[n]-1;++i)
+		{
+			if (dcoord[n]<=discLines[n][i])
+			{
+				if (fabs(dcoord[n]-discLines[n][i])<(fabs(dcoord[n]-discLines[n][i+1])))
+					uicoord[n]=i;
+				else
+					uicoord[n]=i+1;
+				i = numLines[n];
+			}
+		}
+	}
+//	cerr << "Operator::SnapToMesh: " << discLines[0][uicoord[0]] << " " << discLines[1][uicoord[1]] << " " << discLines[2][uicoord[2]] << endl;
+//	cerr << "Operator::SnapToMesh: " << uicoord[0] << " " << uicoord[1] << " " << uicoord[2] << endl;
+	return ok;
+}
+
 void Operator::SetGeometryCSX(ContinuousStructure* geo)
 {
 	if (geo==NULL) return;
@@ -87,3 +119,31 @@ void Operator::CalcGaussianPulsExcitation(double f0, double fc)
 //		cerr << ExciteSignal[n] << endl;
 	}
 }
+
+void Operator::DumpOperator2File(string filename)
+{
+	ofstream file(filename.c_str(),ios_base::out);
+//	file.open;
+	if (file.is_open()==false)
+	{
+		cerr << "Operator::DumpOperator2File: Can't open file: " << filename << endl;
+		return;
+	}
+	file << "########### Operator vv ###########" << endl;
+	file << "ix\tiy\tiz\tvv_x\tvv_y\tvv_z" << endl;
+	Dump_N_3DArray2File(file,vv,numLines);
+
+	file << "########### Operator vi ###########" << endl;
+	file << "ix\tiy\tiz\tvi_x\tvi_y\tvi_z" << endl;
+	Dump_N_3DArray2File(file,vi,numLines);
+
+	file << "########### Operator iv ###########" << endl;
+	file << "ix\tiy\tiz\tiv_x\tiv_y\tiv_z" << endl;
+	Dump_N_3DArray2File(file,iv,numLines);
+
+	file << "########### Operator ii ###########" << endl;
+	file << "ix\tiy\tiz\tii_x\tii_y\tii_z" << endl;
+	Dump_N_3DArray2File(file,ii,numLines);
+
+	file.close();
+}

FDTD/operator.h

@@ -17,7 +17,7 @@ public:
 
 	virtual void SetGeometryCSX(ContinuousStructure* geo);
 
-	virtual int CalcECOperator() {};
+	virtual int CalcECOperator() {return -1;};
 
 	virtual void CalcGaussianPulsExcitation(double f0, double fc);
 
@@ -25,12 +25,17 @@ public:
 	virtual void ApplyMagneticBC(bool* dirs) {};
 
 	double GetTimestep() {return dT;};
+	unsigned int GetNyquistNum(double fmax);
 	double GetNumberCells();
 
 	void ShowSize();
 
+	void DumpOperator2File(string filename);
+
 	virtual void Reset();
 
+	bool SnapToMesh(double* coord, unsigned int* uicoord);
+
 protected:
 	virtual void Init();
 	ContinuousStructure* CSX;
@@ -51,14 +56,14 @@ protected:
 
 	//E-Field Excitation
 	//!	  Calc the electric field excitation.
-	virtual bool CalcEFieldExcitation() {};
+	virtual bool CalcEFieldExcitation() {return false;};
 	unsigned int E_Ex_Count;
 	unsigned int* E_Ex_index[3];
 	FDTD_FLOAT* E_Ex_amp[3]; //represented as edge-voltages!!
 	unsigned int* E_Ex_delay;
 
 	//Calc timestep only internal use
-	virtual double CalcTimestep() {};
+	virtual double CalcTimestep() {return 0;};
 	double dT; //FDTD timestep!
 };
 

FDTD/processing.cpp

@@ -0,0 +1,24 @@
+#include "processing.h"
+
+Processing::Processing(Operator* op, Engine* eng)
+{
+	Op=op;
+	Eng=eng;
+}
+
+Processing::~Processing()
+{
+	file.close();
+}
+
+void Processing::OpenFile(string outfile)
+{
+	if (file.is_open()) file.close();
+
+	file.open(outfile.c_str());
+	if (file.is_open()==false)
+	{
+		cerr << "Can't open file: " << outfile << endl;
+		return;
+	}
+}

FDTD/processing.h

@@ -0,0 +1,25 @@
+#ifndef PROCESSING_H
+#define PROCESSING_H
+
+#include <iostream>
+#include <fstream>
+#include "operator.h"
+#include "engine.h"
+
+class Processing
+{
+public:
+	Processing(Operator* op, Engine* eng);
+	virtual ~Processing();
+
+	virtual void OpenFile(string outfile);
+
+	virtual void Process() {};
+
+protected:
+	Operator* Op;
+	Engine* Eng;
+	ofstream file;
+};
+
+#endif // PROCESSING_H

FDTD/processvoltage.cpp

@@ -0,0 +1,42 @@
+#include "processvoltage.h"
+
+ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng)
+{
+}
+
+ProcessVoltage::~ProcessVoltage()
+{
+}
+
+void ProcessVoltage::DefineStartStopCoord(double* dstart, double* dstop)
+{
+	if (Op->SnapToMesh(dstart,start)==false) cerr << "ProcessVoltage::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
+	if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessVoltage::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
+}
+
+void ProcessVoltage::Process()
+{
+	FDTD_FLOAT voltage=0;
+	unsigned int pos[3]={start[0],start[1],start[2]};
+//	cerr << Eng->volt[1][pos[0]][pos[1]][pos[2]] << endl;
+	for (int n=0;n<3;++n)
+	{
+		if (start[n]<stop[n])
+		{
+			for (;pos[n]<stop[n];++pos[n])
+			{
+				voltage+=Eng->volt[n][pos[0]][pos[1]][pos[2]];
+//				cerr << n << " " << pos[0] << " " << pos[1] << " " <<  pos[2] << " " << Eng->volt[n][pos[0]][pos[1]][pos[2]] << endl;
+			}
+
+		}
+		else
+		{
+			for (;pos[n]>stop[n];--pos[n])
+				voltage-=Eng->volt[n][pos[0]][pos[1]][pos[2]];
+		}
+	}
+//	cerr << voltage << endl;
+	voltages.push_back(voltage);
+	file << Eng->numTS << "\t" << voltage << endl;
+}

FDTD/processvoltage.h

@@ -0,0 +1,23 @@
+#ifndef PROCESSVOLTAGE_H
+#define PROCESSVOLTAGE_H
+
+#include "processing.h"
+
+class ProcessVoltage : public Processing
+{
+public:
+	ProcessVoltage(Operator* op, Engine* eng);
+	virtual ~ProcessVoltage();
+
+	void DefineStartStopCoord(double* dstart, double* dstop);
+
+	virtual void Process();
+
+protected:
+	unsigned int start[3];
+	unsigned int stop[3];
+
+	vector<FDTD_FLOAT> voltages;
+};
+
+#endif // PROCESSVOLTAGE_H

main.cpp

@@ -1,7 +1,11 @@
 #include <iostream>
+#include <fstream>
+#include <sstream>
 #include <time.h>
+#include "tools/array_ops.h"
 #include "FDTD/cartoperator.h"
 #include "FDTD/engine.h"
+#include "FDTD/processvoltage.h"
 
 #include "ContinuousStructure.h"
 
@@ -11,65 +15,85 @@ void BuildMSL(ContinuousStructure &CSX);
 
 int main(int argc, char *argv[])
 {
-	ContinuousStructure CSX;
-
 	time_t startTime=time(NULL);
 
+	//*************** setup/read geometry ************//
+	ContinuousStructure CSX;
 //	CSX.ReadFromXML("csx-files/MSL.xml");
 	BuildMSL(CSX);
 
+	//*************** setup operator ************//
 	CartOperator cop;
 	cop.SetGeometryCSX(&CSX);
 
 	cop.CalcECOperator();
 
-	cop.CalcGaussianPulsExcitation(1e9,1e9);
+	double fmax=1e9;
+	cop.CalcGaussianPulsExcitation(fmax/2,fmax/2);
 
 	time_t OpDoneTime=time(NULL);
 
 	cop.ShowSize();
 
+//	cop.DumpOperator2File("tmp/Operator");
+
+	cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl;
+	unsigned int NrIter = cop.GetNyquistNum(fmax);
+
 	cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl;
 
+	//create FDTD engine
 	Engine eng(&cop);
 
 	time_t currTime = time(NULL);
 
-	unsigned int NrIter = 500;
+	//*************** setup processing ************//
+	ProcessVoltage PV(&cop,&eng);
+	PV.OpenFile("tmp/u1");
+	double start[]={0,0,0};
+	double stop[]={0,50,0};
+	PV.DefineStartStopCoord(start,stop);
+	unsigned int maxIter = 5000;
 
+	//*************** simulate ************//
+	for (unsigned int i=0;i<maxIter;i+=NrIter)
+	{
 		eng.IterateTS(NrIter);
+		PV.Process();
+	}
 
+	//*************** postproc ************//
 	time_t prevTime = currTime;
 	currTime = time(NULL);
 
 	double t_diff = difftime(currTime,prevTime);
 
-	cerr << "Time for " << NrIter << " iterations with " << cop.GetNumberCells() << " cells : " << t_diff << " sec" << endl;
+	cerr << "Time for " << eng.GetNumberOfTimesteps() << " iterations with " << cop.GetNumberCells() << " cells : " << t_diff << " sec" << endl;
-	cerr << "Speed (MCells/s): " << (double)cop.GetNumberCells()*(double)NrIter/t_diff/1e6 << endl;
+	cerr << "Speed (MCells/s): " << (double)cop.GetNumberCells()*(double)eng.GetNumberOfTimesteps()/t_diff/1e6 << endl;
     return 0;
 }
 
 void BuildMSL(ContinuousStructure &CSX)
 {
-	CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
+//	CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
-	mat->SetEpsilon(3.6);
+//	mat->SetEpsilon(3.6);
-	CSX.AddProperty(mat);
+//	CSX.AddProperty(mat);
-
+//
-	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat);
+//	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat);
-	box->SetCoord(0,-200.0);box->SetCoord(1,200.0);
+//	box->SetCoord(0,-200.0);box->SetCoord(1,200.0);
-	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
+//	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
-	box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0);
+//	box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0);
-	CSX.AddPrimitive(box);
+//	CSX.AddPrimitive(box);
-
+//
-	CSPropMaterial* MSL = new CSPropMaterial(CSX.GetParameterSet());
+//	CSPropMaterial* MSL = new CSPropMaterial(CSX.GetParameterSet());
-	MSL->SetKappa(56e6);
+//	MSL->SetKappa(56e6);
-	CSX.AddProperty(MSL);
+//	CSX.AddProperty(MSL);
-
+//
-	box = new CSPrimBox(CSX.GetParameterSet(),MSL);
+//	box = new CSPrimBox(CSX.GetParameterSet(),MSL);
-	box->SetCoord(0,-20.0);box->SetCoord(1,20.0);
+//	box->SetCoord(0,-20.0);box->SetCoord(1,20.0);
-	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
+//	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
-	box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0);
+//	box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0);
-	CSX.AddPrimitive(box);
+//	CSX.AddPrimitive(box);
 
 	CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet());
 	elec->SetExcitation(1,1);
@@ -77,19 +101,19 @@ void BuildMSL(ContinuousStructure &CSX)
 //	elec->SetDelay(2.0e-9);
 	CSX.AddProperty(elec);
 
-	box = new CSPrimBox(CSX.GetParameterSet(),elec);
+	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec);
-	box->SetCoord(0,-20.0);box->SetCoord(1,20.0);
+	box->SetCoord(0,-1.0);box->SetCoord(1,1.0);
 	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
-	box->SetCoord(4,0.0);box->SetCoord(5,0.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)
+	for (int n=-100;n<=100;n+=10)
 		grid->AddDiscLine(2,(double)n);
-	for (int n=-200;n<=200;n+=10)
+	for (int n=-100;n<=100;n+=10)
 		grid->AddDiscLine(0,(double)n);
-	for (int n=0;n<=300;n+=10)
+	for (int n=0;n<=100;n+=10)
 		grid->AddDiscLine(1,(double)n);
 
 	grid->SetDeltaUnit(1e-3);

openEMS.pro

@@ -20,11 +20,15 @@ SOURCES += main.cpp \
     tools/AdrOp.cpp \
     FDTD/engine.cpp \
     FDTD/operator.cpp \
-    tools/array_ops.cpp
+    tools/array_ops.cpp \
+    FDTD/processvoltage.cpp \
+    FDTD/processing.cpp
 HEADERS += FDTD/cartoperator.h \
     tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
     FDTD/engine.h \
     FDTD/operator.h \
-    tools/array_ops.h
+    tools/array_ops.h \
+    FDTD/processvoltage.h \
+    FDTD/processing.h

tools/array_ops.cpp

@@ -1,4 +1,5 @@
 #include "array_ops.h"
+#include <ostream>
 
 FDTD_FLOAT*** Create3DArray(unsigned int* numLines)
 {
@@ -56,3 +57,21 @@ void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines)
 	}
 	delete[] array;
 }
+
+void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, unsigned int* numLines)
+{
+	unsigned int pos[3];
+	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+	{
+		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+		{
+			for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+			{
+				file << pos[0] << "\t" << pos[1] << "\t" << pos[2];
+				for (int n=0;n<3;++n)
+					file << "\t" << array[n][pos[0]][pos[1]][pos[2]];
+				file << endl;
+			}
+		}
+	}
+}

tools/array_ops.h

@@ -9,4 +9,6 @@ void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines);
 FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines);
 void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines);
 
+void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, unsigned int* numLines);
+
 #endif // ARRAY_OPS_H