new frequency domain dump type: local SAR

Common/processfields.cpp

@@ -67,6 +67,8 @@ string ProcessFields::GetFieldNameByType(DumpType type)
 		return "J-Field";
 	case ROTH_FIELD_DUMP:
 		return "RotH-Field";
+	case SAR_LOCAL_DUMP:
+		return "SAR-local";
 	}
 	return "unknown field";
 }
@@ -470,6 +472,55 @@ bool ProcessFields::DumpVectorArray2HDF5(string filename, string groupName, stri
 	return true;
 }
 
+bool ProcessFields:: DumpScalarArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* array, unsigned int const* numLines, string attr_name, float attr_value)
+{
+	const H5std_string FILE_NAME(filename);
+	const H5std_string DATASET_NAME( name );
+
+	H5::H5File file( FILE_NAME, H5F_ACC_RDWR );
+
+	H5::Group group( file.openGroup( groupName ));
+
+	hsize_t dimsf[3];              // dataset dimensions
+
+	dimsf[0] = numLines[2];
+	dimsf[1] = numLines[1];
+	dimsf[2] = numLines[0];
+
+	H5::DataSpace dataspace( 3, dimsf );
+
+	H5::FloatType datatype( H5::PredType::NATIVE_FLOAT );
+//	datatype.setOrder( H5T_ORDER_LE );
+	H5::DataSet dataset = group.createDataSet( DATASET_NAME, datatype, dataspace );
+
+	if (!attr_name.empty())
+	{
+		hsize_t t_dimsf[] = {1};
+		H5::DataSpace t_dataspace( 1, t_dimsf );
+		H5::Attribute attr = dataset.createAttribute(attr_name,H5::PredType::NATIVE_FLOAT,t_dataspace);
+		attr.write( H5::PredType::NATIVE_FLOAT , &attr_value);
+	}
+
+	// I have not the slightest idea why this array-copy action is necessary...  but it's the only way hdf5 does what it is supposed to do anyway!!
+	// at least it is save in case FDTD_FLOAT was defined as double...
+	// why does hdf5 write the dimensions backwards??? or matlab???
+	unsigned long pos = 0;
+	float *hdf5array = new float[numLines[0]*numLines[1]*numLines[2]];
+	for (unsigned int k=0; k<numLines[2]; ++k)
+	{
+		for (unsigned int j=0; j<numLines[1]; ++j)
+		{
+			for (unsigned int i=0; i<numLines[0]; ++i)
+			{
+				hdf5array[pos++] = array[i][j][k];
+			}
+		}
+	}
+	dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
+	delete[] hdf5array;
+	return true;
+}
+
 bool ProcessFields::DumpVectorArray2HDF5(string filename, string groupName, string name, std::complex<float> const* const* const* const* array, unsigned int const* numLines, float weight, float frequency)
 {
 	const H5std_string FILE_NAME(filename);
@@ -551,6 +602,7 @@ 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];

Common/processfields.h

@@ -37,7 +37,7 @@ public:
 	  Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP),
 	  (conduction) electric current density (kappa*E) (J_FIELD_DUMP) and total current density (rotH)
 	  */
-	enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP, J_FIELD_DUMP, ROTH_FIELD_DUMP};
+	enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP, J_FIELD_DUMP, ROTH_FIELD_DUMP, SAR_LOCAL_DUMP};
 
 	virtual string GetProcessingName() const {return "common field processing";}
 
@@ -89,6 +89,9 @@ public:
 	//! Dump a time-domain vector dump to an HDF5 file
 	static bool DumpVectorArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, float time=0);
 
+	//! Dump a scalar field to an HDF5 file
+	static bool DumpScalarArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* array, unsigned int const* numLines, string attr_name=string(), float attr_value=0);
+
 	//! Dump a frequency-domain complex-vector dump to an HDF5 file
 	static bool DumpVectorArray2HDF5(string filename, string groupName, string name, std::complex<float> const* const* const* const* array, unsigned int const* numLines, float weight, float frequency);
 

Common/processfields_fd.h

@@ -34,7 +34,7 @@ public:
 	virtual void PostProcess();
 
 protected:
-	void DumpFDData();
+	virtual void DumpFDData();
 
 	//! frequency domain field storage
 	vector<std::complex<float>****> m_FD_Fields;

Common/processfields_sar.cpp

@@ -0,0 +1,133 @@
+/*
+*	Copyright (C) 2011 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
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "processfields_sar.h"
+#include "operator_base.h"
+
+ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessFieldsFD(eng_if)
+{
+}
+
+ProcessFieldsSAR::~ProcessFieldsSAR()
+{
+}
+
+void ProcessFieldsSAR::InitProcess()
+{
+	if (m_DumpType!=SAR_LOCAL_DUMP)
+	{
+		Enabled=false;
+		cerr << "ProcessFieldsSAR::InitProcess(): Error, wrong dump type... this should not happen... skipping!" << endl;
+		return;
+	}
+	if (m_Eng_Interface->GetInterpolationType()!=Engine_Interface_Base::NODE_INTERPOLATE)
+	{
+		cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to NODE!" << endl;
+		SetDumpMode2Node();
+	}
+	ProcessFieldsFD::InitProcess();
+}
+
+double ProcessFieldsSAR::GetKappaDensityRatio(const unsigned int* pos)
+{
+
+	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()
+{
+#ifdef OUTPUT_IN_DRAWINGUNITS
+	double discLines_scaling = 1;
+#else
+	double discLines_scaling = Op->GetGridDelta();
+#endif
+
+	unsigned int pos[3];
+	FDTD_FLOAT*** SAR = Create3DArray<float>(numLines);
+	std::complex<float>**** field_fd = NULL;
+	double kdRatio = 0;
+
+	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
+	{
+		field_fd = m_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);
+					SAR[pos[0]][pos[1]][pos[2]] = kdRatio*pow(abs(field_fd[0][pos[0]][pos[1]][pos[2]]) , 2);
+					SAR[pos[0]][pos[1]][pos[2]] = kdRatio*pow(abs(field_fd[1][pos[0]][pos[1]][pos[2]]) , 2);
+					SAR[pos[0]][pos[1]][pos[2]] = kdRatio*pow(abs(field_fd[2][pos[0]][pos[1]][pos[2]]) , 2);
+				}
+			}
+		}
+
+		if (m_fileType==VTK_FILETYPE)
+		{
+			stringstream ss;
+			ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << ".vtk";
+			ofstream file(ss.str().c_str());
+			if (file.is_open()==false)
+				cerr << "ProcessFieldsSAR::DumpFDData: can't open file '" << ss.str() << "' for writing... abort! " << endl;
+			DumpScalarArray2VTK(file,GetFieldNameByType(m_DumpType),SAR,discLines,numLines,m_precision,string("Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString(), m_Mesh_Type, discLines_scaling);
+			file.close();
+		}
+		else if (m_fileType==HDF5_FILETYPE)
+		{
+			stringstream ss;
+			ss << "f" << n;
+			DumpScalarArray2HDF5(m_filename.c_str(), "/FieldData/FD", ss.str(), SAR,numLines , "frequency", m_FD_Samples.at(n));
+		}
+		else
+			cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl;
+	}
+
+	Delete3DArray(SAR,numLines);
+}

Common/processfields_sar.h

@@ -0,0 +1,39 @@
+/*
+*	Copyright (C) 2011 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
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PROCESSFIELDS_SAR_H
+#define PROCESSFIELDS_SAR_H
+
+#include "processfields_fd.h"
+
+class ProcessFieldsSAR : public ProcessFieldsFD
+{
+public:
+	ProcessFieldsSAR(Engine_Interface_Base* eng_if);
+	virtual ~ProcessFieldsSAR();
+
+	virtual string GetProcessingName() const {return "SAR dump";}
+
+	virtual void InitProcess();
+
+protected:
+	virtual void DumpFDData();
+
+	double GetKappaDensityRatio(const unsigned int* pos);
+};
+
+#endif // PROCESSFIELDS_SAR_H

openEMS.pro

@@ -94,7 +94,8 @@ SOURCES += Common/operator_base.cpp \
     Common/processfields_td.cpp \
     Common/processcurrent.cpp \
     Common/processfields_fd.cpp \
-	Common/processfieldprobe.cpp
+	Common/processfieldprobe.cpp \
+	Common/processfields_sar.cpp
 
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
@@ -146,7 +147,8 @@ HEADERS += Common/operator_base.h \
     Common/processcurrent.h \
     Common/processmodematch.h \
     Common/processfields_fd.h \
-	Common/processfieldprobe.h
+	Common/processfieldprobe.h \
+	Common/processfields_sar.h
 
 QMAKE_CXXFLAGS_RELEASE = -O3 \
     -g \

openems.cpp

@@ -33,6 +33,7 @@
 #include "Common/processmodematch.h"
 #include "Common/processfields_td.h"
 #include "Common/processfields_fd.h"
+#include "Common/processfields_sar.h"
 #include <sys/time.h>
 #include <time.h>
 #include <H5Cpp.h> // only for H5get_libversion()
@@ -378,6 +379,8 @@ bool openEMS::SetupProcessing()
 					ProcField = new ProcessFieldsTD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
 				else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=13))
 					ProcField = new ProcessFieldsFD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
+				else if (db->GetDumpType()==20)
+					ProcField = new ProcessFieldsSAR(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
 				else
 					cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl;
 				if (ProcField)
@@ -398,8 +401,14 @@ bool openEMS::SetupProcessing()
 					else
 						ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
 
-					if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12)) && Enable_Dumps )
+					if (db->GetDumpType()==20)
-						FDTD_Op->SetMaterialStoreFlags(1,true); //keep kappa material data (prevent cleanup)
+					{
+						ProcField->SetDumpType(ProcessFields::SAR_LOCAL_DUMP);
+					}
+
+					//SetupMaterialStorages() has previewed storage needs... refresh here to prevent cleanup!!!
+					if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12) || (db->GetDumpType()==20)) && Enable_Dumps )
+						FDTD_Op->SetMaterialStoreFlags(1,true);
 
 					ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode());
 					ProcField->SetFileType((ProcessFields::FileType)db->GetFileType());
@@ -437,7 +446,7 @@ bool openEMS::SetupMaterialStorages()
 		if (db->GetQtyPrimitives()==0)
 			continue;
 		//check for current density dump types
-		if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12)) && Enable_Dumps )
+		if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12) || (db->GetDumpType()==20)) && Enable_Dumps )
 			FDTD_Op->SetMaterialStoreFlags(1,true); //tell operator to store kappa material data
 	}
 	return true;