1 /*
2 * Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include "processfields_sar.h"
19 #include "operator_base.h"
20 #include "tools/vtk_file_writer.h"
21 #include "tools/hdf5_file_writer.h"
22 #include "tools/sar_calculation.h"
23
24 #include "CSPropMaterial.h"
25
26 using namespace std;
27
ProcessFieldsSAR(Engine_Interface_Base * eng_if)28 ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessFieldsFD(eng_if)
29 {
30 m_UseCellKappa = true;
31 m_SAR_method = "Simple";
32 }
33
~ProcessFieldsSAR()34 ProcessFieldsSAR::~ProcessFieldsSAR()
35 {
36 for (size_t n = 0; n<m_E_FD_Fields.size(); ++n)
37 Delete_N_3DArray(m_E_FD_Fields.at(n),numLines);
38 m_E_FD_Fields.clear();
39
40 for (size_t n = 0; n<m_J_FD_Fields.size(); ++n)
41 Delete_N_3DArray(m_J_FD_Fields.at(n),numLines);
42 m_J_FD_Fields.clear();
43 }
44
SetDumpType(DumpType type)45 void ProcessFieldsSAR::SetDumpType(DumpType type)
46 {
47 if (type==SAR_RAW_DATA)
48 m_UseCellKappa = true;
49 ProcessFieldsFD::SetDumpType(type);
50 }
51
NeedConductivity() const52 bool ProcessFieldsSAR::NeedConductivity() const
53 {
54 return !m_UseCellKappa;
55 }
56
SetSubSampling(unsigned int subSampleRate,int dir)57 void ProcessFieldsSAR::SetSubSampling(unsigned int subSampleRate, int dir)
58 {
59 UNUSED(subSampleRate);UNUSED(dir);
60 cerr << "ProcessFieldsSAR::SetSubSampling: Warning: Defining a sub-sampling for SAR calculation is not allowed!!! Skipped!" << endl;
61 }
62
SetOptResolution(double optRes,int dir)63 void ProcessFieldsSAR::SetOptResolution(double optRes, int dir)
64 {
65 UNUSED(optRes);UNUSED(dir);
66 cerr << "ProcessFieldsSAR::SetOptResolution: Warning: Defining a sub-sampling for SAR calculation is not allowed!!! Skipped!" << endl;
67 }
68
InitProcess()69 void ProcessFieldsSAR::InitProcess()
70 {
71 if ((m_DumpType!=SAR_LOCAL_DUMP) && (m_DumpType!=SAR_1G_DUMP) && (m_DumpType!=SAR_10G_DUMP) && (m_DumpType!=SAR_RAW_DATA))
72 {
73 Enabled=false;
74 cerr << "ProcessFieldsSAR::InitProcess(): Error, wrong dump type... this should not happen... skipping!" << endl;
75 return;
76 }
77 if (m_Eng_Interface->GetInterpolationType()!=Engine_Interface_Base::CELL_INTERPOLATE)
78 {
79 cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to CELL!" << endl;
80 SetDumpMode2Cell();
81 }
82
83 if ((m_DumpType==SAR_RAW_DATA) && (m_fileType!=HDF5_FILETYPE))
84 {
85 Enabled=false;
86 cerr << "ProcessFieldsSAR::InitProcess(): Error, wrong file type for dumping raw SAR data! skipping" << endl;
87 return;
88
89 }
90
91 ProcessFieldsFD::InitProcess();
92
93 if (Enabled==false) return;
94
95 //create data structures...
96 for (size_t n = 0; n<m_FD_Samples.size(); ++n)
97 {
98 m_E_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
99 if (!m_UseCellKappa)
100 m_J_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
101 }
102 }
103
Process()104 int ProcessFieldsSAR::Process()
105 {
106 if (Enabled==false) return -1;
107 if (CheckTimestep()==false) return GetNextInterval();
108
109 if ((m_FD_Interval==0) || (m_Eng_Interface->GetNumberOfTimesteps()%m_FD_Interval!=0))
110 return GetNextInterval();
111
112 std::complex<float>**** field_fd = NULL;
113 unsigned int pos[3];
114 double T;
115 FDTD_FLOAT**** field_td=NULL;
116
117 //save dump type
118 DumpType save_dump_type = m_DumpType;
119
120 // calc E-field
121 m_DumpType = E_FIELD_DUMP;
122 field_td = CalcField();
123 T = m_Eng_Interface->GetTime(m_dualTime);
124 for (size_t n = 0; n<m_FD_Samples.size(); ++n)
125 {
126 std::complex<float> exp_jwt_2_dt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
127 exp_jwt_2_dt *= 2; // *2 for single-sided spectrum
128 exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval
129 field_fd = m_E_FD_Fields.at(n);
130 for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
131 {
132 for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
133 {
134 for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
135 {
136 field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
137 field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
138 field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
139 }
140 }
141 }
142 }
143 Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
144
145 // calc J-field
146 if (!m_UseCellKappa)
147 {
148 m_DumpType = J_FIELD_DUMP;
149 field_td = CalcField();
150 T = m_Eng_Interface->GetTime(m_dualTime);
151 for (size_t n = 0; n<m_FD_Samples.size(); ++n)
152 {
153 std::complex<float> exp_jwt_2_dt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
154 exp_jwt_2_dt *= 2; // *2 for single-sided spectrum
155 exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval
156 field_fd = m_J_FD_Fields.at(n);
157 for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
158 {
159 for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
160 {
161 for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
162 {
163 field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
164 field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
165 field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
166 }
167 }
168 }
169 }
170 Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
171 }
172
173 //reset dump type
174 m_DumpType = save_dump_type;
175
176 ++m_FD_SampleCount;
177 return GetNextInterval();
178 }
179
DumpFDData()180 void ProcessFieldsSAR::DumpFDData()
181 {
182 if (Enabled==false) return;
183 unsigned int pos[3];
184 unsigned int orig_pos[3];
185 float*** SAR = Create3DArray<float>(numLines);
186 double coord[3];
187 ContinuousStructure* CSX = Op->GetGeometryCSX();
188 CSProperties* prop = NULL;
189 CSPropMaterial* matProp = NULL;
190
191 double power;
192
193 float*** cell_volume = Create3DArray<float>(numLines);
194 float*** cell_density = Create3DArray<float>(numLines);
195 float*** cell_kappa = NULL;
196 if (m_UseCellKappa)
197 cell_kappa = Create3DArray<float>(numLines);
198
199 bool found_UnIsotropic=false;
200
201 // calculate volumes and masses for all cells
202 for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
203 {
204 orig_pos[0] = posLines[0][pos[0]];
205 for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
206 {
207 orig_pos[1] = posLines[1][pos[1]];
208 vector<CSPrimitives*> vPrims = Op->GetPrimitivesBoundBox(orig_pos[0], orig_pos[1], -1, CSProperties::MATERIAL);
209 for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
210 {
211 orig_pos[2] = posLines[2][pos[2]];
212
213 cell_volume[pos[0]][pos[1]][pos[2]] = Op->GetCellVolume(orig_pos);
214 cell_density[pos[0]][pos[1]][pos[2]] = 0.0;
215
216 Op->GetCellCenterMaterialAvgCoord(orig_pos, coord);
217 prop = CSX->GetPropertyByCoordPriority(coord, vPrims);
218 // prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
219 if (prop!=0)
220 {
221 matProp = dynamic_cast<CSPropMaterial*>(prop);
222 if (matProp)
223 {
224 found_UnIsotropic |= !matProp->GetIsotropy();
225 cell_density[pos[0]][pos[1]][pos[2]] = matProp->GetDensityWeighted(coord);
226 if (m_UseCellKappa)
227 cell_kappa[pos[0]][pos[1]][pos[2]] = matProp->GetKappaWeighted(0,coord);
228 }
229 }
230 }
231 }
232 }
233 if (found_UnIsotropic)
234 cerr << "ProcessFieldsSAR::DumpFDData(): Warning, found unisotropic material in SAR calculation... this is unsupported!" << endl;
235
236 float* cellWidth[3];
237 for (int n=0;n<3;++n)
238 {
239 cellWidth[n]=new float[numLines[n]];
240 for (unsigned int i=0;i<numLines[n];++i)
241 cellWidth[n][i]=Op->GetDiscDelta(n,posLines[n][i])*Op->GetGridDelta();
242 }
243
244 if (m_DumpType == SAR_RAW_DATA)
245 {
246 if (m_fileType!=HDF5_FILETYPE)
247 {
248 cerr << "ProcessFieldsSAR::DumpFDData(): Error, wrong file type, this should not happen!!! skipped" << endl;
249 return;
250 }
251
252 size_t datasize[]={numLines[0],numLines[1],numLines[2]};
253 for (size_t n = 0; n<m_FD_Samples.size(); ++n)
254 {
255 stringstream ss;
256 ss << "f" << n;
257 if (m_HDF5_Dump_File->WriteVectorField(ss.str(), m_E_FD_Fields.at(n), datasize)==false)
258 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
259 }
260
261 m_HDF5_Dump_File->SetCurrentGroup("/CellData");
262 if (m_UseCellKappa==false)
263 cerr << "ProcessFieldsSAR::DumpFDData: Error, cell conductivity data not available, this should not happen... skipping! " << endl;
264 else if (m_HDF5_Dump_File->WriteScalarField("Conductivity", cell_kappa, datasize)==false)
265 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
266 if (m_HDF5_Dump_File->WriteScalarField("Density", cell_density, datasize)==false)
267 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
268 if (m_HDF5_Dump_File->WriteScalarField("Volume", cell_volume, datasize)==false)
269 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
270 }
271 else
272 {
273 SAR_Calculation SAR_Calc;
274 SAR_Calc.SetAveragingMethod(m_SAR_method, g_settings.GetVerboseLevel()==0);
275 SAR_Calc.SetDebugLevel(g_settings.GetVerboseLevel());
276 SAR_Calc.SetNumLines(numLines);
277 if (m_DumpType == SAR_LOCAL_DUMP)
278 SAR_Calc.SetAveragingMass(0);
279 else if (m_DumpType == SAR_1G_DUMP)
280 SAR_Calc.SetAveragingMass(1e-3);
281 else if (m_DumpType == SAR_10G_DUMP)
282 SAR_Calc.SetAveragingMass(10e-3);
283 else
284 {
285 cerr << "ProcessFieldsSAR::DumpFDData: unknown SAR dump type...!" << endl;
286 }
287 SAR_Calc.SetCellDensities(cell_density);
288 SAR_Calc.SetCellWidth(cellWidth);
289 SAR_Calc.SetCellVolumes(cell_volume);
290 SAR_Calc.SetCellCondictivity(cell_kappa); // cell_kappa will be NULL if m_UseCellKappa is false
291
292 for (size_t n = 0; n<m_FD_Samples.size(); ++n)
293 {
294 SAR_Calc.SetEField(m_E_FD_Fields.at(n));
295 if (!m_UseCellKappa)
296 SAR_Calc.SetJField(m_J_FD_Fields.at(n));
297 power = SAR_Calc.CalcSARPower();
298 SAR_Calc.CalcSAR(SAR);
299
300 if (m_fileType==VTK_FILETYPE)
301 {
302 stringstream ss;
303 ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n);
304
305 m_Vtk_Dump_File->SetFilename(ss.str());
306 m_Vtk_Dump_File->ClearAllFields();
307 m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR);
308 if (m_Vtk_Dump_File->Write()==false)
309 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
310 }
311 else if (m_fileType==HDF5_FILETYPE)
312 {
313 stringstream ss;
314 ss << "f" << n;
315 size_t datasize[]={numLines[0],numLines[1],numLines[2]};
316 if (m_HDF5_Dump_File->WriteScalarField(ss.str(), SAR, datasize)==false)
317 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
318 float freq[1] = {(float)m_FD_Samples.at(n)};
319 if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false)
320 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
321 float pow[1] = {(float)power};
322 if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false)
323 cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
324 }
325 else
326 cerr << "ProcessFieldsSAR::DumpFDData: unknown File-Type" << endl;
327 }
328 }
329 for (int n=0;n<3;++n)
330 delete[] cellWidth[n];
331 Delete3DArray(cell_volume,numLines);
332 Delete3DArray(cell_density,numLines);
333 Delete3DArray(cell_kappa,numLines);
334 Delete3DArray(SAR,numLines);
335 }
336