Write mesh to output file.
70{
72
73
74 if (
m_f->m_exp[0]->GetNumElmts() == 0)
75 {
76 return;
77 }
78
80 "Need to specify fromfld=file.fld ");
81
82 int nstart, i, j, nfields = 0;
83 bool wssg = false;
85 string fromfld, basename, endname, nstartStr;
86 stringstream filename;
87 vector<string> infiles(nfld);
88 vector<std::shared_ptr<Field>> fromField(nfld);
89
90
91 fromfld =
m_config[
"fromfld"].as<
string>();
92 basename = fromfld.substr(0, fromfld.find_first_of("_") + 1);
93 filename << fromfld.substr(fromfld.find_first_of("_") + 1, fromfld.size());
94 filename >> nstart;
95 filename.str("");
96 filename << nstart;
97 filename >> nstartStr;
98 filename.str("");
99 endname = fromfld.substr(fromfld.find(nstartStr) + nstartStr.size(),
100 fromfld.size());
101
102 for (i = 0; i < nfld; ++i)
103 {
104 stringstream filename;
105 filename << basename << i + nstart << endname;
106 filename >> infiles[i];
107 cout << infiles[i] << endl;
108 }
109
110 for (i = 0; i < nfld; ++i)
111 {
112 fromField[i] = std::shared_ptr<Field>(new Field());
113 fromField[i]->m_session =
m_f->m_session;
114 fromField[i]->m_graph =
m_f->m_graph;
115 }
116
117
118 for (i = 0; i < nfld; ++i)
119 {
120 if (
m_f->m_exp.size())
121 {
122
123 Array<OneD, int> ElementGIDs(
m_f->m_exp[0]->GetExpSize());
124 for (j = 0; j <
m_f->m_exp[0]->GetExpSize(); ++j)
125 {
126 ElementGIDs[j] =
127 m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
128 }
129 m_f->FieldIOForFile(infiles[i])
130 ->Import(infiles[i], fromField[i]->m_fielddef,
131 fromField[i]->m_data,
133 }
134 else
135 {
136 m_f->FieldIOForFile(infiles[i])
137 ->Import(infiles[i], fromField[i]->m_fielddef,
138 fromField[i]->m_data,
140 }
141
142 nfields = fromField[i]->m_fielddef[0]->m_fields.size();
143 int NumHomogeneousDir =
144 fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
145
146 if (nfields == 5 || nfields == 7)
147 {
148 wssg = true;
149 }
150
151
152 fromField[i]->m_graph->SetExpansionInfo(fromField[i]->m_fielddef);
153
154
155 fromField[i]->m_exp.resize(nfields);
156 fromField[i]->m_exp[0] =
157 fromField[i]->SetUpFirstExpList(NumHomogeneousDir, true);
158
159 for (j = 1; j < nfields; ++j)
160 {
161 fromField[i]->m_exp[j] =
m_f->AppendExpList(NumHomogeneousDir);
162 }
163
164 for (j = 0; j < nfields; ++j)
165 {
166 for (int k = 0; k < fromField[i]->m_data.size(); ++k)
167 {
168 fromField[i]->m_exp[j]->ExtractDataToCoeffs(
169 fromField[i]->m_fielddef[k], fromField[i]->m_data[k],
170 fromField[i]->m_fielddef[k]->m_fields[j],
171 fromField[i]->m_exp[j]->UpdateCoeffs());
172 }
173 fromField[i]->m_exp[j]->BwdTrans(
174 fromField[i]->m_exp[j]->GetCoeffs(),
175 fromField[i]->m_exp[j]->UpdatePhys());
176 }
177 }
178
179 int spacedim =
m_f->m_graph->GetSpaceDimension();
180 if ((fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
181 (fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 2)
182 {
183 spacedim = 3;
184 }
185
186 int nout = 5;
187 if (wssg)
188 {
189 nout = 6;
190 }
191
192 int npoints = fromField[0]->m_exp[0]->GetNpoints();
193 Array<OneD, Array<OneD, NekDouble>> normTemporalMeanVec(spacedim),
194 normCrossDir(spacedim), outfield(nout), dtemp(spacedim);
195 Array<OneD, NekDouble> TemporalMeanMag(npoints, 0.0),
196 DotProduct(npoints, 0.0), temp(npoints, 0.0);
197
198 for (i = 0; i < spacedim; ++i)
199 {
200 normTemporalMeanVec[i] = Array<OneD, NekDouble>(npoints);
201 normCrossDir[i] = Array<OneD, NekDouble>(npoints);
202 dtemp[i] = Array<OneD, NekDouble>(npoints);
205 }
206
207 for (i = 0; i < nout; ++i)
208 {
209 outfield[i] = Array<OneD, NekDouble>(npoints, 0.0);
210 }
211
212
213
214
215 for (i = 0; i < nfld; ++i)
216 {
217 for (j = 0; j < spacedim; ++j)
218 {
219 Vmath::Vadd(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
220 normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
221 }
222 }
223
224 for (i = 0; i < spacedim; ++i)
225 {
226 Vmath::Smul(npoints, 1.0 / nfld, normTemporalMeanVec[i], 1,
227 normTemporalMeanVec[i], 1);
228 Vmath::Vvtvp(npoints, normTemporalMeanVec[i], 1, normTemporalMeanVec[i],
229 1, TemporalMeanMag, 1, TemporalMeanMag, 1);
230 }
231
232 Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
233
234 for (i = 0; i < spacedim; ++i)
235 {
236 Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1,
237 normTemporalMeanVec[i], 1);
238 }
239
240
241 if (wssg)
242
243 {
244 Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
245 normTemporalMeanVec[1], 1, normCrossDir[0], 1);
246 Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
247 normTemporalMeanVec[2], 1, normCrossDir[0], 1,
248 normCrossDir[0], 1);
249 Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
250 normTemporalMeanVec[2], 1, normCrossDir[1], 1);
251 Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
252 normTemporalMeanVec[0], 1, normCrossDir[1], 1,
253 normCrossDir[1], 1);
254 Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
255 normTemporalMeanVec[0], 1, normCrossDir[2], 1);
256 Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
257 normTemporalMeanVec[1], 1, normCrossDir[2], 1,
258 normCrossDir[2], 1);
259 }
260
261
262 for (i = 0; i < nfld; ++i)
263 {
264 for (j = 0; j < spacedim; ++j)
265 {
266 Vmath::Vvtvp(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
267 normTemporalMeanVec[j], 1, DotProduct, 1, DotProduct,
268 1);
269 }
270
271 Vmath::Vadd(npoints, fromField[i]->m_exp[spacedim]->GetPhys(), 1,
272 outfield[0], 1, outfield[0], 1);
273
274 Vmath::Vmul(npoints, DotProduct, 1, DotProduct, 1, temp, 1);
275 Vmath::Vvtvm(npoints, fromField[i]->m_exp[spacedim]->GetPhys(), 1,
276 fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1, temp,
277 1);
278
279 for (j = 0; j < npoints; ++j)
280 {
281 if (temp[j] > 0.0)
282 {
283 outfield[1][j] = outfield[1][j] +
sqrt(temp[j]);
284 }
285 }
286
287
289 fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
290 Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
291
292
293 for (j = 0; j < npoints; ++j)
294 {
295 temp[j] = 1 - temp[j] * temp[j];
296 if (temp[j] > 0.0)
297 {
298 outfield[4][j] = outfield[4][j] +
sqrt(temp[j]);
299 }
300 }
301
302
303 if (wssg)
304 {
305
307
308 fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
309 dtemp[2]);
310 for (j = 0; j < spacedim; j++)
311 {
312 Vmath::Vvtvp(npoints, dtemp[j], 1, normTemporalMeanVec[j], 1,
313 temp, 1, temp, 1);
314 }
316
317
319
320 for (j = 0; j < spacedim; ++j)
321 {
322 Vmath::Vvtvp(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
323 normCrossDir[j], 1, DotProduct, 1, DotProduct, 1);
324 }
325 fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
326 dtemp[2]);
328
329 for (j = 0; j < spacedim; j++)
330 {
332 DotProduct, 1, DotProduct, 1);
333 }
334 Vmath::Vvtvp(npoints, DotProduct, 1, DotProduct, 1, temp, 1, temp,
335 1);
336
337
339 Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
340 }
341
343 }
344
345
346 Vmath::Smul(npoints, 1.0 / nfld, outfield[0], 1, outfield[0], 1);
347 Vmath::Smul(npoints, 1.0 / nfld, outfield[1], 1, outfield[1], 1);
348 Vmath::Smul(npoints, 1.0 / nfld, outfield[3], 1, outfield[3], 1);
349 Vmath::Smul(npoints, 1.0 / nfld, outfield[4], 1, outfield[4], 1);
350
351
352 for (i = 0; i < npoints; ++i)
353 {
354 outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i] / outfield[0][i]);
355 }
356
357
358
359
360
361
362
363
364 m_f->m_exp.resize(nout);
365 m_f->m_fielddef = fromField[0]->m_fielddef;
367 m_f->SetUpFirstExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir,
true);
368
369 for (i = 1; i < nout; ++i)
370 {
372 m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
373 }
374
375 m_f->m_fielddef[0]->m_fields.resize(nout);
376 m_f->m_fielddef[0]->m_fields[0] =
"TAWSS";
377 m_f->m_fielddef[0]->m_fields[1] =
"transWSS";
378 m_f->m_fielddef[0]->m_fields[2] =
"OSI";
379 m_f->m_fielddef[0]->m_fields[3] =
"TAAFI";
380 m_f->m_fielddef[0]->m_fields[4] =
"TACFI";
381
382 if (wssg)
383 {
384 m_f->m_fielddef[0]->m_fields[5] =
"|WSSG|";
385 Vmath::Smul(npoints, 1.0 / nfld, outfield[5], 1, outfield[5], 1);
386 }
387
388 m_f->m_variables =
m_f->m_fielddef[0]->m_fields;
389
390 for (i = 0; i < nout; ++i)
391 {
392 m_f->m_exp[i]->FwdTrans(outfield[i],
m_f->m_exp[i]->UpdateCoeffs());
393 m_f->m_exp[i]->BwdTrans(
m_f->m_exp[i]->GetCoeffs(),
394 m_f->m_exp[i]->UpdatePhys());
395 }
396}
#define ASSERTL0(condition, msg)
FieldSharedPtr m_f
Field object.
static FieldMetaDataMap NullFieldMetaDataMap
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
void Zero(int n, T *x, const int incx)
Zero vector.
scalarT< T > sqrt(scalarT< T > in)