40 #include <boost/core/ignore_unused.hpp>
55 "Computes shear stress metrics.");
62 "First fld file. First underscore flags position of id in name.");
74 if (
m_f->m_exp[0]->GetNumElmts() == 0)
80 "Need to specify fromfld=file.fld ");
82 int nstart, i, j, nfields = 0;
85 string fromfld, basename, endname, nstartStr;
86 stringstream filename;
87 vector<string> infiles(nfld);
88 vector<std::shared_ptr<Field>> fromField(nfld);
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());
97 filename >> nstartStr;
99 endname = fromfld.substr(fromfld.find(nstartStr) + nstartStr.size(),
102 for (i = 0; i < nfld; ++i)
104 stringstream filename;
105 filename << basename << i + nstart << endname;
106 filename >> infiles[i];
107 cout << infiles[i] << endl;
110 for (i = 0; i < nfld; ++i)
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;
118 for (i = 0; i < nfld; ++i)
120 if (
m_f->m_exp.size())
124 for (j = 0; j <
m_f->m_exp[0]->GetExpSize(); ++j)
127 m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
129 m_f->FieldIOForFile(infiles[i])
130 ->Import(infiles[i], fromField[i]->m_fielddef,
131 fromField[i]->m_data,
136 m_f->FieldIOForFile(infiles[i])
137 ->Import(infiles[i], fromField[i]->m_fielddef,
138 fromField[i]->m_data,
142 nfields = fromField[i]->m_fielddef[0]->m_fields.size();
143 int NumHomogeneousDir =
144 fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
146 if (nfields == 5 || nfields == 7)
152 fromField[i]->m_graph->SetExpansionInfo(fromField[i]->m_fielddef);
155 fromField[i]->m_exp.resize(nfields);
156 fromField[i]->m_exp[0] =
157 fromField[i]->SetUpFirstExpList(NumHomogeneousDir,
true);
159 for (j = 1; j < nfields; ++j)
161 fromField[i]->m_exp[j] =
m_f->AppendExpList(NumHomogeneousDir);
164 for (j = 0; j < nfields; ++j)
166 for (
int k = 0; k < fromField[i]->m_data.size(); ++k)
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());
173 fromField[i]->m_exp[j]->BwdTrans(
174 fromField[i]->m_exp[j]->GetCoeffs(),
175 fromField[i]->m_exp[j]->UpdatePhys());
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)
192 int npoints = fromField[0]->m_exp[0]->GetNpoints();
194 normCrossDir(spacedim), outfield(nout), dtemp(spacedim);
196 DotProduct(npoints, 0.0), temp(npoints, 0.0);
198 for (i = 0; i < spacedim; ++i)
207 for (i = 0; i < nout; ++i)
215 for (i = 0; i < nfld; ++i)
217 for (j = 0; j < spacedim; ++j)
219 Vmath::Vadd(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
220 normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
224 for (i = 0; i < spacedim; ++i)
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);
232 Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
234 for (i = 0; i < spacedim; ++i)
236 Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1,
237 normTemporalMeanVec[i], 1);
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,
250 Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
251 normTemporalMeanVec[2], 1, normCrossDir[1], 1);
252 Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
253 normTemporalMeanVec[0], 1, normCrossDir[1], 1,
256 Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
257 normTemporalMeanVec[0], 1, normCrossDir[2], 1);
258 Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
259 normTemporalMeanVec[1], 1, normCrossDir[2], 1,
264 for (i = 0; i < nfld; ++i)
266 for (j = 0; j < spacedim; ++j)
268 Vmath::Vvtvp(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
269 normTemporalMeanVec[j], 1, DotProduct, 1, DotProduct,
274 Vmath::Vadd(npoints, fromField[i]->m_exp[spacedim]->GetPhys(), 1,
275 outfield[0], 1, outfield[0], 1);
278 Vmath::Vmul(npoints, DotProduct, 1, DotProduct, 1, temp, 1);
279 Vmath::Vvtvm(npoints, fromField[i]->m_exp[spacedim]->GetPhys(), 1,
280 fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1, temp,
283 for (j = 0; j < npoints; ++j)
287 outfield[1][j] = outfield[1][j] +
sqrt(temp[j]);
293 fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
294 Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
297 for (j = 0; j < npoints; ++j)
299 temp[j] = 1 - temp[j] * temp[j];
302 outfield[4][j] = outfield[4][j] +
sqrt(temp[j]);
312 fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
314 for (j = 0; j < spacedim; j++)
316 Vmath::Vvtvp(npoints, dtemp[j], 1, normTemporalMeanVec[j], 1,
324 for (j = 0; j < spacedim; ++j)
326 Vmath::Vvtvp(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
327 normCrossDir[j], 1, DotProduct, 1, DotProduct, 1);
329 fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
333 for (j = 0; j < spacedim; j++)
336 DotProduct, 1, DotProduct, 1);
338 Vmath::Vvtvp(npoints, DotProduct, 1, DotProduct, 1, temp, 1, temp,
343 Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
350 Vmath::Smul(npoints, 1.0 / nfld, outfield[0], 1, outfield[0], 1);
351 Vmath::Smul(npoints, 1.0 / nfld, outfield[1], 1, outfield[1], 1);
352 Vmath::Smul(npoints, 1.0 / nfld, outfield[3], 1, outfield[3], 1);
353 Vmath::Smul(npoints, 1.0 / nfld, outfield[4], 1, outfield[4], 1);
356 for (i = 0; i < npoints; ++i)
358 outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i] / outfield[0][i]);
368 m_f->m_exp.resize(nout);
369 m_f->m_fielddef = fromField[0]->m_fielddef;
371 m_f->SetUpFirstExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir,
true);
373 for (i = 1; i < nout; ++i)
376 m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
379 m_f->m_fielddef[0]->m_fields.resize(nout);
380 m_f->m_fielddef[0]->m_fields[0] =
"TAWSS";
381 m_f->m_fielddef[0]->m_fields[1] =
"transWSS";
382 m_f->m_fielddef[0]->m_fields[2] =
"OSI";
383 m_f->m_fielddef[0]->m_fields[3] =
"TAAFI";
384 m_f->m_fielddef[0]->m_fields[4] =
"TACFI";
388 m_f->m_fielddef[0]->m_fields[5] =
"|WSSG|";
389 Vmath::Smul(npoints, 1.0 / nfld, outfield[5], 1, outfield[5], 1);
392 m_f->m_variables =
m_f->m_fielddef[0]->m_fields;
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());
#define ASSERTL0(condition, msg)
FieldSharedPtr m_f
Field object.
std::map< std::string, ConfigOption > m_config
List of configuration values.
Abstract base class for processing modules.
virtual ~ProcessMultiShear()
virtual void Process(po::variables_map &vm)
Write mesh to output file.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
std::shared_ptr< Field > FieldSharedPtr
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
static FieldMetaDataMap NullFieldMetaDataMap
The above copyright notice and this permission notice shall be included.
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 plus 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)
Represents a command-line configuration option.