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