A couple of new features

- ProcessCurrent
- Weighted Excitation and Material
- Testing

Analyse/PlotVoltage.m

@@ -2,22 +2,62 @@
 clear all;
 clc
 
-tmp = load('../tmp/u1');
+figure(1);
+tmpu = load('../tmp/u1');
+tmpi = load('../tmp/i1');
 
-t = tmp(:,1);
+t = tmpu(:,1);
-u = tmp(:,2);
+u = tmpu(:,2);
 
-L=numel(t);
+subplot(2,2,1);
-
+title('u_1 TD');
-subplot(2,1,1);
 plot(t,u);
+xlabel('t \rightarrow');
+ylabel('ut_1 \rightarrow');
 grid on;
 
 dt=t(2)-t(1);
+u= [u ; zeros(size(u))];
+L=numel(u);
+t = (1:L)*dt;
 
 f = (1:L)/L/dt;
 fu = fft(u)/L;
-subplot(2,1,2);
+subplot(2,2,2);
+title('u_1 FD');
 plot(f(1:L/2),abs(fu(1:L/2)));
+xlabel('f \rightarrow');
+ylabel('|uf_1| \rightarrow');
 grid on;
 
+
+t = tmpi(:,1);
+i = tmpi(:,2);
+
+subplot(2,2,3);
+title('i_1 TD');
+plot(t,i);
+xlabel('t \rightarrow');
+ylabel('it_1 \rightarrow');
+grid on;
+
+dt=t(2)-t(1);
+i = [i; zeros(size(t))];
+L=numel(i);
+t = (1:L)*dt;
+f = (1:L)/L/dt;
+
+fi = fft(i)/L;
+subplot(2,2,4);
+title('i_1 FD');
+plot(f(1:L/2),abs(fi(1:L/2)));
+xlabel('f \rightarrow');
+ylabel('|if_1| \rightarrow');
+grid on;
+
+figure(2);
+plot(f,real(fu./fi));
+xlim([0 1e9]);
+grid on;
+
+

FDTD/cartoperator.cpp

@@ -204,8 +204,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] = mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
-		inEC[1] = mat->GetKappa(n)*fabs(deltaP*deltaPP);
+		inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
 	}
 	else
 	{
@@ -220,8 +220,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
-		inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP);
+		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
 	}
 	else
 	{
@@ -237,8 +237,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
-		inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP);
+		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
 	}
 	else
 	{
@@ -254,8 +254,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
-		inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP);
+		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
 	}
 	else
 	{
@@ -275,9 +275,9 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[2] = fabs(delta_M) / mat->GetMue(n);
+		inEC[2] = fabs(delta_M) / mat->GetMueWeighted(n,shiftCoord);
 		if (mat->GetSigma(n))
-			inEC[3] = fabs(delta_M) / mat->GetSigma(n);
+			inEC[3] = fabs(delta_M) / mat->GetSigmaWeighted(n,shiftCoord);
 		else
 			inEC[3] = 0;
 	}
@@ -295,8 +295,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
 		inEC[2] += mat->GetMue(n)*fabs(delta);
-		if (mat->GetSigma(n))
+		if (mat->GetSigmaWeighted(n,shiftCoord))
-			inEC[3] += fabs(delta)/mat->GetSigma(n);
+			inEC[3] += fabs(delta)/mat->GetSigmaWeighted(n,shiftCoord);
 		else
 			inEC[3] = 0;
 	}
@@ -454,8 +454,10 @@ bool CartOperator::CalcEFieldExcitation()
 								vExcit[n].push_back(0);
 							if ((elec->GetExcitType()==1) && (elec->GetActiveDir(n))) //hard excite
 							{
-								vv[n][pos[0]][pos[1]][pos[2]] = 0;
+								vv[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0;
-								vi[n][pos[0]][pos[1]][pos[2]] = 0;
+								vi[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0;
+								vv[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0;
+								vi[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0;
 							}
 						}
 					}

FDTD/engine.h

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

FDTD/processcurrent.cpp

@@ -0,0 +1,71 @@
+#include "processcurrent.h"
+
+ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng)
+{
+}
+
+ProcessCurrent::~ProcessCurrent()
+{
+}
+
+void ProcessCurrent::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;
+	}
+}
+
+void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
+{
+	if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl;
+	if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
+}
+
+void ProcessCurrent::Process()
+{
+	if (Enabled==false) return;
+	FDTD_FLOAT current=0;
+//	FDTD_FLOAT help=0;
+	double sign[3]={1,1,1};
+	unsigned int pos[3]={start[0],start[1],start[2]};
+	double loc_start[]={start[0],start[1],start[2]};
+	double loc_stop[]={stop[0],stop[1],stop[2]};
+
+	for (int n=0;n<3;++n)
+	{
+		if (start[n]>stop[n])
+		{
+			unsigned int help=start[n];
+			start[n]=stop[n];
+			stop[n]=help;
+		}
+	}
+
+	//x-current
+	for (int i=start[0];i<stop[0];++i)
+		current+=Eng->curr[0][i][start[1]][start[2]];
+	//y-current
+	for (int i=start[1];i<stop[1];++i)
+		current+=Eng->curr[1][stop[0]][i][start[2]];
+	//z-current
+	for (int i=start[2];i<stop[2];++i)
+		current+=Eng->curr[2][stop[0]][stop[1]][i];
+	//x-current
+	for (int i=start[0];i<stop[0];++i)
+		current-=Eng->curr[0][i][stop[1]][stop[2]];
+	//y-current
+	for (int i=start[1];i<stop[1];++i)
+		current-=Eng->curr[1][start[0]][i][stop[2]];
+	//z-current
+	for (int i=start[2];i<stop[2];++i)
+		current-=Eng->curr[2][start[0]][start[1]][i];
+
+//	cerr << "ts: " << Eng->numTS << " i: " << current << endl;
+	v_current.push_back(current);
+	file << (double)Eng->numTS*Op->GetTimestep() << "\t" << current << endl;
+}

FDTD/processcurrent.h

@@ -0,0 +1,24 @@
+#ifndef PROCESSCURRENT_H
+#define PROCESSCURRENT_H
+
+#include "processing.h"
+
+class ProcessCurrent : public Processing
+{
+public:
+	ProcessCurrent(Operator* op, Engine* eng);
+	virtual ~ProcessCurrent();
+
+	virtual void OpenFile(string outfile);
+
+	virtual void DefineStartStopCoord(double* dstart, double* dstop);
+
+	virtual void Process();
+
+protected:
+	ofstream file;
+
+	vector<FDTD_FLOAT> v_current;
+};
+
+#endif // PROCESSCURRENT_H

main.cpp

@@ -6,6 +6,7 @@
 #include "FDTD/cartoperator.h"
 #include "FDTD/engine.h"
 #include "FDTD/processvoltage.h"
+#include "FDTD/processcurrent.h"
 #include "FDTD/processfields_td.h"
 
 #include "ContinuousStructure.h"
@@ -21,23 +22,26 @@ int main(int argc, char *argv[])
 	time_t startTime=time(NULL);
 
 	//*************** setup/read geometry ************//
+	cerr << "Create Geometry..." << endl;
 	ContinuousStructure CSX;
 
-	BuildPlaneWave(CSX);
+	//BuildPlaneWave(CSX);
+	BuildMSL(CSX);
 
 	//*************** setup operator ************//
+	cerr << "Create Operator..." << endl;
 	CartOperator cop;
 	cop.SetGeometryCSX(&CSX);
 
 	cop.CalcECOperator();
 
 	double fmax=1e9;
-	cop.CalcGaussianPulsExcitation(fmax/2,fmax/8);
+	cop.CalcGaussianPulsExcitation(fmax/2,fmax/2);
 
 	time_t OpDoneTime=time(NULL);
 
 	cop.ShowSize();
-	bool bnd[] = {1,1,0,0,0,0};
+	bool bnd[] = {1,1,0,1,0,1};
 	cop.ApplyMagneticBC(bnd);
 
 	cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl;
@@ -53,16 +57,23 @@ int main(int argc, char *argv[])
 	//*************** setup processing ************//
 	ProcessVoltage PV(&cop,&eng);
 	PV.OpenFile("tmp/u1");
-	double start[]={-100,-250,0};
+	double start[]={0,50,0};
-	double stop[]={-100,250,0};
+	double stop[]={0,0,0};
 	PV.DefineStartStopCoord(start,stop);
-	unsigned int maxIter = 2000;
+
+	ProcessCurrent PCurr(&cop,&eng);
+	PCurr.OpenFile("tmp/i1");
+	start[0]=-50;start[1]=40;start[2]=-0;
+	stop[0]=50;stop[1]=60;stop[2]=-0;
+	PCurr.DefineStartStopCoord(start,stop);
+
+	unsigned int maxIter = 1800;
 
 	ProcessFieldsTD PETD(&cop,&eng);
-	start[0]=-480;start[1]=0;start[2]=-4000;
+	start[0]=-500;start[1]=25;start[2]=-800;
-	stop[0]=480;stop[1]=0;stop[2]=4000;
+	stop[0] = 500;stop[1] =24;stop[2] = 800;
-//	start[0]=-250;start[1]=-250;start[2]=0;
+//	start[0]=-300;start[1]=0;start[2]=-4000;
-//	stop[0]=250;stop[1]=250;stop[2]=0;
+//	stop[0]=300;stop[1]=300;stop[2]=-4000;
 	PETD.SetDumpType(0);
 	PETD.SetFilePattern("tmp/Et_");
 	PETD.DefineStartStopCoord(start,stop);
@@ -78,6 +89,7 @@ int main(int argc, char *argv[])
 	PHTD.SetEnable(false);
 
 	PV.Process();
+	PCurr.Process();
 	PETD.Process();
 	PHTD.Process();	
 
@@ -86,6 +98,7 @@ int main(int argc, char *argv[])
 	{
 		eng.IterateTS(NrIter);
 		PV.Process();
+		PCurr.Process();
 		PETD.Process();
 		PHTD.Process();
 	}
@@ -97,7 +110,7 @@ int main(int argc, char *argv[])
 	double t_diff = difftime(currTime,prevTime);
 
 	cerr << "Time for " << eng.GetNumberOfTimesteps() << " iterations with " << cop.GetNumberCells() << " cells : " << t_diff << " sec" << endl;
-	cerr << "Speed (MCells/s): " << (double)cop.GetNumberCells()*(double)eng.GetNumberOfTimesteps()/t_diff/1e6 << endl;
+	cerr << "Speed: " << (double)cop.GetNumberCells()*(double)eng.GetNumberOfTimesteps()/t_diff/1e6 << " MCells/s " << endl;
     return 0;
 }
 
@@ -192,46 +205,73 @@ void BuildPlaneWave(ContinuousStructure &CSX)
 
 void BuildMSL(ContinuousStructure &CSX)
 {
-//	CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
+	double width = 1000;
-//	mat->SetEpsilon(3.6);
+	double hight = 500;
-//	CSX.AddProperty(mat);
+	double length = 2000;
-//
+	double abs_l = 200;
-//	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat);
-//	box->SetCoord(0,-200.0);box->SetCoord(1,200.0);
-//	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
-//	box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0);
-//	CSX.AddPrimitive(box);
-//
-//	CSPropMaterial* MSL = new CSPropMaterial(CSX.GetParameterSet());
-//	MSL->SetKappa(56e6);
-//	CSX.AddProperty(MSL);
-//
-//	box = new CSPrimBox(CSX.GetParameterSet(),MSL);
-//	box->SetCoord(0,-20.0);box->SetCoord(1,20.0);
-//	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
-//	box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0);
-//	CSX.AddPrimitive(box);
 
+	//substrate....
+	CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
+//	mat->SetEpsilon(3.6);
+	double finalKappa = 3/pow(abs_l,4);
+	mat->SetKappa(finalKappa);
+	std::ostringstream fct;
+	fct << "pow(abs(z)-" << length/2.0-abs_l << ",4)";
+	mat->SetKappaWeightFunction(fct.str(),0);
+	mat->SetKappaWeightFunction(fct.str(),1);
+	mat->SetKappaWeightFunction(fct.str(),2);
+	mat->SetSigma(finalKappa*__MUE0__/__EPS0__);
+	mat->SetSigmaWeightFunction(fct.str(),0);
+	mat->SetSigmaWeightFunction(fct.str(),1);
+	mat->SetSigmaWeightFunction(fct.str(),2);
+	CSX.AddProperty(mat);
+	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat);
+	box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0);
+	box->SetCoord(2,0.0);box->SetCoord(3,hight);
+	box->SetCoord(4,length/2.0-abs_l); box->SetCoord(5,length/2.0);
+	CSX.AddPrimitive(box);
+	box = new CSPrimBox(CSX.GetParameterSet(),mat);
+	box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0);
+	box->SetCoord(2,0.0);box->SetCoord(3,hight);
+	box->SetCoord(4,length/-2.0+abs_l); box->SetCoord(5,length/-2.0);
+	CSX.AddPrimitive(box);
+
+//	double coords[]={0,0,length/2};
+//	cerr << fct.str() << endl;
+//	cerr << finalKappa*pow(abs_l,4) << endl;
+//	cerr << mat->GetKappaWeighted(0,coords) << endl;
+//	exit(0);
+
+	//MSL
+	CSPropMetal* MSL = new CSPropMetal(CSX.GetParameterSet());
+//	MSL->SetKappa(56e6);
+	CSX.AddProperty(MSL);
+	box = new CSPrimBox(CSX.GetParameterSet(),MSL);
+	box->SetCoord(0,-40.0);box->SetCoord(1,40.0);
+	box->SetCoord(2,50.0);box->SetCoord(3,50.0);
+	box->SetCoord(4,length/-2);box->SetCoord(5,length/2.0);
+	CSX.AddPrimitive(box);
+
+	//MSL excite...
 	CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet());
-	elec->SetExcitation(1,1);
+	elec->SetExcitation(-1,1);
-	elec->SetExcitType(0);
+	elec->SetExcitType(1);
 //	elec->SetDelay(2.0e-9);
 	CSX.AddProperty(elec);
-
+	box = new CSPrimBox(CSX.GetParameterSet(),elec);
-	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec);
+	box->SetCoord(0,-40.0);box->SetCoord(1,40.0);
-	box->SetCoord(0,-1.0);box->SetCoord(1,1.0);
 	box->SetCoord(2,0.0);box->SetCoord(3,50.0);
-	box->SetCoord(4,-1.0);box->SetCoord(5,1.0);
+	box->SetCoord(4,length/-2.0);box->SetCoord(5,length/-2.0);
 	CSX.AddPrimitive(box);
 
 	CSRectGrid* grid = CSX.GetGrid();
 
-	for (int n=-300;n<=300;n+=20)
+	for (double n=width/-2.0;n<=width/2;n+=20)
-		grid->AddDiscLine(0,(double)n);
+		grid->AddDiscLine(0,n);
-	for (int n=0;n<=300;n+=20)
+	for (double n=0;n<=500;n+=10)
-		grid->AddDiscLine(1,(double)n);
+		grid->AddDiscLine(1,n);
-	for (int n=-1000;n<=1000;n+=20)
+	for (double n=length/-2.0;n<=length/2.0;n+=20)
-		grid->AddDiscLine(2,(double)n);
+		grid->AddDiscLine(2,n);
 
 	grid->SetDeltaUnit(1e-3);
 

openEMS.pro

@@ -7,7 +7,8 @@ CONFIG += console
 CONFIG -= app_bundle
 TEMPLATE = app
 OBJECTS_DIR = obj
-INCLUDEPATH += ../CSXCAD
+INCLUDEPATH += ../CSXCAD \
+			../fparser
 LIBS += -L../CSXCAD \
     -lCSXCAD \
     -L../fparser \
@@ -24,7 +25,8 @@ SOURCES += main.cpp \
     FDTD/processvoltage.cpp \
     FDTD/processing.cpp \
     FDTD/processfields.cpp \
-    FDTD/processfields_td.cpp
+    FDTD/processfields_td.cpp \
+    FDTD/processcurrent.cpp
 HEADERS += FDTD/cartoperator.h \
     tools/ErrorMsg.h \
     tools/AdrOp.h \
@@ -35,4 +37,5 @@ HEADERS += FDTD/cartoperator.h \
     FDTD/processvoltage.h \
     FDTD/processing.h \
     FDTD/processfields.h \
-    FDTD/processfields_td.h
+    FDTD/processfields_td.h \
+    FDTD/processcurrent.h

tools/constants.h

@@ -3,6 +3,8 @@
 
 #define __EPS0__ 8.85418781762e-12
 #define __MUE0__ 1.256637062e-6
+#define __C0__ 299792458
+#define __Z0__ 376.730313461
 #define PI 3.141592653589793238462643383279
 
 #endif // CONSTANTS_H