From ea2823377d8de43865ce296229f3d9686d4a8115 Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <thorsten.liebig@gmx.de>
Date: Sat, 27 Mar 2010 23:05:45 +0100
Subject: [PATCH] added Curve for excitation and dump excitation in
 operator-dump

---
 FDTD/operator.cpp | 93 ++++++++++++++++++++++++++++++++++++++++++-----
 FDTD/operator.h   |  4 +-
 openems.cpp       | 33 ++++++++++++-----
 3 files changed, 109 insertions(+), 21 deletions(-)

diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index 514f571..7f7006e 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -106,7 +106,7 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower)
 		else if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
 		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=1;i<numLines[n]-1;++i)
+			for (unsigned int i=1;i<numLines[n];++i)
 			{
 				if (dcoord[n]<discLines[n][i])
 				{
@@ -224,9 +224,9 @@ void Operator::ShowSize()
 	cout << "-----------------------------" << endl;
 }
 
-void Operator::CalcGaussianPulsExcitation(double f0, double fc)
+bool Operator::CalcGaussianPulsExcitation(double f0, double fc)
 {
-	if (dT==0) return;
+	if (dT==0) return false;
 
 	ExciteLength = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
 	cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
@@ -238,14 +238,16 @@ void Operator::CalcGaussianPulsExcitation(double f0, double fc)
 		ExciteSignal[n] = cos(2.0*PI*f0*(n*dT-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*n*dT/3.0-3,2));
 //		cerr << ExciteSignal[n] << endl;
 	}
+	return true;
 }
 
-void Operator::CalcSinusExcitation(double f0, int nTS)
+bool Operator::CalcSinusExcitation(double f0, int nTS)
 {
-	if (dT==0) return;
-	if (nTS<=0) return;
+	if (dT==0) return false;
+	if (nTS<=0) return false;
 
 	ExciteLength = (unsigned int)(nTS);
+	cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
 	delete[] ExciteSignal;
 	ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
 	ExciteSignal[0]=0.0;
@@ -254,6 +256,7 @@ void Operator::CalcSinusExcitation(double f0, int nTS)
 		ExciteSignal[n] = sin(2.0*PI*f0*n*dT);
 //		cerr << ExciteSignal[n] << endl;
 	}
+	return true;
 }
 
 void Operator::DumpOperator2File(string filename)
@@ -266,10 +269,18 @@ void Operator::DumpOperator2File(string filename)
 		return;
 	}
 
-	string names[] = {"vv", "vi", "iv" , "ii"};
-	FDTD_FLOAT**** array[] = {vv,vi,iv,ii};
+	FDTD_FLOAT**** exc = Create_N_3DArray(numLines);
+	for (unsigned int n=0;n<E_Exc_Count;++n)
+	{
+		exc[E_Exc_dir[n]][E_Exc_index[0][n]][E_Exc_index[1][n]][E_Exc_index[2][n]] = E_Exc_amp[n];
+	}
 
-	ProcessFields::DumpMultiVectorArray2VTK(file, names , array , 4, discLines, numLines);
+	string names[] = {"vv", "vi", "iv" , "ii", "exc"};
+	FDTD_FLOAT**** array[] = {vv,vi,iv,ii,exc};
+
+	ProcessFields::DumpMultiVectorArray2VTK(file, names , array , 5, discLines, numLines);
+
+	Delete_N_3DArray(exc,numLines);
 
 	file.close();
 }
@@ -774,6 +785,70 @@ bool Operator::CalcEFieldExcitation()
 		}
 	}
 
+	//special treatment for primitives of type curve (treated as wires) see also Calc_PEC
+	double p1[3];
+	double p2[3];
+	double deltaN=0.0;
+	int n;
+	struct Grid_Path path;
+	CSPropElectrode* elec=NULL;
+	CSProperties* prop=NULL;
+	vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::ELECTRODE);
+	for (size_t p=0;p<vec_prop.size();++p)
+	{
+		prop = vec_prop.at(p);
+		elec = prop->ToElectrode();
+		for (size_t n=0;n<prop->GetQtyPrimitives();++n)
+		{
+			CSPrimitives* prim = prop->GetPrimitive(n);
+			CSPrimCurve* curv = prim->ToCurve();
+			if (curv)
+			{
+				for (size_t i=1;i<curv->GetNumberOfPoints();++i)
+				{
+					curv->GetPoint(i-1,p1);
+					curv->GetPoint(i,p2);
+					path = FindPath(p1,p2);
+					for (size_t t=0;t<path.dir.size();++t)
+					{
+						n = path.dir.at(t);
+						pos[0] = path.posPath[0].at(t);
+						pos[1] = path.posPath[1].at(t);
+						pos[2] = path.posPath[2].at(t);
+						MainOp->SetPos(pos[0],pos[1],pos[2]);
+						deltaN=fabs(MainOp->GetIndexDelta(n,pos[n]));
+						coord[0] = discLines[0][pos[0]];
+						coord[1] = discLines[1][pos[1]];
+						coord[2] = discLines[2][pos[2]];
+						coord[n] += 0.5*deltaN;
+//						cerr << n << " " << coord[0] << " " << coord[1] << " " << coord[2] << endl;
+						if (elec!=NULL)
+						{
+							if ((elec->GetActiveDir(n)) && (pos[n]<(int)numLines[n]-1))
+							{
+								amp = elec->GetWeightedExcitation(n,coord)*deltaN*gridDelta;
+								if (amp!=0)
+								{
+									vExcit.push_back(amp);
+									vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
+									vDir.push_back(n);
+									vIndex[0].push_back(pos[0]);
+									vIndex[1].push_back(pos[1]);
+									vIndex[2].push_back(pos[2]);
+								}
+								if (elec->GetExcitType()==1) //hard excite
+								{
+									vv[n][pos[0]][pos[1]][pos[2]] = 0;
+									vi[n][pos[0]][pos[1]][pos[2]] = 0;
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
 	E_Exc_Count = vExcit.size();
 	cerr << "Operator::CalcEFieldExcitation: Found number of excitations points: " << E_Exc_Count << endl;
 	if (E_Exc_Count==0)
diff --git a/FDTD/operator.h b/FDTD/operator.h
index fc3047a..cc25e6e 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -38,9 +38,9 @@ public:
 	virtual int CalcECOperator();
 
 	//! Calculate an excitation with center of f0 and the half bandwidth fc
-	virtual void CalcGaussianPulsExcitation(double f0, double fc);
+	virtual bool CalcGaussianPulsExcitation(double f0, double fc);
 	//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps
-	virtual void CalcSinusExcitation(double f0, int nTS);
+	virtual bool CalcSinusExcitation(double f0, int nTS);
 
 	virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
 	virtual void ApplyMagneticBC(bool* dirs);
diff --git a/openems.cpp b/openems.cpp
index fe9790e..a9eb635 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -185,26 +185,39 @@ int openEMS::SetupFDTD(const char* file)
 	FDTD_Op = new Operator();
 	if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1);
 
-	if (DebugMat)
-	{
-		FDTD_Op->DumpMaterial2File("material_dump.vtk");
-	}
 	FDTD_Op->CalcECOperator();
-	if (DebugOp)
-	{
-		FDTD_Op->DumpOperator2File("operator_dump.vtk");
-	}
 
 	if (Excit_Type==0)
-		FDTD_Op->CalcGaussianPulsExcitation(f0,fc);
+	{
+		if (!FDTD_Op->CalcGaussianPulsExcitation(f0,fc))
+		{
+			cerr << "openEMS: excitation setup failed!!" << endl;
+			exit(2);
+		}
+	}
 	else if (Excit_Type==1)
-		FDTD_Op->CalcSinusExcitation(f0,NrTS);
+	{
+		if (!FDTD_Op->CalcSinusExcitation(f0,NrTS))
+		{
+			cerr << "openEMS: excitation setup failed!!" << endl;
+			exit(2);
+		}
+	}
 	else
 	{
 		cerr << "openEMS: Excitation type is unknown" << endl;
 		exit(-1);
 	}
 
+	if (DebugMat)
+	{
+		FDTD_Op->DumpMaterial2File("material_dump.vtk");
+	}
+	if (DebugOp)
+	{
+		FDTD_Op->DumpOperator2File("operator_dump.vtk");
+	}
+
 	time_t OpDoneTime=time(NULL);
 
 	FDTD_Op->ShowSize();