diff --git a/Common/processfields.cpp b/Common/processfields.cpp index 60c011c..5e9aa0b 100644 --- a/Common/processfields.cpp +++ b/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 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 const* const* const* const* array, unsigned int const* numLines, float weight, float frequency); diff --git a/Common/processfields_fd.h b/Common/processfields_fd.h index 61d192c..d078137 100644 --- a/Common/processfields_fd.h +++ b/Common/processfields_fd.h @@ -34,7 +34,7 @@ public: virtual void PostProcess(); protected: - void DumpFDData(); + virtual void DumpFDData(); //! frequency domain field storage vector****> m_FD_Fields; diff --git a/Common/processfields_sar.cpp b/Common/processfields_sar.cpp new file mode 100644 index 0000000..76a8f0a --- /dev/null +++ b/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 . +*/ + +#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(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(numLines); + std::complex**** field_fd = NULL; + double kdRatio = 0; + + for (size_t n = 0; nGetInterpolationTypeString(), 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); +} diff --git a/Common/processfields_sar.h b/Common/processfields_sar.h new file mode 100644 index 0000000..adc4975 --- /dev/null +++ b/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 . +*/ + +#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 diff --git a/openEMS.pro b/openEMS.pro index 2123466..05669a4 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -92,9 +92,10 @@ SOURCES += Common/operator_base.cpp \ Common/processintegral.cpp \ Common/processfields.cpp \ Common/processfields_td.cpp \ - Common/processcurrent.cpp \ - Common/processfields_fd.cpp \ - Common/processfieldprobe.cpp + Common/processcurrent.cpp \ + Common/processfields_fd.cpp \ + Common/processfieldprobe.cpp \ + Common/processfields_sar.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ @@ -144,9 +145,10 @@ HEADERS += Common/operator_base.h \ Common/processfields.h \ Common/processfields_td.h \ Common/processcurrent.h \ - Common/processmodematch.h \ - Common/processfields_fd.h \ - Common/processfieldprobe.h + Common/processmodematch.h \ + Common/processfields_fd.h \ + Common/processfieldprobe.h \ + Common/processfields_sar.h QMAKE_CXXFLAGS_RELEASE = -O3 \ -g \ diff --git a/openems.cpp b/openems.cpp index dcd5f6f..9223af5 100644 --- a/openems.cpp +++ b/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 #include #include // only for H5get_libversion() @@ -376,8 +377,10 @@ bool openEMS::SetupProcessing() { if ((db->GetDumpType()>=0) && (db->GetDumpType()<=3)) ProcField = new ProcessFieldsTD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng)); - else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=13)) + 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 ) - FDTD_Op->SetMaterialStoreFlags(1,true); //keep kappa material data (prevent cleanup) + if (db->GetDumpType()==20) + { + 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;