nf2ff: add boundary mirroring

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

matlab/CalcNF2FF.m

@@ -30,6 +30,11 @@ function nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, theta, phi, varargin)
 % 'Mue_r':   specify the relative magnetic permeability for the nf2ff
 % 'MPI'    : set true if MPI was used
 %
+% 'Mirror':  Add mirroring in a given direction (dir), with a given 
+%            mirror type (PEC or PMC) and a mirror position in the given
+%            direction.
+%            Example: 'Mirror', {0, 'PMC', +100}
+%
 % See also: CreateNF2FFBox, ReadNF2FF
 %
 % openEMS matlab interface
@@ -56,6 +61,15 @@ for n=1:2:numel(varargin)-1
         mode = varargin{n+1};
     elseif (strcmp(varargin{n},'MPI'))
         MPI = varargin{n+1};
+    elseif (strcmp(varargin{n},'Mirror'))
+        if isfield(nf2ff_xml,'Mirror')
+            pos = length(nf2ff_xml.Mirror)+1;
+        else
+            pos = 1;
+        end
+        nf2ff_xml.Mirror{pos}.ATTRIBUTE.Dir=varargin{n+1}{1};
+        nf2ff_xml.Mirror{pos}.ATTRIBUTE.Type=varargin{n+1}{2};
+        nf2ff_xml.Mirror{pos}.ATTRIBUTE.Pos=varargin{n+1}{3};
     else
         nf2ff_xml.ATTRIBUTE.(varargin{n})=varargin{n+1};
     end

nf2ff/main.cpp

@@ -1,5 +1,5 @@
 /*
-*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*	Copyright (C) 2012-2014 Thorsten Liebig (Thorsten.Liebig@gmx.de)
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
@@ -25,7 +25,7 @@ int main(int argc, char *argv[])
 {
 	cout << " ---------------------------------------------------------------------- " << endl;
 	cout << " | nf2ff, near-field to far-field transformation for openEMS " << endl;
-	cout << " | (C) 2012 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license" << endl;
+	cout << " | (C) 2012-2014 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license" << endl;
 	cout << " ---------------------------------------------------------------------- " << endl;
 
 	if (argc<=1)

nf2ff/nf2ff.cpp

@@ -1,5 +1,5 @@
 /*
-*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*	Copyright (C) 2012-2014 Thorsten Liebig (Thorsten.Liebig@gmx.de)
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
@@ -121,7 +121,14 @@ void nf2ff::SetPermeability(vector<float> permeability)
 	}
 	for (size_t fn=0;fn<m_nf2ff.size();++fn)
 		m_nf2ff.at(fn)->SetPermeability(permeability.at(fn));
+}
 
+void nf2ff::SetMirror(int type, int dir, float pos)
+{
+	if (m_Verbose>0)
+		cerr << "Enable mirror of type: "<< type << " in direction: " << dir << " at: " << pos << endl;
+	for (size_t fn=0;fn<m_nf2ff.size();++fn)
+		m_nf2ff.at(fn)->SetMirror(type, dir, pos);
 }
 
 bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
@@ -224,7 +231,26 @@ bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
 	if (ti_nf2ff->QueryFloatAttribute("Radius",&radius) ==  TIXML_SUCCESS)
 		l_nf2ff->SetRadius(radius);
 
-	TiXmlElement* ti_Planes = ti_nf2ff->FirstChildElement();
+	// read mirrors
+	TiXmlElement* ti_Mirros = ti_nf2ff->FirstChildElement("Mirror");
+	int dir=-1;
+	string type;
+	float pos=0.0;
+	while (ti_Mirros!=NULL)
+	{
+		type = string(ti_Mirros->Attribute("Type"));
+		if (ti_Mirros->QueryIntAttribute("Dir",&dir) != TIXML_SUCCESS)
+			dir = -1;
+		if (ti_Mirros->QueryFloatAttribute("Pos",&pos) != TIXML_SUCCESS)
+			dir = -1;
+		if ((dir>=0) && (strcmp(type.c_str(),"PEC")==0))
+			l_nf2ff->SetMirror(MIRROR_PEC, dir, pos);
+		else if ((dir>=0) && (strcmp(type.c_str(),"PMC")==0))
+			l_nf2ff->SetMirror(MIRROR_PMC, dir, pos);
+		ti_Mirros = ti_Mirros->NextSiblingElement("Mirror");
+	}
+	
+	TiXmlElement* ti_Planes = ti_nf2ff->FirstChildElement("Planes");
 	string E_name;
 	string H_name;
 	while (ti_Planes!=NULL)

nf2ff/nf2ff.h

@@ -1,5 +1,5 @@
 /*
-*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*	Copyright (C) 2012-2014 Thorsten Liebig (Thorsten.Liebig@gmx.de)
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
@@ -40,6 +40,8 @@ public:
 	void SetRadius(float radius);
 	void SetPermittivity(vector<float> permittivity);
 	void SetPermeability(vector<float> permeability);
+	
+	void SetMirror(int type, int dir, float pos);
 
 	double GetTotalRadPower(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetTotalRadPower();}
 	double GetMaxDirectivity(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetMaxDirectivity();}

nf2ff/nf2ff_calc.cpp

@@ -1,5 +1,5 @@
 /*
-*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*	Copyright (C) 2012-2014 Thorsten Liebig (Thorsten.Liebig@gmx.de)
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
@@ -200,6 +200,12 @@ nf2ff_calc::nf2ff_calc(float freq, vector<float> theta, vector<float> phi, vecto
 	m_maxDir = 0;
 	m_radius = 1;
 
+	for (int n=0;n<3;++n)
+	{
+		m_MirrorType[n] = MIRROR_OFF;
+		m_MirrorPos[n]  = 0.0;
+	}
+
 	m_Barrier = NULL;
 	m_numThreads = boost::thread::hardware_concurrency();
 }
@@ -227,9 +233,8 @@ nf2ff_calc::~nf2ff_calc()
 	m_Barrier = NULL;
 }
 
-bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType)
+int nf2ff_calc::GetNormalDir(unsigned int* numLines)
 {
-	//find normal direction
 	int ny = -1;
 	int nP,nPP;
 	for (int n=0;n<3;++n)
@@ -239,13 +244,116 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 		if ((numLines[n]==1) && (numLines[nP]>2) && (numLines[nPP]>2))
 			ny=n;
 	}
-	nP = (ny+1)%3;
+	return ny;
-	nPP = (ny+2)%3;
+}
+
+void nf2ff_calc::SetMirror(int type, int dir, float pos)
+{
+	if ((dir<0) || (dir>3))
+	{
+		cerr << "nf2ff_calc::SetMirror: Error, invalid direction!" << endl;
+		return;
+	}
+	if ((type!=MIRROR_PEC) && (type!=MIRROR_PMC))
+	{
+		cerr << "nf2ff_calc::SetMirror: Error, invalid type!" << endl;
+		return;
+	}
+	m_MirrorType[dir] = type;
+	m_MirrorPos[dir] = pos;
+}
+
+bool nf2ff_calc::AddMirrorPlane(int n, float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType)
+{
+	float E_factor[3] = {1,1,1};
+	float H_factor[3] = {1,1,1};
+
+	int nP  = (n+1)%3;
+	int nPP = (n+2)%3;
+	
+	// mirror in ny direction
+	for (unsigned int i=0;i<numLines[n];++i)
+		lines[n][i] = 2.0*m_MirrorPos[n] - lines[n][i];
+	if (m_MirrorType[n]==MIRROR_PEC)
+	{
+		H_factor[n]  =-1.0;
+		E_factor[nP] =-1.0;
+		E_factor[nPP]=-1.0;
+	}
+	else if (m_MirrorType[n]==MIRROR_PMC)
+	{
+		E_factor[n]  = -1.0;
+		H_factor[nP] = -1.0;
+		H_factor[nPP]= -1.0;
+	}
+
+	for (int d=0;d<3;++d)
+		for (unsigned int i=0;i<numLines[0];++i)
+			for (unsigned int j=0;j<numLines[1];++j)
+				for (unsigned int k=0;k<numLines[2];++k)
+				{
+					E_field[d][i][j][k] *= E_factor[d];
+					H_field[d][i][j][k] *= H_factor[d];
+				}
+
+	return this->AddSinglePlane(lines, numLines, E_field, H_field, MeshType);
+}
+
+bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType)
+{
+	this->AddSinglePlane(lines, numLines, E_field, H_field, MeshType);
+
+	for (int n=0;n<3;++n)
+	{
+		int nP  = (n+1)%3;
+		int nPP = (n+2)%3;
+		// check if a single mirror plane is on
+		if ((m_MirrorType[n]!=MIRROR_OFF) && (m_MirrorType[nP]==MIRROR_OFF) && (m_MirrorType[nPP]==MIRROR_OFF))
+		{
+			cerr << "single plane in " << n << endl;
+			this->AddMirrorPlane(n, lines, numLines, E_field, H_field, MeshType);
+			break;
+		}
+		//check if two planes are on 
+		else if ((m_MirrorType[n]==MIRROR_OFF) && (m_MirrorType[nP]!=MIRROR_OFF) && (m_MirrorType[nPP]!=MIRROR_OFF))
+		{
+			cerr << "two planes in " << nP << " and " << nPP << endl;
+			this->AddMirrorPlane(nP, lines, numLines, E_field, H_field, MeshType);
+			this->AddMirrorPlane(nPP, lines, numLines, E_field, H_field, MeshType);
+			this->AddMirrorPlane(nP, lines, numLines, E_field, H_field, MeshType);
+			break;
+		}
+	}
+	// check if all planes are on
+	if ((m_MirrorType[0]!=MIRROR_OFF) && (m_MirrorType[1]!=MIRROR_OFF) && (m_MirrorType[2]!=MIRROR_OFF))
+	{
+		cerr << "all three planes on " << endl;
+		this->AddMirrorPlane(0, lines, numLines, E_field, H_field, MeshType);
+		this->AddMirrorPlane(1, lines, numLines, E_field, H_field, MeshType);
+		this->AddMirrorPlane(0, lines, numLines, E_field, H_field, MeshType);
+		this->AddMirrorPlane(2, lines, numLines, E_field, H_field, MeshType);
+		this->AddMirrorPlane(0, lines, numLines, E_field, H_field, MeshType);
+		this->AddMirrorPlane(1, lines, numLines, E_field, H_field, MeshType);
+		this->AddMirrorPlane(0, lines, numLines, E_field, H_field, MeshType);
+	}
+
+	//cleanup E- & H-Fields
+	Delete_N_3DArray(E_field,numLines);
+	Delete_N_3DArray(H_field,numLines);
+	return true;
+}
+
+bool nf2ff_calc::AddSinglePlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType)
+{
+	//find normal direction
+	int ny = this->GetNormalDir(numLines);
 	if (ny<0)
 	{
 		cerr << "nf2ff_calc::AddPlane: Error can't determine normal direction..." << endl;
 		return false;
 	}
+	int nP  = (ny+1)%3;
+	int nPP = (ny+2)%3;
 
 	complex<float>**** Js = Create_N_3DArray<complex<float> >(numLines);
 	complex<float>**** Ms = Create_N_3DArray<complex<float> >(numLines);
@@ -259,15 +367,15 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 
 	float edge_length_P[numLines[nP]];
 	for (unsigned int n=1;n<numLines[nP]-1;++n)
-		edge_length_P[n]=0.5*(lines[nP][n+1]-lines[nP][n-1]);
+		edge_length_P[n]=0.5*fabs(lines[nP][n+1]-lines[nP][n-1]);
-	edge_length_P[0]=0.5*(lines[nP][1]-lines[nP][0]);
+	edge_length_P[0]=0.5*fabs(lines[nP][1]-lines[nP][0]);
-	edge_length_P[numLines[nP]-1]=0.5*(lines[nP][numLines[nP]-1]-lines[nP][numLines[nP]-2]);
+	edge_length_P[numLines[nP]-1]=0.5*fabs(lines[nP][numLines[nP]-1]-lines[nP][numLines[nP]-2]);
 
 	float edge_length_PP[numLines[nPP]];
 	for (unsigned int n=1;n<numLines[nPP]-1;++n)
-		edge_length_PP[n]=0.5*(lines[nPP][n+1]-lines[nPP][n-1]);
+		edge_length_PP[n]=0.5*fabs(lines[nPP][n+1]-lines[nPP][n-1]);
-	edge_length_PP[0]=0.5*(lines[nPP][1]-lines[nPP][0]);
+	edge_length_PP[0]=0.5*fabs(lines[nPP][1]-lines[nPP][0]);
-	edge_length_PP[numLines[nPP]-1]=0.5*(lines[nPP][numLines[nPP]-1]-lines[nPP][numLines[nPP]-2]);
+	edge_length_PP[numLines[nPP]-1]=0.5*fabs(lines[nPP][numLines[nPP]-1]-lines[nPP][numLines[nPP]-2]);
 
 	//check for cylindrical mesh
 	if (MeshType==1)
@@ -342,10 +450,6 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 
 	m_Barrier->wait(); //combine all thread local Nt,Np,Lt and Lp
 
-	//cleanup E- & H-Fields
-	Delete_N_3DArray(E_field,numLines);
-	Delete_N_3DArray(H_field,numLines);
-
 	complex<float>** Nt = Create2DArray<complex<float> >(numAngles);
 	complex<float>** Np = Create2DArray<complex<float> >(numAngles);
 	complex<float>** Lt = Create2DArray<complex<float> >(numAngles);

nf2ff/nf2ff_calc.h

@@ -1,5 +1,5 @@
 /*
-*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*	Copyright (C) 2012-2014 Thorsten Liebig (Thorsten.Liebig@gmx.de)
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
@@ -29,6 +29,10 @@ using namespace std;
 
 class nf2ff_calc;
 
+#define MIRROR_OFF 0
+#define MIRROR_PEC 1
+#define MIRROR_PMC 2
+
 // data structure to exchange data between thread-controller and worker-threads
 typedef struct
 {
@@ -89,6 +93,8 @@ public:
 	unsigned int GetNumThreads() const {return m_numThreads;}
 	void SetNumThreads(unsigned int n) {m_numThreads=n;}
 
+	void SetMirror(int type, int dir, float pos);
+
 	bool AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType=0);
 
 protected:
@@ -113,6 +119,15 @@ protected:
 	float* m_theta;
 	float* m_phi;
 
+	//mirror settings
+	bool m_EnableMirror;
+	int m_MirrorType[3];
+	float m_MirrorPos[3];
+
+	int GetNormalDir(unsigned int* numLines);
+	bool AddSinglePlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType=0);
+	bool AddMirrorPlane(int n, float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType=0);
+
 	//boost multi-threading
 	unsigned int m_numThreads;
 	boost::thread_group m_thread_group;

tools/array_ops.h

@@ -114,6 +114,19 @@ T*** Create3DArray(const unsigned int* numLines)
 	return array;
 }
 
+template <typename T>
+T*** Copy3DArray(T*** array_in, T*** array_out, const unsigned int* numLines)
+{
+	if (array_out==NULL)
+		array_out = Create3DArray<T>(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])
+				array_out[pos[0]][pos[1]][pos[2]] = array_in[pos[0]][pos[1]][pos[2]];
+	return array_out;
+}
+
 template <typename T>
 T**** Create_N_3DArray(const unsigned int* numLines)
 {
@@ -126,6 +139,16 @@ T**** Create_N_3DArray(const unsigned int* numLines)
 	return array;
 }
 
+template <typename T>
+T**** Copy_N_3DArray(T**** array_in, T**** array_out, const unsigned int* numLines)
+{
+	if (array_out==NULL)
+		array_out = Create_N_3DArray<T>(numLines);
+	for (int n=0; n<3; ++n)
+		array_out[n]=Copy3DArray<T>( array_in[n], array_out[n], numLines);
+	return array_out;
+}
+
 template <typename T>
 void Delete3DArray(T*** array, const unsigned int* numLines)
 {