1 /*
2 * Copyright (C) 2010 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 <iomanip>
19 #include "tools/global.h"
20 #include "tools/vtk_file_writer.h"
21 #include "tools/hdf5_file_writer.h"
22 #include "processfields.h"
23 #include "FDTD/engine_interface_fdtd.h"
24
ProcessFields(Engine_Interface_Base * eng_if)25 ProcessFields::ProcessFields(Engine_Interface_Base* eng_if) : Processing(eng_if)
26 {
27 m_DumpType = E_FIELD_DUMP;
28 // vtk-file is default
29 m_fileType = VTK_FILETYPE;
30 m_SampleType = NONE;
31 m_Vtk_Dump_File = NULL;
32 m_HDF5_Dump_File = NULL;
33 SetPrecision(6);
34 m_dualTime = false;
35
36 // dump box should be always inside the snapped lines
37 m_SnapMethod = 1;
38
39 for (int n=0; n<3; ++n)
40 {
41 numLines[n]=0;
42 posLines[n]=NULL;
43 discLines[n]=NULL;
44 subSample[n]=1;
45 optResolution[n]=0;
46 }
47 }
48
~ProcessFields()49 ProcessFields::~ProcessFields()
50 {
51 delete m_Vtk_Dump_File;
52 m_Vtk_Dump_File = NULL;
53 for (int n=0; n<3; ++n)
54 {
55 delete[] posLines[n];
56 posLines[n]=NULL;
57 delete[] discLines[n];
58 discLines[n]=NULL;
59 }
60 }
61
GetFieldNameByType(DumpType type)62 string ProcessFields::GetFieldNameByType(DumpType type)
63 {
64 switch (type)
65 {
66 case E_FIELD_DUMP:
67 return "E-Field";
68 case H_FIELD_DUMP:
69 return "H-Field";
70 case J_FIELD_DUMP:
71 return "J-Field";
72 case ROTH_FIELD_DUMP:
73 return "RotH-Field";
74 case D_FIELD_DUMP:
75 return "D-Field";
76 case B_FIELD_DUMP:
77 return "B-Field";
78 case SAR_LOCAL_DUMP:
79 return "SAR-local";
80 case SAR_1G_DUMP:
81 return "SAR_1g";
82 case SAR_10G_DUMP:
83 return "SAR_10g";
84 case SAR_RAW_DATA:
85 return "SAR_raw_data";
86 }
87 return "unknown field";
88 }
89
NeedConductivity() const90 bool ProcessFields::NeedConductivity() const
91 {
92 switch (m_DumpType)
93 {
94 case J_FIELD_DUMP:
95 return true;
96 default:
97 return false;
98 }
99 return false;
100 }
101
NeedPermittivity() const102 bool ProcessFields::NeedPermittivity() const
103 {
104 switch (m_DumpType)
105 {
106 case D_FIELD_DUMP:
107 return true;
108 default:
109 return false;
110 }
111 return false;
112 }
113
NeedPermeability() const114 bool ProcessFields::NeedPermeability() const
115 {
116 switch (m_DumpType)
117 {
118 case B_FIELD_DUMP:
119 return true;
120 default:
121 return false;
122 }
123 return false;
124 }
125
InitProcess()126 void ProcessFields::InitProcess()
127 {
128 if (Enabled==false) return;
129
130 CalcMeshPos();
131
132 if (m_fileType==VTK_FILETYPE)
133 {
134 delete m_Vtk_Dump_File;
135 m_Vtk_Dump_File = new VTK_File_Writer(m_filename,(int)m_Mesh_Type);
136
137 #ifdef OUTPUT_IN_DRAWINGUNITS
138 double discScaling = 1;
139 #else
140 double discScaling = Op->GetGridDelta();
141 #endif
142 m_Vtk_Dump_File->SetMeshLines(discLines,numLines,discScaling);
143 m_Vtk_Dump_File->SetNativeDump(g_settings.NativeFieldDumps());
144 }
145 if (m_fileType==HDF5_FILETYPE)
146 {
147 delete m_HDF5_Dump_File;
148 m_HDF5_Dump_File = new HDF5_File_Writer(m_filename+".h5");
149
150 #ifdef OUTPUT_IN_DRAWINGUNITS
151 double discScaling = 1;
152 #else
153 double discScaling = Op->GetGridDelta();
154 #endif
155 m_HDF5_Dump_File->WriteRectMesh(numLines,discLines,(int)m_Mesh_Type,discScaling);
156
157 m_HDF5_Dump_File->WriteAtrribute("/","openEMS_HDF5_version",0.2);
158 }
159 }
160
SetDumpMode(Engine_Interface_Base::InterpolationType mode)161 void ProcessFields::SetDumpMode(Engine_Interface_Base::InterpolationType mode)
162 {
163 m_Eng_Interface->SetInterpolationType(mode);
164 if (mode==Engine_Interface_Base::CELL_INTERPOLATE)
165 m_dualMesh=true;
166 else if (mode==Engine_Interface_Base::NODE_INTERPOLATE)
167 m_dualMesh=false;
168 //else keep the preset/user defined case
169 }
170
DefineStartStopCoord(double * dstart,double * dstop)171 void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
172 {
173 Processing::DefineStartStopCoord(dstart,dstop);
174
175 // normalize order of start and stop
176 for (int n=0; n<3; ++n)
177 {
178 if (start[n]>stop[n])
179 {
180 unsigned int help = start[n];
181 start[n]=stop[n];
182 stop[n]=help;
183 }
184 }
185 }
186
CalcTotalEnergyEstimate() const187 double ProcessFields::CalcTotalEnergyEstimate() const
188 {
189 return m_Eng_Interface->CalcFastEnergy();
190 }
191
SetSubSampling(unsigned int subSampleRate,int dir)192 void ProcessFields::SetSubSampling(unsigned int subSampleRate, int dir)
193 {
194 if (dir>2) return;
195 if (dir<0)
196 {
197 subSample[0]=subSampleRate;
198 subSample[1]=subSampleRate;
199 subSample[2]=subSampleRate;
200 }
201 else subSample[dir]=subSampleRate;
202 m_SampleType = SUBSAMPLE;
203 }
204
SetOptResolution(double optRes,int dir)205 void ProcessFields::SetOptResolution(double optRes, int dir)
206 {
207 if (dir>2) return;
208 if (dir<0)
209 {
210 optResolution[0]=optRes;
211 optResolution[1]=optRes;
212 optResolution[2]=optRes;
213 }
214 else optResolution[dir]=optRes;
215 m_SampleType = OPT_RESOLUTION;
216 }
217
CalcMeshPos()218 void ProcessFields::CalcMeshPos()
219 {
220 if ((m_SampleType==SUBSAMPLE) || (m_SampleType==NONE))
221 {
222 vector<unsigned int> tmp_pos;
223
224 for (int n=0; n<3; ++n)
225 {
226 // construct new discLines
227 tmp_pos.clear();
228 for (unsigned int i=start[n]; i<=stop[n]; i+=subSample[n])
229 tmp_pos.push_back(i);
230
231 numLines[n] = tmp_pos.size();
232 delete[] discLines[n];
233 discLines[n] = new double[numLines[n]];
234 delete[] posLines[n];
235 posLines[n] = new unsigned int[numLines[n]];
236 for (unsigned int i=0; i<numLines[n]; ++i)
237 {
238 posLines[n][i] = tmp_pos.at(i);
239 discLines[n][i] = Op->GetDiscLine(n,tmp_pos.at(i),m_dualMesh);
240 }
241 }
242 }
243 if ((m_SampleType==OPT_RESOLUTION))
244 {
245 vector<unsigned int> tmp_pos;
246 double oldPos=0;
247 for (int n=0; n<3; ++n)
248 {
249 // construct new discLines
250 tmp_pos.clear();
251 tmp_pos.push_back(start[n]);
252 oldPos=Op->GetDiscLine(n,start[n],m_dualMesh);
253 if (stop[n]==0)
254 tmp_pos.push_back(stop[n]);
255 else
256 for (unsigned int i=start[n]+1; i<=stop[n]-1; ++i)
257 {
258 if ( (Op->GetDiscLine(n,i+1,m_dualMesh)-oldPos) >= optResolution[n])
259 {
260 tmp_pos.push_back(i);
261 oldPos=Op->GetDiscLine(n,i,m_dualMesh);
262 }
263 }
264 if (start[n]!=stop[n])
265 tmp_pos.push_back(stop[n]);
266 numLines[n] = tmp_pos.size();
267 delete[] discLines[n];
268 discLines[n] = new double[numLines[n]];
269 delete[] posLines[n];
270 posLines[n] = new unsigned int[numLines[n]];
271 for (unsigned int i=0; i<numLines[n]; ++i)
272 {
273 posLines[n][i] = tmp_pos.at(i);
274 discLines[n][i] = Op->GetDiscLine(n,tmp_pos.at(i),m_dualMesh);
275 }
276 }
277 }
278 }
279
CalcField()280 FDTD_FLOAT**** ProcessFields::CalcField()
281 {
282 unsigned int pos[3];
283 double out[3];
284 //create array
285 FDTD_FLOAT**** field = Create_N_3DArray<FDTD_FLOAT>(numLines);
286 switch (m_DumpType)
287 {
288 case E_FIELD_DUMP:
289 for (unsigned int i=0; i<numLines[0]; ++i)
290 {
291 pos[0]=posLines[0][i];
292 for (unsigned int j=0; j<numLines[1]; ++j)
293 {
294 pos[1]=posLines[1][j];
295 for (unsigned int k=0; k<numLines[2]; ++k)
296 {
297 pos[2]=posLines[2][k];
298
299 m_Eng_Interface->GetEField(pos,out);
300 field[0][i][j][k] = out[0];
301 field[1][i][j][k] = out[1];
302 field[2][i][j][k] = out[2];
303 }
304 }
305 }
306 return field;
307 case H_FIELD_DUMP:
308 for (unsigned int i=0; i<numLines[0]; ++i)
309 {
310 pos[0]=posLines[0][i];
311 for (unsigned int j=0; j<numLines[1]; ++j)
312 {
313 pos[1]=posLines[1][j];
314 for (unsigned int k=0; k<numLines[2]; ++k)
315 {
316 pos[2]=posLines[2][k];
317
318 m_Eng_Interface->GetHField(pos,out);
319 field[0][i][j][k] = out[0];
320 field[1][i][j][k] = out[1];
321 field[2][i][j][k] = out[2];
322 }
323 }
324 }
325 return field;
326 case J_FIELD_DUMP:
327 for (unsigned int i=0; i<numLines[0]; ++i)
328 {
329 pos[0]=posLines[0][i];
330 for (unsigned int j=0; j<numLines[1]; ++j)
331 {
332 pos[1]=posLines[1][j];
333 for (unsigned int k=0; k<numLines[2]; ++k)
334 {
335 pos[2]=posLines[2][k];
336
337 m_Eng_Interface->GetJField(pos,out);
338 field[0][i][j][k] = out[0];
339 field[1][i][j][k] = out[1];
340 field[2][i][j][k] = out[2];
341 }
342 }
343 }
344 return field;
345 case ROTH_FIELD_DUMP:
346 for (unsigned int i=0; i<numLines[0]; ++i)
347 {
348 pos[0]=posLines[0][i];
349 for (unsigned int j=0; j<numLines[1]; ++j)
350 {
351 pos[1]=posLines[1][j];
352 for (unsigned int k=0; k<numLines[2]; ++k)
353 {
354 pos[2]=posLines[2][k];
355
356 m_Eng_Interface->GetRotHField(pos,out);
357 field[0][i][j][k] = out[0];
358 field[1][i][j][k] = out[1];
359 field[2][i][j][k] = out[2];
360 }
361 }
362 }
363 return field;
364 case D_FIELD_DUMP:
365 for (unsigned int i=0; i<numLines[0]; ++i)
366 {
367 pos[0]=posLines[0][i];
368 for (unsigned int j=0; j<numLines[1]; ++j)
369 {
370 pos[1]=posLines[1][j];
371 for (unsigned int k=0; k<numLines[2]; ++k)
372 {
373 pos[2]=posLines[2][k];
374
375 m_Eng_Interface->GetDField(pos,out);
376 field[0][i][j][k] = out[0];
377 field[1][i][j][k] = out[1];
378 field[2][i][j][k] = out[2];
379 }
380 }
381 }
382 return field;
383 case B_FIELD_DUMP:
384 for (unsigned int i=0; i<numLines[0]; ++i)
385 {
386 pos[0]=posLines[0][i];
387 for (unsigned int j=0; j<numLines[1]; ++j)
388 {
389 pos[1]=posLines[1][j];
390 for (unsigned int k=0; k<numLines[2]; ++k)
391 {
392 pos[2]=posLines[2][k];
393
394 m_Eng_Interface->GetBField(pos,out);
395 field[0][i][j][k] = out[0];
396 field[1][i][j][k] = out[1];
397 field[2][i][j][k] = out[2];
398 }
399 }
400 }
401 return field;
402 default:
403 cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;
404 return field;
405 }
406 }
407
408