Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
Nektar::FieldUtils::ProcessMultiShear Class Reference

This processing module calculates the shear stress metrics and writes it to a surface output file. More...

#include <ProcessMultiShear.h>

Inheritance diagram for Nektar::FieldUtils::ProcessMultiShear:
Inheritance graph
[legend]
Collaboration diagram for Nektar::FieldUtils::ProcessMultiShear:
Collaboration graph
[legend]

Public Member Functions

 ProcessMultiShear (FieldSharedPtr f)
 
virtual ~ProcessMultiShear ()
 
virtual void Process (po::variables_map &vm)
 Write mesh to output file. More...
 
virtual std::string GetModuleName ()
 
- Public Member Functions inherited from Nektar::FieldUtils::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
- Public Member Functions inherited from Nektar::FieldUtils::Module
FIELD_UTILS_EXPORT Module (FieldSharedPtr p_f)
 
FIELD_UTILS_EXPORT void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
FIELD_UTILS_EXPORT void PrintConfig ()
 Print out all configuration options for a module. More...
 
FIELD_UTILS_EXPORT void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
FIELD_UTILS_EXPORT bool GetRequireEquiSpaced (void)
 
FIELD_UTILS_EXPORT void SetRequireEquiSpaced (bool pVal)
 
FIELD_UTILS_EXPORT void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 

Static Public Member Functions

static boost::shared_ptr< Modulecreate (FieldSharedPtr f)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::FieldUtils::Module
 Module ()
 
- Protected Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 

Detailed Description

This processing module calculates the shear stress metrics and writes it to a surface output file.

Definition at line 50 of file ProcessMultiShear.h.

Constructor & Destructor Documentation

Nektar::FieldUtils::ProcessMultiShear::ProcessMultiShear ( FieldSharedPtr  f)

Definition at line 58 of file ProcessMultiShear.cpp.

References ASSERTL0, Nektar::FieldUtils::Module::m_config, and Nektar::FieldUtils::Module::m_f.

58  : ProcessModule(f)
59 {
60  m_config["N"] = ConfigOption(false, "1", "Number of chk or fld files");
61  m_config["fromfld"] = ConfigOption(
62  false, "NotSet",
63  "First fld file. First underscore flags position of id in name.");
64 
65  ASSERTL0(m_config["fromfld"].as<string>().compare("NotSet") != 0,
66  "Need to specify fromfld=file.fld ");
67 
68  m_f->m_fldToBnd = false;
69 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
FieldSharedPtr m_f
Field object.
Nektar::FieldUtils::ProcessMultiShear::~ProcessMultiShear ( )
virtual

Definition at line 71 of file ProcessMultiShear.cpp.

72 {
73 }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::FieldUtils::ProcessMultiShear::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 54 of file ProcessMultiShear.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

55  {
57  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual std::string Nektar::FieldUtils::ProcessMultiShear::GetModuleName ( )
inlinevirtual

Implements Nektar::FieldUtils::Module.

Definition at line 66 of file ProcessMultiShear.h.

67  {
68  return "ProcessMultiShear";
69  }
void Nektar::FieldUtils::ProcessMultiShear::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Implements Nektar::FieldUtils::Module.

Definition at line 75 of file ProcessMultiShear.cpp.

References Nektar::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, Nektar::LibUtilities::NullFieldMetaDataMap, Vmath::Smul(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vvtvm(), Vmath::Vvtvp(), and Vmath::Zero().

76 {
77  if (m_f->m_verbose)
78  {
79  if (m_f->m_comm->TreatAsRankZero())
80  {
81  cout << "ProcessMultiShear: Calculating shear stress metrics..."
82  << endl;
83  }
84  }
85 
86  int nstart, i, j, nfields=0;
87  bool wssg = false;
88  NekDouble nfld = m_config["N"].as<NekDouble>();
89  string fromfld, basename, endname, nstartStr;
90  stringstream filename;
91  vector<string> infiles(nfld);
92  vector<boost::shared_ptr<Field> > m_fromField(nfld);
93 
94  // Set up list of input fld files.
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());
98  filename >> nstart;
99  filename.str("");
100  filename << nstart;
101  filename >> nstartStr;
102  filename.str("");
103  endname = fromfld.substr(fromfld.find(nstartStr) + nstartStr.size(),
104  fromfld.size());
105 
106  for (i = 0; i < nfld; ++i)
107  {
108  stringstream filename;
109  filename << basename << i + nstart << endname;
110  filename >> infiles[i];
111  cout << infiles[i] << endl;
112  }
113 
114  for (i = 0; i < nfld; ++i)
115  {
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;
119  }
120 
121  // Import all fld files.
122  for (i = 0; i < nfld; ++i)
123  {
124  if (m_f->m_exp.size())
125  {
126  // Set up ElementGIDs in case of parallel processing
127  Array<OneD, int> ElementGIDs(m_f->m_exp[0]->GetExpSize());
128  for (j = 0; j < m_f->m_exp[0]->GetExpSize(); ++j)
129  {
130  ElementGIDs[j] =
131  m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
132  }
133  m_f->FieldIOForFile(infiles[i])->Import(
134  infiles[i], m_fromField[i]->m_fielddef, m_fromField[i]->m_data,
136  }
137  else
138  {
139  m_f->FieldIOForFile(infiles[i])->Import(
140  infiles[i], m_fromField[i]->m_fielddef, m_fromField[i]->m_data,
142  }
143 
144  nfields = m_fromField[i]->m_fielddef[0]->m_fields.size();
145  int NumHomogeneousDir =
146  m_fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
147 
148  if (nfields == 5 || nfields == 7)
149  {
150  wssg = true;
151  }
152 
153  // Set up Expansion information to use mode order from field
154  m_fromField[i]->m_graph->SetExpansions(m_fromField[i]->m_fielddef);
155 
156  // Set up expansions, and extract data.
157  m_fromField[i]->m_exp.resize(nfields);
158  m_fromField[i]->m_exp[0] =
159  m_fromField[i]->SetUpFirstExpList(NumHomogeneousDir, true);
160 
161  for (j = 1; j < nfields; ++j)
162  {
163  m_fromField[i]->m_exp[j] = m_f->AppendExpList(NumHomogeneousDir);
164  }
165 
166  for (j = 0; j < nfields; ++j)
167  {
168  for (int k = 0; k < m_fromField[i]->m_data.size(); ++k)
169  {
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());
174  }
175  m_fromField[i]->m_exp[j]->BwdTrans(
176  m_fromField[i]->m_exp[j]->GetCoeffs(),
177  m_fromField[i]->m_exp[j]->UpdatePhys());
178  }
179  }
180 
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)
184  {
185  spacedim = 3;
186  }
187 
188  int nout = 5; // TAWSS, OSI, transWSS, AFI, CFI (optional WSSG)
189  if (wssg)
190  {
191  nout = 6;
192  }
193 
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);
199 
200  for (i = 0; i < spacedim; ++i)
201  {
202  normTemporalMeanVec[i] = Array<OneD, NekDouble>(npoints);
203  normCrossDir[i] = Array<OneD, NekDouble>(npoints);
204  dtemp[i] = Array<OneD, NekDouble>(npoints);
205  Vmath::Zero(npoints, normTemporalMeanVec[i], 1);
206  Vmath::Zero(npoints, normCrossDir[i], 1);
207  }
208 
209  for (i = 0; i < nout; ++i)
210  {
211  outfield[i] = Array<OneD, NekDouble>(npoints, 0.0);
212  }
213 
214  // -----------------------------------------------------
215  // Compute temporal average wall shear stress vector,
216  // it's spatial average, and normalise it.
217  for (i = 0; i < nfld; ++i)
218  {
219  for (j = 0; j < spacedim; ++j)
220  {
221  Vmath::Vadd(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1,
222  normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
223  }
224  }
225 
226  for (i = 0; i < spacedim; ++i)
227  {
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);
232  }
233 
234  Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
235 
236  for (i = 0; i < spacedim; ++i)
237  {
238  Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1,
239  normTemporalMeanVec[i], 1);
240  }
241  //---------------------------------------------------
242 
243  if (wssg) // cross product with normals to obtain direction parallel to
244  // temporal mean vector.
245  {
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,
250  normCrossDir[0], 1);
251 
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,
256  normCrossDir[1], 1);
257 
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,
262  normCrossDir[2], 1);
263  }
264 
265  // Compute tawss, trs, osi, taafi, tacfi, WSSG.
266  for (i = 0; i < nfld; ++i)
267  {
268  for (j = 0; j < spacedim; ++j)
269  {
270  Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1,
271  normTemporalMeanVec[j], 1, DotProduct, 1, DotProduct,
272  1);
273  }
274 
275  // TAWSS
276  Vmath::Vadd(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1,
277  outfield[0], 1, outfield[0], 1);
278 
279  // transWSS
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,
283  temp, 1);
284 
285  for (j = 0; j < npoints; ++j)
286  {
287  if (temp[j] > 0.0)
288  {
289  outfield[1][j] = outfield[1][j] + sqrt(temp[j]);
290  }
291  }
292 
293  // TAAFI
294  Vmath::Vdiv(npoints, DotProduct, 1,
295  m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
296  Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
297 
298  // TACFI
299  for (j = 0; j < npoints; ++j)
300  {
301  temp[j] = 1 - temp[j] * temp[j];
302  if (temp[j] > 0.0)
303  {
304  outfield[4][j] = outfield[4][j] + sqrt(temp[j]);
305  }
306  }
307 
308  // WSSG
309  if (wssg)
310  {
311  // parallel component:
312  Vmath::Zero(npoints, temp, 1);
313 
314  m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
315  dtemp[2]);
316  for (j = 0; j < spacedim; j++)
317  {
318  Vmath::Vvtvp(npoints, dtemp[j], 1, normTemporalMeanVec[j], 1,
319  temp, 1, temp, 1);
320  }
321  Vmath::Vmul(npoints, temp, 1, temp, 1, temp, 1);
322 
323  // perpendicular component:
324  Vmath::Zero(npoints, DotProduct, 1);
325 
326  for (j = 0; j < spacedim; ++j)
327  {
328  Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1,
329  normCrossDir[j], 1, DotProduct, 1, DotProduct, 1);
330  }
331  m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
332  dtemp[2]);
333  Vmath::Zero(npoints, DotProduct, 1);
334 
335  for (j = 0; j < spacedim; j++)
336  {
337  Vmath::Vvtvp(npoints, dtemp[j], 1, normCrossDir[j], 1,
338  DotProduct, 1, DotProduct, 1);
339  }
340  Vmath::Vvtvp(npoints, DotProduct, 1, DotProduct, 1, temp, 1, temp,
341  1);
342 
343  // Final wssg computation
344  Vmath::Vsqrt(npoints, temp, 1, temp, 1);
345  Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
346  }
347 
348  Vmath::Zero(npoints, DotProduct, 1);
349  }
350 
351  // Divide by nfld
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);
356 
357  // OSI
358  for (i = 0; i < npoints; ++i)
359  {
360  outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i] / outfield[0][i]);
361  }
362 
363  /* TAWSS = sum(wss)/nfld
364  * transWSS = sum( sqrt( wss^2 - (wss . normTempMean)^2) )/nfld.
365  * OSI = 0.5*(1-TemporalMeanMag/TAWSS)
366  * TAAFI = sum(cos)/nfld
367  * TACFI = sum(sin)/nfld = sum( sqrt(1-cos^2) )/nfld.
368  */
369 
370  m_f->m_exp.resize(nout);
371  m_f->m_fielddef = m_fromField[0]->m_fielddef;
372  m_f->m_exp[0] =
373  m_f->SetUpFirstExpList(m_f->m_fielddef[0]->m_numHomogeneousDir, true);
374 
375  for (i = 1; i < nout; ++i)
376  {
377  m_f->m_exp[i] =
378  m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
379  }
380 
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";
387 
388  if (wssg)
389  {
390  m_f->m_fielddef[0]->m_fields[5] = "|WSSG|";
391  Vmath::Smul(npoints, 1.0 / nfld, outfield[5], 1, outfield[5], 1);
392  }
393 
394  for (i = 0; i < nout; ++i)
395  {
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());
399  }
400 
401  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
402  m_fromField[0]->m_exp[0]->GetFieldDefinitions();
403  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
404 
405  for (i = 0; i < nout; ++i)
406  {
407  for (j = 0; j < FieldDef.size(); ++j)
408  {
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]);
411  }
412  }
413 
414  m_f->m_fielddef = FieldDef;
415  m_f->m_data = FieldData;
416 }
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)
Definition: Vmath.cpp:408
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
Definition: Vmath.cpp:442
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.
Definition: Vmath.cpp:241
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
double NekDouble
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
Definition: Vmath.cpp:465
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:55
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.
Definition: Vmath.cpp:299
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.
Definition: Vmath.cpp:183
FieldSharedPtr m_f
Field object.

Member Data Documentation

ModuleKey Nektar::FieldUtils::ProcessMultiShear::className
static
Initial value:

Definition at line 58 of file ProcessMultiShear.h.