From 84ba0104824da8295afc7eafa46d53b0982322fe Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <thorsten.liebig@gmx.de>
Date: Tue, 2 Mar 2010 22:55:50 +0100
Subject: [PATCH] PEC support added

---
 FDTD/cartoperator.cpp | 37 +++++++++++++++++++++++++++++++++++++
 FDTD/cartoperator.h   |  1 +
 main.cpp              | 30 ++++++++++++++++++++----------
 3 files changed, 58 insertions(+), 10 deletions(-)

diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp
index bc402e7..f3ab6a6 100644
--- a/FDTD/cartoperator.cpp
+++ b/FDTD/cartoperator.cpp
@@ -115,6 +115,7 @@ int CartOperator::CalcECOperator()
 	ApplyElectricBC(PEC);
 
 	if (CalcEFieldExcitation()==false) return -1;
+	CalcPEC();
 
     return 0;
 }
@@ -483,3 +484,39 @@ bool CartOperator::CalcEFieldExcitation()
 	}
 	return true;
 }
+
+bool CartOperator::CalcPEC()
+{
+	unsigned int pos[3];
+	double coord[3];
+	double delta;
+
+	for (int n=0;n<3;++n)
+	{
+		for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+		{
+			coord[2] = discLines[2][pos[2]];
+			for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+			{
+				coord[1] = discLines[1][pos[1]];
+				for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+				{
+					coord[0] = discLines[0][pos[0]];
+					MainOp->SetPos(pos[0],pos[1],pos[2]);
+					delta=MainOp->GetIndexDelta(n,pos[n]);
+					coord[n]= discLines[n][pos[n]] + delta*0.5;
+					CSProperties* prop = CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
+					if (prop)
+					{
+						if (prop->GetType()==CSProperties::METAL) //set to PEC
+						{
+							vv[n][pos[0]][pos[1]][pos[2]] = 0;
+							vi[n][pos[0]][pos[1]][pos[2]] = 0;
+//							cerr << "CartOperator::CalcPEC: PEC found at " << pos[0] << " ; "  << pos[1] << " ; " << pos[2] << endl;
+						}
+					}
+				}
+			}
+		}
+	}
+}
diff --git a/FDTD/cartoperator.h b/FDTD/cartoperator.h
index 911a014..3a0f63a 100644
--- a/FDTD/cartoperator.h
+++ b/FDTD/cartoperator.h
@@ -25,6 +25,7 @@ protected:
     AdrOp* DualOp;
 
 	virtual bool CalcEFieldExcitation();
+	virtual bool CalcPEC();
 	virtual double CalcTimestep();
 
 	//EC elements, internal only!
diff --git a/main.cpp b/main.cpp
index 53fcf5a..d372e71 100644
--- a/main.cpp
+++ b/main.cpp
@@ -22,6 +22,7 @@ int main(int argc, char *argv[])
 
 	//*************** setup/read geometry ************//
 	ContinuousStructure CSX;
+
 	BuildPlaneWave(CSX);
 
 	//*************** setup operator ************//
@@ -58,8 +59,8 @@ int main(int argc, char *argv[])
 	unsigned int maxIter = 2000;
 
 	ProcessFieldsTD PETD(&cop,&eng);
-	start[0]=-250;start[1]=0;start[2]=-2000;
-	stop[0]=250;stop[1]=0;stop[2]=4000;
+	start[0]=-480;start[1]=0;start[2]=-4000;
+	stop[0]=480;stop[1]=0;stop[2]=4000;
 //	start[0]=-250;start[1]=-250;start[2]=0;
 //	stop[0]=250;stop[1]=250;stop[2]=0;
 	PETD.SetDumpType(0);
@@ -73,7 +74,7 @@ int main(int argc, char *argv[])
 	PHTD.SetFilePattern("tmp/Ht_");
 	PHTD.DefineStartStopCoord(start,stop);
 
-	PETD.SetEnable(false);
+	PETD.SetEnable(true);
 	PHTD.SetEnable(false);
 
 	PV.Process();
@@ -149,7 +150,7 @@ void BuildPlaneWave(ContinuousStructure &CSX)
 //	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,-2000.0);matbox->SetCoord(5,4000.0);
+//	matbox->SetCoord(4,-4000.0);matbox->SetCoord(5,4000.0);
 //	CSX.AddPrimitive(matbox);
 
 	CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet());
@@ -161,18 +162,27 @@ void BuildPlaneWave(ContinuousStructure &CSX)
 	CSX.AddProperty(elec);
 
 	CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec);
-	box->SetCoord(0,-250.0);box->SetCoord(1,250.0);
-	box->SetCoord(2,-250.0);box->SetCoord(3,250.0);
-	box->SetCoord(4,-0000.0);box->SetCoord(5,-0000.0);
+	box->SetCoord(0,-500.0);box->SetCoord(1,500.0);
+	box->SetCoord(2,-500.0);box->SetCoord(3,500.0);
+	box->SetCoord(4,-4000.0);box->SetCoord(5,-4000.0);
 	CSX.AddPrimitive(box);
 
+	CSPropMetal* metal = new CSPropMetal(CSX.GetParameterSet());
+	CSX.AddProperty(metal);
+	CSPrimCylinder* cyl = new CSPrimCylinder(CSX.GetParameterSet(),metal);
+	cyl->SetRadius(100);
+	cyl->SetCoord(0,0.0);cyl->SetCoord(1,0.0);
+	cyl->SetCoord(2,-250.0);cyl->SetCoord(3,250.0);
+	cyl->SetCoord(4,-0000.0);cyl->SetCoord(5,-0000.0);
+	CSX.AddPrimitive(cyl);
+
 	CSRectGrid* grid = CSX.GetGrid();
 
-	for (int n=-250;n<=250;n+=10)
+	for (int n=-500;n<=500;n+=20)
 		grid->AddDiscLine(0,(double)n);
-	for (int n=-250;n<=250;n+=10)
+	for (int n=-500;n<=500;n+=20)
 		grid->AddDiscLine(1,(double)n);
-	for (int n=-4000;n<=4000;n+=10)
+	for (int n=-4000;n<=4000;n+=20)
 		grid->AddDiscLine(2,(double)n);
 
 	grid->SetDeltaUnit(1e-3);