55 ProcessMultiShear::create,
"Computes shear stress metrics.");
61 "First fld file. First underscore flags position of id in name.");
64 "Need to specify fromfld=file.fld ");
66 m_f->m_fldToBnd =
false;
77 cout <<
"ProcessMultiShear: Calculating shear stress metrics..." << endl;
80 int nstart, i, j, nfields;
83 string fromfld, basename, endname, nstartStr;
84 stringstream filename;
85 vector<string> infiles(nfld);
86 vector< boost::shared_ptr<Field> > m_fromField(nfld);
89 fromfld =
m_config[
"fromfld"].as<
string>();
90 basename = fromfld.substr(0, fromfld.find_first_of(
"_")+1);
91 filename << fromfld.substr(fromfld.find_first_of(
"_")+1, fromfld.size());
95 filename >> nstartStr;
97 endname = fromfld.substr(fromfld.find(nstartStr)+nstartStr.size(), fromfld.size());
99 for (i=0; i<nfld; ++i)
101 stringstream filename;
102 filename << basename << i+nstart << endname;
103 filename >> infiles[i];
104 cout << infiles[i]<<endl;
107 for ( i = 0; i<nfld; ++i)
109 m_fromField[i] = boost::shared_ptr<Field>(
new Field());
110 m_fromField[i]->m_session =
m_f->m_session;
111 m_fromField[i]->m_graph =
m_f->m_graph;
117 for (i = 0; i < nfld; ++i)
119 if(
m_f->m_exp.size())
123 for (j = 0; j <
m_f->m_exp[0]->GetExpSize(); ++j)
125 ElementGIDs[j] =
m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
127 m_fromField[i]->m_fld->Import(infiles[i],m_fromField[i]->m_fielddef,
128 m_fromField[i]->m_data,
134 m_fromField[i]->m_fld->Import(infiles[i],m_fromField[i]->m_fielddef,
135 m_fromField[i]->m_data,
139 nfields = m_fromField[i]->m_fielddef[0]->m_fields.size();
140 int NumHomogeneousDir = m_fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
142 if (nfields == 5 || nfields == 7)
148 m_fromField[i]->m_graph->SetExpansions(m_fromField[i]->m_fielddef);
151 m_fromField[i]->m_exp.resize(nfields);
152 m_fromField[i]->m_exp[0] = m_fromField[i]->SetUpFirstExpList(NumHomogeneousDir,
true);
154 for(j = 1; j < nfields; ++j)
156 m_fromField[i]->m_exp[j] =
m_f->AppendExpList(NumHomogeneousDir);
159 for (j = 0; j < nfields; ++j)
161 for (
int k = 0; k < m_fromField[i]->m_data.size(); ++k)
163 m_fromField[i]->m_exp[j]->ExtractDataToCoeffs(
164 m_fromField[i]->m_fielddef[k],
165 m_fromField[i]->m_data[k],
166 m_fromField[i]->m_fielddef[k]->m_fields[j],
167 m_fromField[i]->m_exp[j]->UpdateCoeffs());
169 m_fromField[i]->m_exp[j]->BwdTrans(m_fromField[i]->m_exp[j]->GetCoeffs(),
170 m_fromField[i]->m_exp[j]->UpdatePhys());
175 int spacedim =
m_f->m_graph->GetSpaceDimension();
176 if ((m_fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
177 (m_fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 2)
184 if (wssg) { nout = 6; }
186 int npoints = m_fromField[0]->m_exp[0]->GetNpoints();
191 for (i = 0; i < spacedim; ++i)
200 for (i = 0; i < nout; ++i)
208 for (i = 0; i < nfld; ++i)
210 for (j = 0; j < spacedim; ++j)
212 Vmath::Vadd(npoints,m_fromField[i]->m_exp[j]->GetPhys(), 1, normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
216 for (i = 0; i < spacedim; ++i)
218 Vmath::Smul(npoints, 1.0/nfld, normTemporalMeanVec[i], 1, normTemporalMeanVec[i], 1);
219 Vmath::Vvtvp(npoints, normTemporalMeanVec[i], 1, normTemporalMeanVec[i], 1,
220 TemporalMeanMag, 1, TemporalMeanMag, 1);
223 Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
225 for (i = 0; i < spacedim; ++i)
227 Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1, normTemporalMeanVec[i], 1);
233 Vmath::Vmul(npoints,m_fromField[0]->m_exp[nfields-1]->GetPhys(),1,normTemporalMeanVec[1],1,normCrossDir[0],1);
234 Vmath::Vvtvm(npoints,m_fromField[0]->m_exp[nfields-2]->GetPhys(),1,normTemporalMeanVec[2],1,normCrossDir[0],1,
237 Vmath::Vmul(npoints,m_fromField[0]->m_exp[nfields-3]->GetPhys(),1,normTemporalMeanVec[2],1,normCrossDir[1],1);
238 Vmath::Vvtvm(npoints,m_fromField[0]->m_exp[nfields-1]->GetPhys(),1,normTemporalMeanVec[0],1,normCrossDir[1],1,
241 Vmath::Vmul(npoints,m_fromField[0]->m_exp[nfields-2]->GetPhys(),1,normTemporalMeanVec[0],1,normCrossDir[2],1);
242 Vmath::Vvtvm(npoints,m_fromField[0]->m_exp[nfields-3]->GetPhys(),1,normTemporalMeanVec[1],1,normCrossDir[2],1,
248 for (i = 0; i < nfld; ++i)
250 for (j = 0; j < spacedim; ++j)
252 Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1, normTemporalMeanVec[j], 1,
253 DotProduct, 1, DotProduct, 1);
257 Vmath::Vadd(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, outfield[0], 1, outfield[0], 1);
260 Vmath::Vmul(npoints, DotProduct, 1, DotProduct, 1, temp,1);
261 Vmath::Vvtvm(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1,
262 m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1, temp, 1);
264 for (j = 0; j < npoints; ++j)
268 outfield[1][j] = outfield[1][j] + sqrt(temp[j]);
273 Vmath::Vdiv(npoints, DotProduct, 1, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
274 Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
277 for (j = 0; j < npoints; ++j)
279 temp[j] = 1 - temp[j]*temp[j];
282 outfield[4][j] = outfield[4][j] + sqrt(temp[j]);
292 m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct,dtemp[0],dtemp[1],dtemp[2]);
293 for(j=0; j<spacedim; j++)
295 Vmath::Vvtvp(npoints,dtemp[j],1,normTemporalMeanVec[j],1,temp,1,temp,1);
303 for(j = 0; j < spacedim; ++j)
305 Vmath::Vvtvp(npoints,m_fromField[i]->m_exp[j]->GetPhys(),1,normCrossDir[j],1,
306 DotProduct, 1, DotProduct, 1);
309 m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct,dtemp[0],dtemp[1],dtemp[2]);
312 for(j=0; j<spacedim; j++)
314 Vmath::Vvtvp(npoints,dtemp[j],1,normCrossDir[j],1,DotProduct,1,DotProduct,1);
316 Vmath::Vvtvp(npoints,DotProduct,1,DotProduct,1,temp,1, temp,1);
320 Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
328 Vmath::Smul(npoints, 1.0/nfld, outfield[0], 1, outfield[0], 1);
329 Vmath::Smul(npoints, 1.0/nfld, outfield[1], 1, outfield[1], 1);
330 Vmath::Smul(npoints, 1.0/nfld, outfield[3], 1, outfield[3], 1);
331 Vmath::Smul(npoints, 1.0/nfld, outfield[4], 1, outfield[4], 1);
334 for (i = 0; i < npoints; ++i)
336 outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i]/outfield[0][i]);
346 m_f->m_exp.resize(nout);
347 m_f->m_fielddef = m_fromField[0]->m_fielddef;
348 m_f->m_exp[0] =
m_f->SetUpFirstExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir,
true);
350 for(i = 1; i < nout; ++i)
352 m_f->m_exp[i] =
m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
355 m_f->m_fielddef[0]->m_fields.resize(nout);
356 m_f->m_fielddef[0]->m_fields[0] =
"TAWSS";
357 m_f->m_fielddef[0]->m_fields[1] =
"transWSS";
358 m_f->m_fielddef[0]->m_fields[2] =
"OSI";
359 m_f->m_fielddef[0]->m_fields[3] =
"TAAFI";
360 m_f->m_fielddef[0]->m_fields[4] =
"TACFI";
364 m_f->m_fielddef[0]->m_fields[5] =
"|WSSG|";
365 Vmath::Smul(npoints, 1.0/nfld, outfield[5], 1, outfield[5], 1);
369 for(i = 0; i < nout; ++i)
371 m_f->m_exp[i]->FwdTrans(outfield[i],
372 m_f->m_exp[i]->UpdateCoeffs());
373 m_f->m_exp[i]->BwdTrans(
m_f->m_exp[i]->GetCoeffs(),
374 m_f->m_exp[i]->UpdatePhys());
378 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
379 = m_fromField[0]->m_exp[0]->GetFieldDefinitions();
380 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
382 for( i = 0; i < nout; ++i)
384 for ( j = 0; j < FieldDef.size(); ++j)
386 FieldDef[j]->m_fields.push_back(
m_f->m_fielddef[0]->m_fields[i]);
387 m_f->m_exp[i]->AppendFieldData(FieldDef[j], FieldData[j]);
391 m_f->m_fielddef = FieldDef;
392 m_f->m_data = FieldData;
#define ASSERTL0(condition, msg)
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
map< string, ConfigOption > m_config
List of configuration values.
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.
FieldSharedPtr m_f
Field object.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
boost::shared_ptr< Field > FieldSharedPtr
Represents a command-line configuration option.
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
virtual ~ProcessMultiShear()
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.
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.