Field dumps (E-field) TD, hard E-source

Analyse/PlotVoltage.m

@@ -11,6 +11,7 @@ L=numel(t);
 
 subplot(2,1,1);
 plot(t,u);
+grid on;
 
 dt=t(2)-t(1);
 
@@ -18,4 +19,5 @@ f = (1:L)/L/dt;
 fu = fft(u)/L;
 subplot(2,1,2);
 plot(f(1:L/2),abs(fu(1:L/2)));
+grid on;
 

FDTD/cartoperator.cpp

@@ -432,14 +432,22 @@ bool CartOperator::CalcEFieldExcitation()
 				if (prop)
 				{
 					CSPropElectrode* elec = prop->ToElectrode();
-					if (elec->GetExcitType()==0)
+					if ((elec->GetExcitType()==0) || (elec->GetExcitType()==1)) //soft or hard E-Field excite!
 					{
 						vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
 						for (int n=0;n<3;++n)
 						{
 							vIndex[n].push_back(pos[n]);
 							double delta=MainOp->GetIndexDelta(n,pos[n])*gridDelta;
-							vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta);
+							if (elec->GetActiveDir(n))
+								vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta);
+							else
+								vExcit[n].push_back(0);
+							if ((elec->GetExcitType()==1) && (elec->GetActiveDir(n))) //hard excite
+							{
+								vv[n][pos[0]][pos[1]][pos[2]] = 0;
+								vi[n][pos[0]][pos[1]][pos[2]] = 0;
+							}
 						}
 					}
 				}

FDTD/engine.cpp

@@ -62,7 +62,8 @@ bool Engine::IterateTS(unsigned int iterTS)
 		for (unsigned int n=0;n<Op->E_Ex_Count;++n)
 		{
 			exc_pos = (int)numTS - (int)Op->E_Ex_delay[n];
-			exc_pos*= (exc_pos>0 && exc_pos<Op->ExciteLength);
+			exc_pos*= (exc_pos>0 && exc_pos<(int)Op->ExciteLength);
+//			if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
 			volt[0][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[0][n]*Op->ExciteSignal[exc_pos];
 			volt[1][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[1][n]*Op->ExciteSignal[exc_pos];
 			volt[2][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[2][n]*Op->ExciteSignal[exc_pos];

FDTD/engine.h

@@ -6,6 +6,7 @@
 class Engine
 {
 	friend class ProcessVoltage;
+	friend class ProcessFieldsTD;
 public:
 	Engine(Operator* op);
 	virtual ~Engine();

FDTD/operator.cpp

@@ -54,27 +54,33 @@ unsigned int Operator::GetNyquistNum(double fmax)
 	return floor(T0/2/dT);
 }
 
-bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord)
+bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower)
 {
 	bool ok=true;
 	for (int n=0;n<3;++n)
 	{
-		if (dcoord[n]<discLines[n][0]) {ok=false;uicoord[n]=0;};
+		uicoord[n]=0;
-		if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1;};
+		if (dcoord[n]<discLines[n][0]) {ok=false;uicoord[n]=0;}
-		for (unsigned int i=0;i<numLines[n]-1;++i)
+		else if (dcoord[n]==discLines[n][0]) {uicoord[n]=0;}
-		{
+		else if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
-			if (dcoord[n]<=discLines[n][i])
+		else if (dcoord[n]==discLines[n][numLines[n]-1]) {uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
+		else
+			for (unsigned int i=0;i<numLines[n]-1;++i)
 			{
-				if (fabs(dcoord[n]-discLines[n][i])<(fabs(dcoord[n]-discLines[n][i+1])))
+				if (dcoord[n]<=discLines[n][i])
-					uicoord[n]=i;
+				{
-				else
+					if (fabs(dcoord[n]-discLines[n][i])<(fabs(dcoord[n]-discLines[n][i+1])))
-					uicoord[n]=i+1;
+						uicoord[n]=i;
-				i = numLines[n];
+					else
+						uicoord[n]=i+1;
+					if (lower) uicoord[n]=i;
+					i = numLines[n];
+				}
 			}
-		}
 	}
-//	cerr << "Operator::SnapToMesh: " << discLines[0][uicoord[0]] << " " << discLines[1][uicoord[1]] << " " << discLines[2][uicoord[2]] << endl;
+//	cerr << "Operator::SnapToMesh Wish: " << dcoord[0] << " " << dcoord[1] << " " << dcoord[2] << endl;
-//	cerr << "Operator::SnapToMesh: " << uicoord[0] << " " << uicoord[1] << " " << uicoord[2] << endl;
+//	cerr << "Operator::SnapToMesh Found: " << discLines[0][uicoord[0]] << " " << discLines[1][uicoord[1]] << " " << discLines[2][uicoord[2]] << endl;
+//	cerr << "Operator::SnapToMesh Index: " << uicoord[0] << " " << uicoord[1] << " " << uicoord[2] << endl;
 	return ok;
 }
 
@@ -99,9 +105,9 @@ void Operator::ShowSize()
 	unsigned int FieldSize = 6*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
 	double MBdiff = 1024*1024;
 
-	cout << "FDTD Operator Size:" << endl;
+	cout << "---- Stat: FDTD Operator ----" << endl;
-	cout << "Size of Operator in Byte  : " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MB) " << endl;
+	cout << "Size of Operator  : " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MB) " << endl;
-	cout << "Size of Field-Data in Byte: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MB) " << endl;
+	cout << "Size of Field-Data: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MB) " << endl;
 }
 
 void Operator::CalcGaussianPulsExcitation(double f0, double fc)

FDTD/operator.h

@@ -11,6 +11,8 @@
 class Operator
 {
 	friend class Engine;
+	friend class ProcessFields;
+	friend class ProcessFieldsTD;
 public:
     Operator();
 	virtual ~Operator();
@@ -34,7 +36,7 @@ public:
 
 	virtual void Reset();
 
-	bool SnapToMesh(double* coord, unsigned int* uicoord);
+	bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false);
 
 protected:
 	virtual void Init();

FDTD/processfields.cpp

@@ -0,0 +1,89 @@
+#include "processfields.h"
+
+ProcessFields::ProcessFields(Operator* op, Engine* eng) : Processing(op, eng)
+{
+	DumpMode=0;
+	DumpType = 0;
+
+	for (int n=0;n<3;++n)
+	{
+		discDLines[n]=NULL;
+	}
+}
+
+ProcessFields::~ProcessFields()
+{
+	for (int n=0;n<3;++n)
+	{
+		delete[] discDLines[n];
+		discDLines[n]=NULL;
+	}
+}
+
+void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
+{
+	if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
+	if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
+
+	for (int n=0;n<3;++n)
+	{
+//		cerr << "start " << start[n] << "stop " << stop[n];
+		if (start[n]>stop[n])
+		{
+			unsigned int help = start[n];
+			start[n]=stop[n];
+			stop[n]=help;
+		}
+		++stop[n];
+		numDLines[n]=stop[n]-start[n];
+//		cerr << " number of lines " << numDLines[n] << endl;
+		delete[] discDLines[n];
+		discDLines[n] = new double[numDLines[n]];
+		for (unsigned int i=0;i<numDLines[n];++i)
+		{
+			discDLines[n][i] = 0.5*(Op->discLines[n][start[n]+i+1] +  Op->discLines[n][start[n]+i]);
+//			cerr << n << " : " << discDLines[n][i] << endl;
+		}
+	}
+}
+
+bool ProcessFields::DumpFieldArray2VTK(ofstream &file, string name, FDTD_FLOAT**** array, double** discLines, unsigned int* numLines)
+{
+	file << "# vtk DataFile Version 2.0" << endl;
+	file << "Rectilinear Grid openEMS_ProcessFields" << endl;
+	file << "ASCII" << endl;
+	file << "DATASET RECTILINEAR_GRID " << endl;
+	file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl;
+	file << "X_COORDINATES " << numLines[0] << " float" << endl;
+	for (unsigned int i=0;i<numLines[0];++i)
+		file << discLines[0][i] << " ";
+	file << endl;
+	file << "Y_COORDINATES " << numLines[1] << " float" << endl;
+	for (unsigned int i=0;i<numLines[1];++i)
+		file << discLines[1][i] << " ";
+	file << endl;
+	file << "Z_COORDINATES " << numLines[2] << " float" << endl;
+	for (unsigned int i=0;i<numLines[2];++i)
+		file << discLines[2][i] << " ";
+	file << endl << endl;
+
+	file << "POINT_DATA " << numLines[0]*numLines[1]*numLines[2] << endl;
+	file << "VECTORS " << name << " float " << endl;
+
+	unsigned int pos[3];
+	for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+	{
+		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+		{
+			for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+			{
+				//in x
+				file << array[0][pos[0]][pos[1]][pos[2]] << " ";
+				//in y
+				file << array[1][pos[0]][pos[1]][pos[2]] << " ";
+				//in z
+				file << array[2][pos[0]][pos[1]][pos[2]] << endl;
+			}
+		}
+	}
+}

FDTD/processfields.h

@@ -0,0 +1,37 @@
+#ifndef PROCESSFIELDS_H
+#define PROCESSFIELDS_H
+
+#include "processing.h"
+#include "../tools/array_ops.h"
+
+class ProcessFields : public Processing
+{
+public:
+	virtual ~ProcessFields();
+
+	virtual void DefineStartStopCoord(double* dstart, double* dstop);
+
+	//! Used file pattern e.g. pattern="tmp/efield_" --> "tmp/efield_000045.vtk" for timestep 45 or "tmp/efield_2.40000e9.vtk" for 2.4GHz E-field dump.
+	void SetFilePattern(string fp) {filePattern=fp;}
+
+	//! This methode will dump all fields in the center of a main cell (dual-node) using 4 E-field and 2 H-fields per direction. (default)
+	void SetDumpMode2Cell() {DumpMode=0;}
+
+	//! Set dump type: 0 for E-fields, 1 for H-fields, 2 for D-fields, 3 for B-fields, 4 for J-fields, etc...
+	void SetDumpType(int type) {DumpType=type;}
+
+//	virtual void Process();
+protected:
+	ProcessFields(Operator* op, Engine* eng);
+
+	bool DumpFieldArray2VTK(ofstream &file, string name, FDTD_FLOAT**** array, double** discLines, unsigned int* numLines);
+
+	int DumpMode;
+	int DumpType;
+	string filePattern;
+
+	unsigned int numDLines[3];
+	double* discDLines[3];
+};
+
+#endif // PROCESSFIELDS_H

FDTD/processfields_td.cpp

@@ -0,0 +1,61 @@
+#include "processfields_td.h"
+#include <iomanip>
+#include <sstream>
+#include <string>
+
+ProcessFieldsTD::ProcessFieldsTD(Operator* op, Engine* eng) : ProcessFields(op, eng)
+{
+	pad_length = 8;
+}
+
+ProcessFieldsTD::~ProcessFieldsTD()
+{
+}
+
+void ProcessFieldsTD::Process()
+{
+	stringstream ss;
+	ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS;
+
+	string filename = filePattern + ss.str() + ".vtk";
+	ofstream file(filename.c_str());
+	if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl; return;};
+
+	if (DumpType==0)
+	{
+		//create array
+		FDTD_FLOAT**** E_T = Create_N_3DArray(numDLines);
+		unsigned int pos[3] = {start[0],start[1],start[2]};
+		unsigned int OpPos[3];
+		double delta;
+//		cerr << "processing e-fields... " << endl;
+		for (pos[0]=0;pos[0]<numDLines[0];++pos[0])
+		{
+			OpPos[0]=start[0]+pos[0];
+			for (pos[1]=0;pos[1]<numDLines[1];++pos[1])
+			{
+				OpPos[1]=start[1]+pos[1];
+				for (pos[2]=0;pos[2]<numDLines[2];++pos[2])
+				{
+					OpPos[2]=start[2]+pos[2];
+					//in x
+					delta  = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
+					E_T[0][pos[0]][pos[1]][pos[2]] = Eng->volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + Eng->volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + Eng->volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1];
+					E_T[0][pos[0]][pos[1]][pos[2]] /= delta;
+					//in y
+					delta  = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
+					E_T[1][pos[0]][pos[1]][pos[2]] = Eng->volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + Eng->volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + Eng->volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1];
+					E_T[1][pos[0]][pos[1]][pos[2]] /= delta;
+					//in z
+					delta  = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
+					E_T[2][pos[0]][pos[1]][pos[2]] = Eng->volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + Eng->volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + Eng->volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]];
+					E_T[2][pos[0]][pos[1]][pos[2]] /= delta;
+				}
+			}
+		}
+		DumpFieldArray2VTK(file,string("E-Field"),E_T,discDLines,numDLines);
+		Delete_N_3DArray(E_T,numDLines);
+		E_T = NULL;
+	}
+	file.close();
+}

FDTD/processfields_td.h

@@ -0,0 +1,21 @@
+#ifndef PROCESSFIELDS_TD_H
+#define PROCESSFIELDS_TD_H
+
+#include "processfields.h"
+
+class ProcessFieldsTD : public ProcessFields
+{
+public:
+	ProcessFieldsTD(Operator* op, Engine* eng);
+	virtual ~ProcessFieldsTD();
+
+	virtual void Process();
+
+	//! Set the length of the filename timestep pad filled with zeros (default is 8)
+	void SetPadLength(int val) {pad_length=val;};
+
+protected:
+	int pad_length;
+};
+
+#endif // PROCESSFIELDS_TD_H

FDTD/processing.cpp

@@ -8,17 +8,12 @@ Processing::Processing(Operator* op, Engine* eng)
 
 Processing::~Processing()
 {
-	file.close();
 }
 
-void Processing::OpenFile(string outfile)
+
+void Processing::DefineStartStopCoord(double* dstart, double* dstop)
 {
-	if (file.is_open()) file.close();
+	if (Op->SnapToMesh(dstart,start)==false) cerr << "Processing::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
-
+	if (Op->SnapToMesh(dstop,stop)==false) cerr << "Processing::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
-	file.open(outfile.c_str());
-	if (file.is_open()==false)
-	{
-		cerr << "Can't open file: " << outfile << endl;
-		return;
-	}
 }
+

FDTD/processing.h

@@ -9,17 +9,19 @@
 class Processing
 {
 public:
-	Processing(Operator* op, Engine* eng);
 	virtual ~Processing();
 
-	virtual void OpenFile(string outfile);
+	virtual void DefineStartStopCoord(double* dstart, double* dstop);
 
 	virtual void Process() {};
 
 protected:
+	Processing(Operator* op, Engine* eng);
 	Operator* Op;
 	Engine* Eng;
-	ofstream file;
+
+	unsigned int start[3];
+	unsigned int stop[3];
 };
 
 #endif // PROCESSING_H

FDTD/processvoltage.cpp

@@ -6,12 +6,19 @@ ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng)
 
 ProcessVoltage::~ProcessVoltage()
 {
+	file.close();
 }
 
-void ProcessVoltage::DefineStartStopCoord(double* dstart, double* dstop)
+void ProcessVoltage::OpenFile(string outfile)
 {
-	if (Op->SnapToMesh(dstart,start)==false) cerr << "ProcessVoltage::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
+	if (file.is_open()) file.close();
-	if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessVoltage::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl;
+
+	file.open(outfile.c_str());
+	if (file.is_open()==false)
+	{
+		cerr << "Can't open file: " << outfile << endl;
+		return;
+	}
 }
 
 void ProcessVoltage::Process()

FDTD/processvoltage.h

@@ -3,19 +3,19 @@
 
 #include "processing.h"
 
+//! Process voltage along a line from start to stop coordinates. ATM integration along the axis e.g.: in x, then y then z direction (Future: diagonal integration)
 class ProcessVoltage : public Processing
 {
 public:
 	ProcessVoltage(Operator* op, Engine* eng);
 	virtual ~ProcessVoltage();
 
-	void DefineStartStopCoord(double* dstart, double* dstop);
+	virtual void OpenFile(string outfile);
 
 	virtual void Process();
 
 protected:
-	unsigned int start[3];
+	ofstream file;
-	unsigned int stop[3];
 
 	vector<FDTD_FLOAT> voltages;
 };

main.cpp

@@ -6,6 +6,7 @@
 #include "FDTD/cartoperator.h"
 #include "FDTD/engine.h"
 #include "FDTD/processvoltage.h"
+#include "FDTD/processfields_td.h"
 
 #include "ContinuousStructure.h"
 
@@ -35,8 +36,6 @@ int main(int argc, char *argv[])
 
 	cop.ShowSize();
 
-//	cop.DumpOperator2File("tmp/Operator");
-
 	cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl;
 	unsigned int NrIter = cop.GetNyquistNum(fmax)/3;
 
@@ -50,18 +49,26 @@ int main(int argc, char *argv[])
 	//*************** setup processing ************//
 	ProcessVoltage PV(&cop,&eng);
 	PV.OpenFile("tmp/u1");
-	double start[]={-500,-150,0};
+	double start[]={-0,-75,0};
-	double stop[]={-500,150,0};
+	double stop[]={-0,75,0};
 	PV.DefineStartStopCoord(start,stop);
 	unsigned int maxIter = 5000;
 
+	ProcessFieldsTD PETD(&cop,&eng);
+	start[0]=-1000;start[1]=0;start[2]=-1000;
+	stop[0]=1000;stop[1]=0;stop[2]=1000;
+	PETD.SetFilePattern("tmp/Et_");
+	PETD.DefineStartStopCoord(start,stop);
+
 	PV.Process();
-//	NrIter=200;
+	PETD.Process();
+
 	//*************** simulate ************//
 	for (unsigned int i=0;i<maxIter;i+=NrIter)
 	{
 		eng.IterateTS(NrIter);
 		PV.Process();
+		PETD.Process();
 	}
 
 	//*************** postproc ************//
@@ -89,7 +96,9 @@ void BuildDipol(ContinuousStructure &CSX)
 
 	CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet());
 	elec->SetExcitation(1,1);
-	elec->SetExcitType(0);
+	elec->SetExcitType(1);
+	elec->SetActiveDir(0,0);//disable x
+	elec->SetActiveDir(0,2);//disable z
 //	elec->SetDelay(2.0e-9);
 	CSX.AddProperty(elec);
 

openEMS.pro

@@ -22,7 +22,9 @@ SOURCES += main.cpp \
     FDTD/operator.cpp \
     tools/array_ops.cpp \
     FDTD/processvoltage.cpp \
-    FDTD/processing.cpp
+    FDTD/processing.cpp \
+    FDTD/processfields.cpp \
+    FDTD/processfields_td.cpp
 HEADERS += FDTD/cartoperator.h \
     tools/ErrorMsg.h \
     tools/AdrOp.h \
@@ -31,4 +33,6 @@ HEADERS += FDTD/cartoperator.h \
     FDTD/operator.h \
     tools/array_ops.h \
     FDTD/processvoltage.h \
-    FDTD/processing.h
+    FDTD/processing.h \
+    FDTD/processfields.h \
+    FDTD/processfields_td.h

tools/array_ops.cpp

@@ -3,7 +3,7 @@
 
 FDTD_FLOAT*** Create3DArray(unsigned int* numLines)
 {
-	FDTD_FLOAT*** array;
+	FDTD_FLOAT*** array=NULL;
 	unsigned int pos[3];
 	array = new FDTD_FLOAT**[numLines[0]];
 	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
@@ -38,7 +38,7 @@ void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines)
 
 FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines)
 {
-	FDTD_FLOAT**** array;
+	FDTD_FLOAT**** array=NULL;
 	array = new FDTD_FLOAT***[3];
 	for (int n=0;n<3;++n)
 	{