Operator: get access to cell center material averaging coordinates

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

Common/operator_base.cpp

@@ -85,6 +85,18 @@ void Operator_Base::SetMaterialStoreFlags(int type, bool val)
 	m_StoreMaterial[type]=val;
 }
 
+bool Operator_Base::GetCellCenterMaterialAvgCoord(const int pos[3], double coord[3]) const
+{
+	unsigned int l_pos[3];
+	for (int n=0;n<3;++n)
+	{
+		if (pos[n]<0)
+			return false;
+		else
+			l_pos[n] = pos[n];
+	}
+	return GetCellCenterMaterialAvgCoord(l_pos, coord);
+}
 
 void Operator_Base::SetBackgroundMaterial(double epsR, double mueR, double kappa, double sigma, double density)
 {

Common/operator_base.h

@@ -108,6 +108,12 @@ public:
 	//! Get stored discrete material (if storage is enabled).
 	virtual double GetDiscMaterial(int type, int ny, const unsigned int pos[3]) const = 0;
 
+	//! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center)
+	virtual bool GetCellCenterMaterialAvgCoord(const int pos[3], double coord[3]) const;
+
+	//! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center)
+	virtual bool GetCellCenterMaterialAvgCoord(const unsigned int pos[3], double coord[3]) const = 0;
+
 	//! Set the background material (default is vacuum)
 	virtual void SetBackgroundMaterial(double epsR=0, double mueR=0, double kappa=0, double sigma=0, double density=0);
 

Common/processfields_sar.cpp

@@ -200,19 +200,17 @@ void ProcessFieldsSAR::DumpFDData()
 	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 	{
 		orig_pos[0] = posLines[0][pos[0]];
-		coord[0] = Op->GetDiscLine(0,orig_pos[0],true);
 		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
 		{
 			orig_pos[1] = posLines[1][pos[1]];
-			coord[1] = Op->GetDiscLine(1,orig_pos[1],true);
 			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
 			{
 				orig_pos[2] = posLines[2][pos[2]];
-				coord[2] = Op->GetDiscLine(2,orig_pos[2],true);
 
 				cell_volume[pos[0]][pos[1]][pos[2]] = Op->GetCellVolume(orig_pos);
 				cell_density[pos[0]][pos[1]][pos[2]] = 0.0;
 
+				Op->GetCellCenterMaterialAvgCoord(orig_pos, coord);
 				prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
 				if (prop!=0)
 				{

FDTD/operator.cpp

@@ -192,7 +192,7 @@ bool Operator::GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool du
 	return true;
 }
 
-bool Operator::GetNodeCoords(unsigned int pos[3], double* coords, bool dualMesh, CoordinateSystem c_system) const
+bool Operator::GetNodeCoords(const unsigned int pos[3], double* coords, bool dualMesh, CoordinateSystem c_system) const
 {
 	for (int n=0;n<3;++n)
 		coords[n]=GetDiscLine(n,pos[n],dualMesh);
@@ -1210,6 +1210,17 @@ double Operator::GetRawDiscDelta(int ny, const int pos) const
 	return (discLines[ny][pos+1] - discLines[ny][pos]);
 }
 
+bool Operator::GetCellCenterMaterialAvgCoord(const unsigned int pos[3], double coord[3]) const
+{
+	for (int n=0;n<3;++n)
+	{
+		if (pos[n]>=numLines[n])
+			return false;
+	}
+	GetNodeCoords(pos, coord, true);
+	return true;
+}
+
 double Operator::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const
 {
 	CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed);

FDTD/operator.h

@@ -107,7 +107,7 @@ public:
 	//! Get the coordinates for a given node index and component, according to the yee-algorithm. Returns true if inside the FDTD domain.
 	virtual bool GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const;
 
-	virtual bool GetNodeCoords(unsigned int pos[3], double* coords, bool dualMesh=false, CoordinateSystem c_system=UNDEFINED_CS) const;
+	virtual bool GetNodeCoords(const unsigned int pos[3], double* coords, bool dualMesh=false, CoordinateSystem c_system=UNDEFINED_CS) const;
 
 	//! Get the node width for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,pos,!dualMesh);}
@@ -160,6 +160,9 @@ public:
 
 	virtual double GetDiscMaterial(int type, int ny, const unsigned int pos[3]) const;
 
+	//! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center)
+	virtual bool GetCellCenterMaterialAvgCoord(const unsigned int pos[3], double coord[3]) const;
+
 	virtual void SetExcitationSignal(Excitation* exc);
 	virtual Excitation* GetExcitationSignal() const {return m_Exc;}
 

FDTD/operator_cylindermultigrid.cpp

@@ -239,6 +239,15 @@ void Operator_CylinderMultiGrid::FillMissingDataStorage()
 	}
 }
 
+bool Operator_CylinderMultiGrid::GetCellCenterMaterialAvgCoord(const unsigned int pos[3], double coord[3]) const
+{
+	if (pos[0]>(m_Split_Pos-1))
+		return Operator_Cylinder::GetCellCenterMaterialAvgCoord(pos, coord);
+
+	unsigned int bg_pos[3] = {pos[0],pos[1]/2,pos[2]};
+	return m_InnerOp->GetCellCenterMaterialAvgCoord(bg_pos, coord);
+}
+
 int Operator_CylinderMultiGrid::CalcECOperator( DebugFlags debugFlags )
 {
 	int retCode=0;

FDTD/operator_cylindermultigrid.h

@@ -59,6 +59,9 @@ public:
 
 	virtual void ShowStat() const;
 
+	//! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center)
+	virtual bool GetCellCenterMaterialAvgCoord(const unsigned int pos[3], double coord[3]) const;
+
 #ifdef MPI_SUPPORT
 	virtual void SetTag(int tag);
 	virtual void SetNeighborUp(int ny, int id);