operator: alternative material averaging method

If constant cell material is activated, material probing is performed
only in the center of a primary cell.

This should improve and simplify SAR calculation if all materials
are assumed as constant within a primary YEE cell.

Usage from Matlab/Octave:
FDTD = InitFDTD('CellConstantMaterial',1);

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

FDTD/operator.cpp

@@ -95,6 +95,7 @@ void Operator::Init()
 
 	m_Exc = 0;
 	m_TimeStepFactor = 1;
+	m_MatCellShiftFaktor = 0.25;
 }
 
 void Operator::Delete()
@@ -1220,8 +1221,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	//******************************* epsilon,kappa averaging *****************************//
 	//shift up-right
 	shiftCoord[n] = coord[n]+delta*0.5;
-	shiftCoord[nP] = coord[nP]+deltaP*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*m_MatCellShiftFaktor;
-	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*m_MatCellShiftFaktor;
 	A_n = GetNodeArea(ny,loc_pos,true);
 	EffMat[0] = GetMaterial(n, shiftCoord, 0)*A_n;
 	EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n;
@@ -1229,8 +1230,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 
 	//shift up-left
 	shiftCoord[n] = coord[n]+delta*0.5;
-	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
+	shiftCoord[nP] = coord[nP]-deltaP_M*m_MatCellShiftFaktor;
-	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*m_MatCellShiftFaktor;
 
 	--loc_pos[nP];
 	A_n = GetNodeArea(ny,loc_pos,true);
@@ -1240,8 +1241,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 
 	//shift down-right
 	shiftCoord[n] = coord[n]+delta*0.5;
-	shiftCoord[nP] = coord[nP]+deltaP*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*m_MatCellShiftFaktor;
-	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
+	shiftCoord[nPP] = coord[nPP]-deltaPP_M*m_MatCellShiftFaktor;
 	++loc_pos[nP];
 	--loc_pos[nPP];
 	A_n = GetNodeArea(ny,loc_pos,true);
@@ -1251,8 +1252,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 
 	//shift down-left
 	shiftCoord[n] = coord[n]+delta*0.5;
-	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
+	shiftCoord[nP] = coord[nP]-deltaP_M*m_MatCellShiftFaktor;
-	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
+	shiftCoord[nPP] = coord[nPP]-deltaPP_M*m_MatCellShiftFaktor;
 	--loc_pos[nP];
 	A_n = GetNodeArea(ny,loc_pos,true);
 	EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n;
@@ -1269,7 +1270,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	double length=0;
 
 	//shift down
-	shiftCoord[n] = coord[n]-delta_M*0.25;
+	shiftCoord[n] = coord[n]-delta_M*m_MatCellShiftFaktor;
 	shiftCoord[nP] = coord[nP]+deltaP*0.5;
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
 	--loc_pos[n];
@@ -1283,7 +1284,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	length=delta_ny;
 
 	//shift up
-	shiftCoord[n] = coord[n]+delta*0.25;
+	shiftCoord[n] = coord[n]+delta*m_MatCellShiftFaktor;
 	shiftCoord[nP] = coord[nP]+deltaP*0.5;
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
 	++loc_pos[n];

FDTD/operator.h

@@ -81,6 +81,9 @@ public:
 	virtual void SetTimestepFactor(double factor);
 	bool GetTimestepValid() const {return !m_InvaildTimestep;}
 
+	//! Set operator to assume a constant material inside a cell (material probing in the cell center)
+	void SetCellConstantMaterial() {m_MatCellShiftFaktor=0.5;}
+
 	virtual double GetNumberCells() const;
 
 	virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();}
@@ -204,6 +207,8 @@ protected:
 	//! Get the material at a given coordinate, direction and type from CSX (internal use only)
 	virtual double GetMaterial(int ny, const double coords[3], int MatType, bool markAsUsed=true) const;
 
+	double m_MatCellShiftFaktor;
+
 	//! Calculate the effective/averaged material properties at the given position and direction ny.
 	virtual bool Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const;
 

matlab/InitFDTD.m

@@ -12,6 +12,8 @@ function FDTD = InitFDTD(varargin)
 %   CoordSystem:    choose coordinate system (0 Cartesian, 1 Cylindrical)
 %   TimeStep:       force to use a given timestep (dangerous!)
 %   TimeStepFactor: reduce the timestep by a given factor (>0 to <=1)
+%   CellConstantMaterial: set to 1 to assume a material is constant inside
+%                         a cell (material probing in cell center)
 %
 % examples:
 %   %default init with 1e9 max. timesteps and -50dB end-criteria

openems.cpp

@@ -647,6 +647,15 @@ int openEMS::SetupFDTD(const char* file)
 	if (SetupOperator(FDTD_Opts)==false)
 		return 2;
 
+	int cellConstantMaterial=0;
+	FDTD_Opts->QueryIntAttribute("CellConstantMaterial",&cellConstantMaterial);
+	if (cellConstantMaterial==1)
+	{
+		FDTD_Op->SetCellConstantMaterial();
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Enabling constant cell material assumption." << endl;
+	}
+
 	m_Exc = new Excitation();
 	FDTD_Op->SetExcitationSignal(m_Exc);
 	FDTD_Op->AddExtension(new Operator_Ext_Excitation(FDTD_Op));