rewritten SAR calculation

Common/operator_base.h

@@ -69,6 +69,9 @@ public:
 	*/
 	virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const =0;
 
+	//! Get the volume of an FDTD cell
+	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const =0;
+
 	//! Snap the given coodinates to mesh indices, return box dimension
 	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool* inside=NULL) const =0;
 

Common/processfields.cpp

@@ -241,7 +241,6 @@ FDTD_FLOAT**** ProcessFields::CalcField()
 	switch (m_DumpType)
 	{
 	case E_FIELD_DUMP:
-	case SAR_LOCAL_DUMP:
 		for (unsigned int i=0; i<numLines[0]; ++i)
 		{
 			pos[0]=posLines[0][i];
@@ -317,8 +316,9 @@ FDTD_FLOAT**** ProcessFields::CalcField()
 			}
 		}
 		return field;
+	default:
+		cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;
+		return field;
 	}
-	cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;
-	return field;
 }
 

Common/processfields_sar.cpp

@@ -26,6 +26,25 @@ ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessField
 
 ProcessFieldsSAR::~ProcessFieldsSAR()
 {
+	for (size_t n = 0; n<m_FD_Fields.size(); ++n)
+	{
+		Delete_N_3DArray(m_E_FD_Fields.at(n),numLines);
+		Delete_N_3DArray(m_J_FD_Fields.at(n),numLines);
+	}
+	m_E_FD_Fields.clear();
+	m_J_FD_Fields.clear();
+}
+
+void ProcessFieldsSAR::SetSubSampling(unsigned int subSampleRate, int dir)
+{
+	cerr << "ProcessFieldsSAR::SetSubSampling: Warning: Defining a sub-sampling for SAR calculation is not recommended!!! Expect false results!" << endl;
+	ProcessFieldsFD::SetSubSampling(subSampleRate,dir);
+}
+
+void ProcessFieldsSAR::SetOptResolution(double optRes, int dir)
+{
+	cerr << "ProcessFieldsSAR::SetOptResolution: Warning: Defining a sub-sampling for SAR calculation is not recommended!!! Expect false results!" << endl;
+	ProcessFieldsFD::SetOptResolution(optRes,dir);
 }
 
 void ProcessFieldsSAR::InitProcess()
@@ -36,71 +55,143 @@ void ProcessFieldsSAR::InitProcess()
 		cerr << "ProcessFieldsSAR::InitProcess(): Error, wrong dump type... this should not happen... skipping!" << endl;
 		return;
 	}
-	if (m_Eng_Interface->GetInterpolationType()!=Engine_Interface_Base::NODE_INTERPOLATE)
+	if (m_Eng_Interface->GetInterpolationType()!=Engine_Interface_Base::CELL_INTERPOLATE)
 	{
-		cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to NODE!" << endl;
+		cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to CELL!" << endl;
-		SetDumpMode2Node();
+		SetDumpMode2Cell();
 	}
 	ProcessFieldsFD::InitProcess();
-}
 
-double ProcessFieldsSAR::GetKappaDensityRatio(const unsigned int* pos)
+	if (Enabled==false) return;
-{
-
-	double coord[3] = {discLines[0][pos[0]],discLines[1][pos[1]],discLines[2][pos[2]]};
-	unsigned int OpPos[] = {posLines[0][pos[0]],posLines[1][pos[1]],posLines[2][pos[2]]};
-	ContinuousStructure* CSX = Op->GetGeometryCSX();
-	CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
-
-	if (prop==0)
-		return 0.0;
-	CSPropMaterial* matProp = dynamic_cast<CSPropMaterial*>(prop);
-	double density = matProp->GetDensityWeighted(coord);
-
-	if (density==0)
-		return 0.0;
-
-	double kappa = 0;
-	double max_kappa = 0;
-
-	for (int n=0;n<3;++n)
-	{
-		kappa = Op->GetDiscMaterial(1,n,OpPos);
-		if (kappa>max_kappa)
-			max_kappa = kappa;
-		if (OpPos[n]>0)
-		{
-			--OpPos[n];
-			kappa = Op->GetDiscMaterial(1,n,OpPos);
-			if (kappa>max_kappa)
-				max_kappa = kappa;
-			++OpPos[n];
-		}
-	}
-
-	return max_kappa/density;
-}
-
-void ProcessFieldsSAR::DumpFDData()
-{
-	unsigned int pos[3];
-	FDTD_FLOAT*** SAR = Create3DArray<float>(numLines);
-	std::complex<float>**** field_fd = NULL;
-	double kdRatio = 0;
 
+	//create data structures...
 	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
 	{
-		field_fd = m_FD_Fields.at(n);
+		m_E_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
+		m_J_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
+	}
+}
+
+int ProcessFieldsSAR::Process()
+{
+	if (Enabled==false) return -1;
+	if (CheckTimestep()==false) return GetNextInterval();
+
+	if ((m_FD_Interval==0) || (m_Eng_Interface->GetNumberOfTimesteps()%m_FD_Interval!=0))
+		return GetNextInterval();
+
+	std::complex<float>**** field_fd = NULL;
+	unsigned int pos[3];
+	double T;
+	FDTD_FLOAT**** field_td=NULL;
+
+	// calc E-field
+	m_DumpType = E_FIELD_DUMP;
+	field_td = CalcField();
+	T = m_Eng_Interface->GetTime(m_dualTime);
+	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
+	{
+		std::complex<float> exp_jwt_2_dt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
+		exp_jwt_2_dt *= 2; // *2 for single-sided spectrum
+		exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval
+		field_fd = m_E_FD_Fields.at(n);
 		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])
 				{
-					kdRatio = GetKappaDensityRatio(pos);
+					field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
-					SAR[pos[0]][pos[1]][pos[2]] = kdRatio*pow(abs(field_fd[0][pos[0]][pos[1]][pos[2]]) , 2);
+					field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
-					SAR[pos[0]][pos[1]][pos[2]] = kdRatio*pow(abs(field_fd[1][pos[0]][pos[1]][pos[2]]) , 2);
+					field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
-					SAR[pos[0]][pos[1]][pos[2]] = kdRatio*pow(abs(field_fd[2][pos[0]][pos[1]][pos[2]]) , 2);
+				}
+			}
+		}
+	}
+	Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
+
+	// calc J-field
+	m_DumpType = J_FIELD_DUMP;
+	field_td = CalcField();
+	T = m_Eng_Interface->GetTime(m_dualTime);
+	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
+	{
+		std::complex<float> exp_jwt_2_dt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
+		exp_jwt_2_dt *= 2; // *2 for single-sided spectrum
+		exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval
+		field_fd = m_J_FD_Fields.at(n);
+		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])
+				{
+					field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+					field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+					field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+				}
+			}
+		}
+	}
+	Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
+
+	//reset dump type
+	m_DumpType = SAR_LOCAL_DUMP;
+
+	++m_FD_SampleCount;
+	return GetNextInterval();
+}
+
+void ProcessFieldsSAR::DumpFDData()
+{
+	unsigned int pos[3];
+	unsigned int orig_pos[3];
+	FDTD_FLOAT*** SAR = Create3DArray<float>(numLines);
+	std::complex<float>**** E_field_fd = NULL;
+	std::complex<float>**** J_field_fd = NULL;
+	double coord[3];
+	double density;
+	ContinuousStructure* CSX = Op->GetGeometryCSX();
+	CSProperties* prop = NULL;
+	CSPropMaterial* matProp = NULL;
+
+	double power;
+
+	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
+	{
+		E_field_fd = m_E_FD_Fields.at(n);
+		J_field_fd = m_J_FD_Fields.at(n);
+		power = 0;
+		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);
+					prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
+					SAR[pos[0]][pos[1]][pos[2]] = 0.0;
+					density=0.0;
+					if (prop!=0)
+					{
+						matProp = dynamic_cast<CSPropMaterial*>(prop);
+						density = matProp->GetDensityWeighted(coord);
+						if (density>0)
+						{
+							SAR[pos[0]][pos[1]][pos[2]]  = abs(E_field_fd[0][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[0][pos[0]][pos[1]][pos[2]]);
+							SAR[pos[0]][pos[1]][pos[2]] += abs(E_field_fd[1][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[1][pos[0]][pos[1]][pos[2]]);
+							SAR[pos[0]][pos[1]][pos[2]] += abs(E_field_fd[2][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[2][pos[0]][pos[1]][pos[2]]);
+							SAR[pos[0]][pos[1]][pos[2]] *= 0.5/density;
+						}
+						else
+							density=0.0;
+					}
+					power+=SAR[pos[0]][pos[1]][pos[2]]*Op->GetCellVolume(orig_pos)*density;
 				}
 			}
 		}
@@ -126,6 +217,9 @@ void ProcessFieldsSAR::DumpFDData()
 			float freq[1]={m_FD_Samples.at(n)};
 			if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false)
 				cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
+			float pow[1]={power};
+			if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false)
+				cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
 		}
 		else
 			cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl;

Common/processfields_sar.h

@@ -30,10 +30,19 @@ public:
 
 	virtual void InitProcess();
 
+	virtual int Process();
+
+	virtual void SetSubSampling(unsigned int subSampleRate, int dir=-1);
+
+	virtual void SetOptResolution(double optRes, int dir=-1);
+
 protected:
 	virtual void DumpFDData();
 
-	double GetKappaDensityRatio(const unsigned int* pos);
+	//! frequency domain electric field storage
+	vector<std::complex<float>****> m_E_FD_Fields;
+	//! frequency domain current density storage
+	vector<std::complex<float>****> m_J_FD_Fields;
 };
 
 #endif // PROCESSFIELDS_SAR_H

FDTD/operator.cpp

@@ -181,6 +181,14 @@ double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) co
 	}
 }
 
+double Operator::GetCellVolume(const unsigned int pos[3], bool dualMesh) const
+{
+	double vol=1;
+	for (int n=0;n<3;++n)
+		vol*=GetEdgeLength(n,pos,dualMesh);
+	return vol;
+}
+
 double Operator::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
 {
 	if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
@@ -1338,6 +1346,7 @@ bool Operator::Calc_LumpedElements()
 								EC_C[ny][ipos] = epsilon * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
 							if (R>0)
 								EC_G[ny][ipos] = kappa * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
+
 							if (R==0) //make lumped element a PEC if resistance is zero
 							{
 								SetVV(ny,pos[0],pos[1],pos[2], 0 );

FDTD/operator.h

@@ -111,6 +111,9 @@ public:
 	//! Get the length of an FDTD edge (unit is meter).
 	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const;
 
+	//! Get the volume of an FDTD cell
+	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const;
+
 	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
 	/*!
 		This will return the area around an edge with a given direction, measured at the middle of the edge.

FDTD/operator_cylinder.cpp

@@ -176,6 +176,11 @@ double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool
 	return length * GetDiscLine(0,pos[0],dualMesh);
 }
 
+double Operator_Cylinder::GetCellVolume(const unsigned int pos[3], bool dualMesh) const
+{
+	return GetEdgeArea(2,pos,dualMesh)*GetEdgeLength(2,pos,dualMesh);
+}
+
 double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh) const
 {
 	if (ny!=0)

FDTD/operator_cylinder.h

@@ -58,6 +58,9 @@ public:
 	//! Get the length of an FDTD edge, including radius corrected alpha-mesh width.
 	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const;
 
+	//! Get the volume of an FDTD cell
+	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const;
+
 	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
 	/*!
 		This will return the area around an edge with a given direction, measured at the middle of the edge.