Nektar++
ProcessMultiShear.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessMultiShear.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Computes tawss, osi, transwss, afi, cfi fields.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 #include <sstream>
39 using namespace std;
40 
41 #include "ProcessMultiShear.h"
42 
45 #include <MultiRegions/ExpList.h>
46 
47 namespace Nektar
48 {
49 namespace Utilities
50 {
51 
52 ModuleKey ProcessMultiShear::className =
54  ModuleKey(eProcessModule, "shear"),
55  ProcessMultiShear::create, "Computes shear stress metrics.");
56 
57 ProcessMultiShear::ProcessMultiShear(FieldSharedPtr f) : 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 }
68 
70 {
71 }
72 
73 void ProcessMultiShear::Process(po::variables_map &vm)
74 {
75  if (m_f->m_verbose)
76  {
77  cout << "ProcessMultiShear: Calculating shear stress metrics..." << endl;
78  }
79 
80  int nstart, i, j, nfields;
81  bool wssg = false;
82  NekDouble nfld = m_config["N"].as<NekDouble>();
83  string fromfld, basename, endname, nstartStr;
84  stringstream filename;
85  vector<string> infiles(nfld);
86  vector< boost::shared_ptr<Field> > m_fromField(nfld);
87 
88  // Set up list of input fld files.
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());
92  filename >> nstart;
93  filename.str("");
94  filename << nstart;
95  filename >> nstartStr;
96  filename.str("");
97  endname = fromfld.substr(fromfld.find(nstartStr)+nstartStr.size(), fromfld.size());
98 
99  for (i=0; i<nfld; ++i)
100  {
101  stringstream filename;
102  filename << basename << i+nstart << endname;
103  filename >> infiles[i];
104  cout << infiles[i]<<endl;
105  }
106 
107  for ( i = 0; i<nfld; ++i)
108  {
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;
112  m_fromField[i]->m_fld = MemoryManager<LibUtilities::FieldIO>
113  ::AllocateSharedPtr(m_fromField[0]->m_session->GetComm());
114  }
115 
116  //Import all fld files.
117  for (i = 0; i < nfld; ++i)
118  {
119  if(m_f->m_exp.size())
120  {
121  // Set up ElementGIDs in case of parallel processing
122  Array<OneD,int> ElementGIDs(m_f->m_exp[0]->GetExpSize());
123  for (j = 0; j < m_f->m_exp[0]->GetExpSize(); ++j)
124  {
125  ElementGIDs[j] = m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
126  }
127  m_fromField[i]->m_fld->Import(infiles[i],m_fromField[i]->m_fielddef,
128  m_fromField[i]->m_data,
130  ElementGIDs);
131  }
132  else
133  {
134  m_fromField[i]->m_fld->Import(infiles[i],m_fromField[i]->m_fielddef,
135  m_fromField[i]->m_data,
137  }
138 
139  nfields = m_fromField[i]->m_fielddef[0]->m_fields.size();
140  int NumHomogeneousDir = m_fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
141 
142  if (nfields == 5 || nfields == 7)
143  {
144  wssg = true;
145  }
146 
147  // Set up Expansion information to use mode order from field
148  m_fromField[i]->m_graph->SetExpansions(m_fromField[i]->m_fielddef);
149 
150  //Set up expansions, and extract data.
151  m_fromField[i]->m_exp.resize(nfields);
152  m_fromField[i]->m_exp[0] = m_fromField[i]->SetUpFirstExpList(NumHomogeneousDir,true);
153 
154  for(j = 1; j < nfields; ++j)
155  {
156  m_fromField[i]->m_exp[j] = m_f->AppendExpList(NumHomogeneousDir);
157  }
158 
159  for (j = 0; j < nfields; ++j)
160  {
161  for (int k = 0; k < m_fromField[i]->m_data.size(); ++k)
162  {
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());
168  }
169  m_fromField[i]->m_exp[j]->BwdTrans(m_fromField[i]->m_exp[j]->GetCoeffs(),
170  m_fromField[i]->m_exp[j]->UpdatePhys());
171  }
172  }
173 
174 
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)
178  {
179  spacedim = 3;
180  }
181 
182 
183  int nout = 5; // TAWSS, OSI, transWSS, AFI, CFI (optional WSSG)
184  if (wssg) { nout = 6; }
185 
186  int npoints = m_fromField[0]->m_exp[0]->GetNpoints();
187  Array<OneD, Array<OneD, NekDouble> > normTemporalMeanVec(spacedim),normCrossDir(spacedim), outfield(nout), dtemp(spacedim);
188  Array<OneD, NekDouble> TemporalMeanMag(npoints,0.0), DotProduct(npoints,0.0), temp(npoints,0.0);
189 
190 
191  for (i = 0; i < spacedim; ++i)
192  {
193  normTemporalMeanVec[i] = Array<OneD, NekDouble>(npoints);
194  normCrossDir[i] = Array<OneD, NekDouble>(npoints);
195  dtemp[i] = Array<OneD, NekDouble>(npoints);
196  Vmath::Zero(npoints, normTemporalMeanVec[i],1);
197  Vmath::Zero(npoints, normCrossDir[i],1);
198  }
199 
200  for (i = 0; i < nout; ++i)
201  {
202  outfield[i] = Array<OneD, NekDouble>(npoints, 0.0);
203  }
204 
205  // -----------------------------------------------------
206  // Compute temporal average wall shear stress vector,
207  // it's spatial average, and normalise it.
208  for (i = 0; i < nfld; ++i)
209  {
210  for (j = 0; j < spacedim; ++j)
211  {
212  Vmath::Vadd(npoints,m_fromField[i]->m_exp[j]->GetPhys(), 1, normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
213  }
214  }
215 
216  for (i = 0; i < spacedim; ++i)
217  {
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);
221  }
222 
223  Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
224 
225  for (i = 0; i < spacedim; ++i)
226  {
227  Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1, normTemporalMeanVec[i], 1);
228  }
229  //---------------------------------------------------
230 
231  if (wssg) //cross product with normals to obtain direction parallel to temporal mean vector.
232  {
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,
235  normCrossDir[0],1);
236 
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,
239  normCrossDir[1],1);
240 
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,
243  normCrossDir[2],1);
244  }
245 
246 
247  // Compute tawss, trs, osi, taafi, tacfi, WSSG.
248  for (i = 0; i < nfld; ++i)
249  {
250  for (j = 0; j < spacedim; ++j)
251  {
252  Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1, normTemporalMeanVec[j], 1,
253  DotProduct, 1, DotProduct, 1);
254  }
255 
256  //TAWSS
257  Vmath::Vadd(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, outfield[0], 1, outfield[0], 1);
258 
259  //transWSS
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);
263 
264  for (j = 0; j < npoints; ++j)
265  {
266  if(temp[j] > 0.0)
267  {
268  outfield[1][j] = outfield[1][j] + sqrt(temp[j]);
269  }
270  }
271 
272  //TAAFI
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);
275 
276  //TACFI
277  for (j = 0; j < npoints; ++j)
278  {
279  temp[j] = 1 - temp[j]*temp[j];
280  if(temp[j] > 0.0)
281  {
282  outfield[4][j] = outfield[4][j] + sqrt(temp[j]);
283  }
284  }
285 
286  // WSSG
287  if (wssg)
288  {
289  //parallel component:
290  Vmath::Zero(npoints, temp,1);
291 
292  m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct,dtemp[0],dtemp[1],dtemp[2]);
293  for(j=0; j<spacedim; j++)
294  {
295  Vmath::Vvtvp(npoints,dtemp[j],1,normTemporalMeanVec[j],1,temp,1,temp,1);
296  }
297  Vmath::Vmul(npoints,temp,1,temp,1,temp,1);
298 
299 
300  //perpendicular component:
301  Vmath::Zero(npoints, DotProduct,1);
302 
303  for(j = 0; j < spacedim; ++j)
304  {
305  Vmath::Vvtvp(npoints,m_fromField[i]->m_exp[j]->GetPhys(),1,normCrossDir[j],1,
306  DotProduct, 1, DotProduct, 1);
307 
308  }
309  m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct,dtemp[0],dtemp[1],dtemp[2]);
310  Vmath::Zero(npoints, DotProduct,1);
311 
312  for(j=0; j<spacedim; j++)
313  {
314  Vmath::Vvtvp(npoints,dtemp[j],1,normCrossDir[j],1,DotProduct,1,DotProduct,1);
315  }
316  Vmath::Vvtvp(npoints,DotProduct,1,DotProduct,1,temp,1, temp,1);
317 
318  //Final wssg computation
319  Vmath::Vsqrt(npoints,temp,1,temp,1);
320  Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
321  }
322 
323 
324  Vmath::Zero(npoints, DotProduct,1);
325  }
326 
327  //Divide by nfld
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);
332 
333  //OSI
334  for (i = 0; i < npoints; ++i)
335  {
336  outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i]/outfield[0][i]);
337  }
338 
339  /* TAWSS = sum(wss)/nfld
340  * transWSS = sum( sqrt( wss^2 - (wss . normTempMean)^2) )/nfld.
341  * OSI = 0.5*(1-TemporalMeanMag/TAWSS)
342  * TAAFI = sum(cos)/nfld
343  * TACFI = sum(sin)/nfld = sum( sqrt(1-cos^2) )/nfld.
344  */
345 
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);
349 
350  for(i = 1; i < nout; ++i)
351  {
352  m_f->m_exp[i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
353  }
354 
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";
361 
362  if (wssg)
363  {
364  m_f->m_fielddef[0]->m_fields[5] = "|WSSG|";
365  Vmath::Smul(npoints, 1.0/nfld, outfield[5], 1, outfield[5], 1);
366  }
367 
368 
369  for(i = 0; i < nout; ++i)
370  {
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());
375  }
376 
377 
378  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
379  = m_fromField[0]->m_exp[0]->GetFieldDefinitions();
380  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
381 
382  for( i = 0; i < nout; ++i)
383  {
384  for ( j = 0; j < FieldDef.size(); ++j)
385  {
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]);
388  }
389  }
390 
391  m_f->m_fielddef = FieldDef;
392  m_f->m_data = FieldData;
393 }
394 
395 }
396 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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)
Definition: Vmath.cpp:394
virtual void Process()=0
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
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
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
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:677
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
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:66
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
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215