Write mesh to output file.
79 if (
m_f->m_comm->TreatAsRankZero())
81 cout <<
"ProcessMultiShear: Calculating shear stress metrics..."
86 int nstart, i, j, nfields=0;
89 string fromfld, basename, endname, nstartStr;
90 stringstream filename;
91 vector<string> infiles(nfld);
92 vector<boost::shared_ptr<Field> > m_fromField(nfld);
95 fromfld =
m_config[
"fromfld"].as<
string>();
96 basename = fromfld.substr(0, fromfld.find_first_of(
"_") + 1);
97 filename << fromfld.substr(fromfld.find_first_of(
"_") + 1, fromfld.size());
101 filename >> nstartStr;
103 endname = fromfld.substr(fromfld.find(nstartStr) + nstartStr.size(),
106 for (i = 0; i < nfld; ++i)
108 stringstream filename;
109 filename << basename << i + nstart << endname;
110 filename >> infiles[i];
111 cout << infiles[i] << endl;
114 for (i = 0; i < nfld; ++i)
116 m_fromField[i] = boost::shared_ptr<Field>(
new Field());
117 m_fromField[i]->m_session =
m_f->m_session;
118 m_fromField[i]->m_graph =
m_f->m_graph;
122 for (i = 0; i < nfld; ++i)
124 if (
m_f->m_exp.size())
127 Array<OneD, int> ElementGIDs(
m_f->m_exp[0]->GetExpSize());
128 for (j = 0; j <
m_f->m_exp[0]->GetExpSize(); ++j)
131 m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
133 m_f->FieldIOForFile(infiles[i])->Import(
134 infiles[i], m_fromField[i]->m_fielddef, m_fromField[i]->m_data,
139 m_f->FieldIOForFile(infiles[i])->Import(
140 infiles[i], m_fromField[i]->m_fielddef, m_fromField[i]->m_data,
144 nfields = m_fromField[i]->m_fielddef[0]->m_fields.size();
145 int NumHomogeneousDir =
146 m_fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
148 if (nfields == 5 || nfields == 7)
154 m_fromField[i]->m_graph->SetExpansions(m_fromField[i]->m_fielddef);
157 m_fromField[i]->m_exp.resize(nfields);
158 m_fromField[i]->m_exp[0] =
159 m_fromField[i]->SetUpFirstExpList(NumHomogeneousDir,
true);
161 for (j = 1; j < nfields; ++j)
163 m_fromField[i]->m_exp[j] =
m_f->AppendExpList(NumHomogeneousDir);
166 for (j = 0; j < nfields; ++j)
168 for (
int k = 0; k < m_fromField[i]->m_data.size(); ++k)
170 m_fromField[i]->m_exp[j]->ExtractDataToCoeffs(
171 m_fromField[i]->m_fielddef[k], m_fromField[i]->m_data[k],
172 m_fromField[i]->m_fielddef[k]->m_fields[j],
173 m_fromField[i]->m_exp[j]->UpdateCoeffs());
175 m_fromField[i]->m_exp[j]->BwdTrans(
176 m_fromField[i]->m_exp[j]->GetCoeffs(),
177 m_fromField[i]->m_exp[j]->UpdatePhys());
181 int spacedim =
m_f->m_graph->GetSpaceDimension();
182 if ((m_fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
183 (m_fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 2)
194 int npoints = m_fromField[0]->m_exp[0]->GetNpoints();
195 Array<OneD, Array<OneD, NekDouble> > normTemporalMeanVec(spacedim),
196 normCrossDir(spacedim), outfield(nout), dtemp(spacedim);
197 Array<OneD, NekDouble> TemporalMeanMag(npoints, 0.0),
198 DotProduct(npoints, 0.0), temp(npoints, 0.0);
200 for (i = 0; i < spacedim; ++i)
202 normTemporalMeanVec[i] = Array<OneD, NekDouble>(npoints);
203 normCrossDir[i] = Array<OneD, NekDouble>(npoints);
204 dtemp[i] = Array<OneD, NekDouble>(npoints);
209 for (i = 0; i < nout; ++i)
211 outfield[i] = Array<OneD, NekDouble>(npoints, 0.0);
217 for (i = 0; i < nfld; ++i)
219 for (j = 0; j < spacedim; ++j)
221 Vmath::Vadd(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1,
222 normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
226 for (i = 0; i < spacedim; ++i)
228 Vmath::Smul(npoints, 1.0 / nfld, normTemporalMeanVec[i], 1,
229 normTemporalMeanVec[i], 1);
230 Vmath::Vvtvp(npoints, normTemporalMeanVec[i], 1, normTemporalMeanVec[i],
231 1, TemporalMeanMag, 1, TemporalMeanMag, 1);
234 Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
236 for (i = 0; i < spacedim; ++i)
238 Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1,
239 normTemporalMeanVec[i], 1);
246 Vmath::Vmul(npoints, m_fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
247 normTemporalMeanVec[1], 1, normCrossDir[0], 1);
248 Vmath::Vvtvm(npoints, m_fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
249 normTemporalMeanVec[2], 1, normCrossDir[0], 1,
252 Vmath::Vmul(npoints, m_fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
253 normTemporalMeanVec[2], 1, normCrossDir[1], 1);
254 Vmath::Vvtvm(npoints, m_fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
255 normTemporalMeanVec[0], 1, normCrossDir[1], 1,
258 Vmath::Vmul(npoints, m_fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
259 normTemporalMeanVec[0], 1, normCrossDir[2], 1);
260 Vmath::Vvtvm(npoints, m_fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
261 normTemporalMeanVec[1], 1, normCrossDir[2], 1,
266 for (i = 0; i < nfld; ++i)
268 for (j = 0; j < spacedim; ++j)
270 Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1,
271 normTemporalMeanVec[j], 1, DotProduct, 1, DotProduct,
276 Vmath::Vadd(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1,
277 outfield[0], 1, outfield[0], 1);
280 Vmath::Vmul(npoints, DotProduct, 1, DotProduct, 1, temp, 1);
281 Vmath::Vvtvm(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1,
282 m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1,
285 for (j = 0; j < npoints; ++j)
289 outfield[1][j] = outfield[1][j] + sqrt(temp[j]);
295 m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
296 Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
299 for (j = 0; j < npoints; ++j)
301 temp[j] = 1 - temp[j] * temp[j];
304 outfield[4][j] = outfield[4][j] + sqrt(temp[j]);
314 m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
316 for (j = 0; j < spacedim; j++)
318 Vmath::Vvtvp(npoints, dtemp[j], 1, normTemporalMeanVec[j], 1,
326 for (j = 0; j < spacedim; ++j)
328 Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1,
329 normCrossDir[j], 1, DotProduct, 1, DotProduct, 1);
331 m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
335 for (j = 0; j < spacedim; j++)
338 DotProduct, 1, DotProduct, 1);
340 Vmath::Vvtvp(npoints, DotProduct, 1, DotProduct, 1, temp, 1, temp,
345 Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
352 Vmath::Smul(npoints, 1.0 / nfld, outfield[0], 1, outfield[0], 1);
353 Vmath::Smul(npoints, 1.0 / nfld, outfield[1], 1, outfield[1], 1);
354 Vmath::Smul(npoints, 1.0 / nfld, outfield[3], 1, outfield[3], 1);
355 Vmath::Smul(npoints, 1.0 / nfld, outfield[4], 1, outfield[4], 1);
358 for (i = 0; i < npoints; ++i)
360 outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i] / outfield[0][i]);
370 m_f->m_exp.resize(nout);
371 m_f->m_fielddef = m_fromField[0]->m_fielddef;
373 m_f->SetUpFirstExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir,
true);
375 for (i = 1; i < nout; ++i)
378 m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
381 m_f->m_fielddef[0]->m_fields.resize(nout);
382 m_f->m_fielddef[0]->m_fields[0] =
"TAWSS";
383 m_f->m_fielddef[0]->m_fields[1] =
"transWSS";
384 m_f->m_fielddef[0]->m_fields[2] =
"OSI";
385 m_f->m_fielddef[0]->m_fields[3] =
"TAAFI";
386 m_f->m_fielddef[0]->m_fields[4] =
"TACFI";
390 m_f->m_fielddef[0]->m_fields[5] =
"|WSSG|";
391 Vmath::Smul(npoints, 1.0 / nfld, outfield[5], 1, outfield[5], 1);
394 for (i = 0; i < nout; ++i)
396 m_f->m_exp[i]->FwdTrans(outfield[i],
m_f->m_exp[i]->UpdateCoeffs());
397 m_f->m_exp[i]->BwdTrans(
m_f->m_exp[i]->GetCoeffs(),
398 m_f->m_exp[i]->UpdatePhys());
401 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
402 m_fromField[0]->m_exp[0]->GetFieldDefinitions();
403 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
405 for (i = 0; i < nout; ++i)
407 for (j = 0; j < FieldDef.size(); ++j)
409 FieldDef[j]->m_fields.push_back(
m_f->m_fielddef[0]->m_fields[i]);
410 m_f->m_exp[i]->AppendFieldData(FieldDef[j], FieldData[j]);
414 m_f->m_fielddef = FieldDef;
415 m_f->m_data = FieldData;
map< string, ConfigOption > m_config
List of configuration values.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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 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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*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 plus vector): z = w*x - y
void Zero(int n, T *x, const int incx)
Zero vector.
static FieldMetaDataMap NullFieldMetaDataMap
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 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.
FieldSharedPtr m_f
Field object.