Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
Nektar::Utilities::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::Utilities::ProcessMultiShear:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::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::Utilities::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
 ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
 
void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
void PrintConfig ()
 Print out all configuration options for a module. More...
 
void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
bool GetRequireEquiSpaced (void)
 
void SetRequireEquiSpaced (bool pVal)
 
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 
 Module (MeshSharedPtr p_m)
 
virtual void Process ()=0
 
void RegisterConfig (std::string key, std::string value)
 
void PrintConfig ()
 
void SetDefaults ()
 
MeshSharedPtr GetMesh ()
 
virtual void ProcessVertices ()
 Extract element vertices. More...
 
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual void ProcessElements ()
 Generate element IDs. More...
 
virtual void ProcessComposites ()
 Generate composites. More...
 
virtual void ClearElementLinks ()
 

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::Utilities::Module
 Module ()
 
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
void PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
MeshSharedPtr m_mesh
 Mesh object. More...
 
std::map< std::string,
ConfigOption
m_config
 List of configuration values. More...
 

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::Utilities::ProcessMultiShear::ProcessMultiShear ( FieldSharedPtr  f)

Definition at line 57 of file ProcessMultiShear.cpp.

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

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

Definition at line 69 of file ProcessMultiShear.cpp.

70 {
71 }

Member Function Documentation

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

Creates an instance of this class.

Definition at line 54 of file ProcessMultiShear.h.

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

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

Implements Nektar::Utilities::Module.

Definition at line 65 of file ProcessMultiShear.h.

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

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 73 of file ProcessMultiShear.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, Nektar::LibUtilities::NullFieldMetaDataMap, Vmath::Smul(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vvtvm(), Vmath::Vvtvp(), and Vmath::Zero().

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

Member Data Documentation

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

Definition at line 57 of file ProcessMultiShear.h.