Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FldAddMultiShear.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 #include <sstream>
4 #include <iostream>
5 
6 #include <MultiRegions/ExpList.h>
7 #include <MultiRegions/ExpList1D.h>
8 #include <MultiRegions/ExpList2D.h>
9 #include <MultiRegions/ExpList3D.h>
10 #include <MultiRegions/ContField3D.h>
11 #include <MultiRegions/ContField2D.h>
12 
13 using namespace Nektar;
14 
15 int main(int argc, char *argv[])
16 {
17  int i,j;
18  int nfiles, nStart, surfID;
19  if(argc != 5)
20  {
21  fprintf(stderr,"Usage: FldAddMultiShear meshfile nfiles FirstInFile BoundaryID\n");
22  exit(1);
23  }
24 
25  nfiles = boost::lexical_cast<int>(argv[2]);
26  surfID = boost::lexical_cast<int>(argv[4]);
27 
28  vector<string> infiles(nfiles);
29  string outfile;
30 
31  // Output file name: multishear.fld
32  stringstream filename2;
33  filename2 << "multishear.fld";
34  filename2 >> outfile;
35 
36  // starting checkpoint file name: name_tn_wss.fld. n can be any number.
37  string basename = argv[3];
38  basename = basename.substr(basename.find_last_of("t")+1, basename.find_last_of(".")-basename.find_last_of("t"));
39  stringstream filename3;
40  filename3 << basename;
41  filename3 >> nStart;
42 
43  for (i = 0; i< nfiles; ++i)
44  {
45  basename = argv[3];
46  string extension = ".fld";
47  basename = basename.substr(0, basename.find_first_of("_"));
48  stringstream filename;
49  filename << basename << "_t" << i+nStart << "_wss.fld";
50  filename >> infiles[i];
51  cout << infiles[i]<<endl;
52  }
53 
54  argv[2] = argv[3];
55  argv[4] = argv[3];
56 
57  LibUtilities::SessionReaderSharedPtr vSession
58  = LibUtilities::SessionReader::CreateInstance(argc, argv);
59 
60 
61  //----------------------------------------------
62  // Read in mesh from input file
63  string meshfile(argv[1]);
64  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);
65  //----------------------------------------------
66 
67 
68  //----------------------------------------------
69  // Import first field file.
70  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
71  vector<vector<NekDouble> > fielddata;
72  LibUtilities::Import(infiles[0],fielddef,fielddata);
73  //----------------------------------------------
74 
75 
76  //----------------------------------------------
77  // Define Expansion
78  //----------------------------------------------
79  int expdim = graphShPt->GetMeshDimension();
80  int nfields = fielddef[0]->m_fields.size();
81  int addfields = 3;
82  int sfields = nfields - expdim;
83  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields+addfields), shear(sfields);
84  MultiRegions::AssemblyMapCGSharedPtr m_locToGlobalMap;
85 
86  switch(expdim)
87  {
88  case 1:
89  {
90  ASSERTL0(false,"Expansion dimension not recognised");
91  }
92  break;
93  case 2:
94  {
95  ASSERTL0(false,"Expansion dimension not recognised");
96  }
97  break;
98  case 3:
99  {
100  i = 0;
101  MultiRegions::ContField3DSharedPtr firstfield =
102  MemoryManager<MultiRegions::ContField3D>
103  ::AllocateSharedPtr(vSession, graphShPt,
104  vSession->GetVariable(i));
105 
106  m_locToGlobalMap = firstfield->GetLocalToGlobalMap();
107 
108  Exp[0] = firstfield;
109  for(i = 1; i < expdim; ++i)
110  {
111  Exp[i] = MemoryManager<MultiRegions::ContField3D>
112  ::AllocateSharedPtr(*firstfield, graphShPt,
113  vSession->GetVariable(i));
114  }
115 
116  for(i = 0; i < sfields; ++i)
117  {
118  Exp[i+expdim] = MemoryManager<MultiRegions::ContField3D>
119  ::AllocateSharedPtr(*firstfield, graphShPt,
120  vSession->GetVariable(i));
121  }
122 
123  for(i = 0; i < addfields; ++i)
124  {
125  Exp[i+nfields] = MemoryManager<MultiRegions::ContField3D>
126  ::AllocateSharedPtr(*firstfield, graphShPt,
127  vSession->GetVariable(i));
128  }
129 
130  }
131  break;
132  default:
133  ASSERTL0(false,"Expansion dimension not recognised");
134  break;
135  }
136 
137 
138  //-------------------------------------
139  // FIRST FILE - INITIALISING EVERYTHING
140  //-------------------------------------
141 
142  // Set up mapping from Boundary condition to element details.
143  StdRegions::StdExpansionSharedPtr elmt;
144  StdRegions::StdExpansion2DSharedPtr bc;
145  Array<OneD, int> BoundarytoElmtID;
146  Array<OneD, int> BoundarytoTraceID;
147  Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> > BndExp(nfields+addfields);
148 
149  // Copy data from field file
150  for(i = 0; i < fielddata.size(); ++i)
151  {
152  Exp[0]->ExtractDataToCoeffs(fielddef[i],
153  fielddata[i],
154  fielddef[i]->m_fields[0],
155  Exp[0]->UpdateCoeffs());
156  }
157  Exp[0]->BwdTrans(Exp[0]->GetCoeffs(),Exp[0]->UpdatePhys());
158 
159  Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
160  BndExp[0] = Exp[0]->GetBndCondExpansions();
161 
162  // Get face 2D expansion from element expansion
163  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][surfID]->GetExp(0));
164 
165  int nfq= bc->GetTotPoints();
166  int nelem = BndExp[0][surfID]->GetExpSize();
167  int nbq = nelem*nfq;
168 
169  int nt = Exp[0]->GetNpoints();
170  Array<OneD, const NekDouble> testing(nt);
171 
172 
173  // Define local arrays for wss components, and outputs (TAwss, osi, trs)
174  int n, cnt, elmtid, offset, boundary, bndOffset;
175  Array<OneD, NekDouble> Sx(nbq), Sy(nbq), Sz(nbq), S(nbq), Sxr(nbq), Syr(nbq), Szr(nbq), sx(nbq), sy(nbq), sz(nbq);
176  Array<OneD, NekDouble> Save(nbq), temp2(nbq), trs(nbq), TAwss(nbq), osi(nbq);
177  Array<OneD, Array<OneD, NekDouble> > temp(sfields), values(nfields);
178 
179  for (i = 0; i < sfields; i++)
180  {
181  temp[i]= Array<OneD, NekDouble>(nfq);
182  }
183 
184  Vmath::Zero (nbq, Sx, 1);
185  Vmath::Zero (nbq, Sy, 1);
186  Vmath::Zero (nbq, Sz, 1);
187 
188  // -----------------------------------------------------
189  // Compute temporal average wall shear stress vector,
190  // and the spatial average of the temporal average.
191  // -----------------------------------------------------
192 
193  for (int fileNo = 0; fileNo < nfiles ; ++fileNo)
194  {
195  // Import field file.
196  fielddef.clear();
197  fielddata.clear();
198  LibUtilities::Import(infiles[fileNo],fielddef,fielddata);
199 
200  // Copy data from field file
201  for(j = 0; j < nfields; ++j)
202  {
203  for(int i = 0; i < fielddata.size(); ++i)
204  {
205  Exp[j]->ExtractDataToCoeffs(fielddef[i],
206  fielddata[i],
207  fielddef[i]->m_fields[j],
208  Exp[j]->UpdateCoeffs());
209  }
210  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
211  }
212 
213  Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
214  BndExp[0] = Exp[0]->GetBndCondExpansions();
215 
216  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
217  {
218  if(n == surfID)
219  {
220  for(i = 0; i < nelem; ++i, cnt++)
221  {
222  // find element and face of this expansion.
223  elmtid = BoundarytoElmtID[cnt];
224  elmt = Exp[0]->GetExp(elmtid);
225  offset = Exp[0]->GetPhys_Offset(elmtid);
226  bndOffset = nfq*i;
227 
228  if(expdim == 2)
229  {
230  // Not implemented in 2D.
231  }
232  else
233  {
234  // Get face 2D expansion from element expansion
235  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
236 
237  //identify boundary of element looking at.
238  boundary = BoundarytoTraceID[cnt];
239 
240  // Get face stress values.
241  for (int t = 0; t< expdim; t++)
242  {
243  elmt->GetFacePhysVals(boundary,bc, Exp[t+expdim]->GetPhys() + offset, temp[t]);
244  }
245 
246  for (int m = 0; m < nfq; m++)
247  {
248  Sxr[bndOffset + m] = temp[0][m];
249  Syr[bndOffset + m] = temp[1][m];
250  Szr[bndOffset + m] = temp[2][m];
251  }
252  }
253  }
254 
255  // Sx = Sx + Sxr;
256  Vmath::Vadd (nbq, Sxr, 1, Sx, 1, Sx, 1);
257  Vmath::Vadd (nbq, Syr, 1, Sy, 1, Sy, 1);
258  Vmath::Vadd (nbq, Szr, 1, Sz, 1, Sz, 1);
259 
260  Vmath::Zero (nbq, Sxr, 1);
261  Vmath::Zero (nbq, Syr, 1);
262  Vmath::Zero (nbq, Szr, 1);
263  }
264  else
265  {
266  cnt += BndExp[0][n]->GetExpSize();
267  }
268  }
269  }
270 
271  // Temporal averages of each wss component: Sx = Sx / nfiles. I.e. mean temporal wss
272  Vmath::Smul(nbq, 1.0/nfiles , Sx, 1, Sx, 1);
273  Vmath::Smul(nbq, 1.0/nfiles , Sy, 1, Sy, 1);
274  Vmath::Smul(nbq, 1.0/nfiles , Sz, 1, Sz, 1);
275 
276  // Store temporal mean vectors sx, sy, sz
277  Vmath::Vcopy(nbq, Sx, 1, sx, 1);
278  Vmath::Vcopy(nbq, Sy, 1, sy, 1);
279  Vmath::Vcopy(nbq, Sz, 1, sz, 1);
280 
281  // Spatial average of the temporal averaged wss vector: Save = sqrt(sx^2 + Sy^2 + Sz^2);
282  // i.e magnitude of mean temporal wss.
283  Vmath::Vvtvvtp(nbq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, Save, 1);
284  Vmath::Vvtvp(nbq, Sz, 1, Sz, 1, Save, 1, Save, 1);
285  Vmath::Vsqrt(nbq, Save, 1, Save, 1);
286 
287  //temporal averaged vector / spatial average (sort of normalisation): Sxr = Sx/Save;
288  //unit vector of mean temporal wss (t).
289  Vmath::Vdiv(nbq, Sx, 1, Save, 1, Sxr, 1);
290  Vmath::Vdiv(nbq, Sy, 1, Save, 1, Syr, 1);
291  Vmath::Vdiv(nbq, Sz, 1, Save, 1, Szr, 1);
292 
293  Vmath::Zero (nbq, Sx, 1);
294  Vmath::Zero (nbq, Sy, 1);
295  Vmath::Zero (nbq, Sz, 1);
296  Vmath::Zero (nbq, trs, 1);
297  Vmath::Zero (nbq, TAwss, 1);
298 
299  // -----------------------------------------------------
300  // Loop through files again, and compute transverse wss
301  // -----------------------------------------------------
302 
303  for (int fileNo = 0; fileNo < nfiles ; ++fileNo)
304  {
305  Vmath::Zero (nbq, Sx, 1);
306  Vmath::Zero (nbq, Sy, 1);
307  Vmath::Zero (nbq, Sz, 1);
308 
309  // Import field file.
310  fielddef.clear();
311  fielddata.clear();
312  LibUtilities::Import(infiles[fileNo],fielddef,fielddata);
313 
314  // Copy data from field file
315  for(j = 0; j < nfields; ++j)
316  {
317  for(int i = 0; i < fielddata.size(); ++i)
318  {
319  Exp[j]->ExtractDataToCoeffs(fielddef[i],
320  fielddata[i],
321  fielddef[i]->m_fields[j],
322  Exp[j]->UpdateCoeffs());
323  }
324  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
325  }
326 
327  Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
328  BndExp[0] = Exp[0]->GetBndCondExpansions();
329 
330  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
331  {
332  if(n == surfID)
333  {
334  for(i = 0; i < nelem; ++i, cnt++)
335  {
336  // find element and face of this expansion.
337  elmtid = BoundarytoElmtID[cnt];
338  elmt = Exp[0]->GetExp(elmtid);
339  offset = Exp[0]->GetPhys_Offset(elmtid);
340  bndOffset = nfq*i;
341 
342  if(expdim == 2)
343  {
344  }
345  else
346  {
347  // Get face 2D expansion from element expansion
348  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
349 
350  //identify boundary of element looking at.
351  boundary = BoundarytoTraceID[cnt];
352 
353  // Get face stress values.
354  for (int t = 0; t< sfields; t++)
355  {
356  elmt->GetFacePhysVals(boundary,bc,Exp[t+expdim]->GetPhys() + offset, temp[t]);
357  }
358 
359  for (int m = 0; m < nfq; m++)
360  {
361  Sx[bndOffset + m] = temp[0][m];
362  Sy[bndOffset + m] = temp[1][m];
363  Sz[bndOffset + m] = temp[2][m];
364  //S [bndOffset + m] = temp[3][m];
365  }
366  }
367  }
368 
369  Vmath::Vvtvvtp(nbq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, S, 1);
370  Vmath::Vvtvp(nbq, Sz, 1, Sz, 1, S, 1, S, 1);
371  Vmath::Vsqrt(nbq, S, 1, S, 1);
372 
373  // transverse wall shear stress
374  // trs = trs + sqrt(s^2 - (Sx*Sxr + Sy*Syr + Sz*Szr)^2)
375  Vmath::Vvtvvtp(nbq, Sx, 1, Sxr, 1, Sy, 1, Syr, 1, temp2, 1);
376  Vmath::Vvtvp(nbq, Sz, 1, Szr, 1, temp2, 1, temp2, 1);
377  Vmath::Vmul(nbq, temp2, 1, temp2, 1, temp2, 1);
378  Vmath::Vvtvm(nbq, S, 1, S, 1, temp2, 1, temp2, 1);
379  for (int m = 0; m < nbq; m++)
380  {
381  if (temp2[m]> 0.0)
382  {
383  trs[m] = trs[m] + sqrt(temp2[m]);
384  }
385  else
386  {
387  trs[m] = 0.0;
388  }
389  }
390 
391  Vmath::Zero (nbq, temp2, 1);
392 
393  //Time averaged wss TAwss = TAwss + S;
394  Vmath::Vadd (nbq, S, 1, TAwss, 1, TAwss, 1);
395  }
396  else
397  {
398  cnt += BndExp[0][n]->GetExpSize();
399  }
400  }
401  }
402 
403 
404  // Final step in computing trs and TAwss: trs = trs/nfiles;
405  Vmath::Smul (nbq, (1.0/nfiles), trs, 1, trs, 1);
406  Vmath::Smul (nbq, (1.0/nfiles), TAwss, 1, TAwss, 1);
407 
408 
409  //Compute osi = 0.5*(1- Save/TAwss)
410  Vmath::Vdiv(nbq, Save, 1, TAwss, 1, osi, 1);
411  Vmath::Smul(nbq, -0.5, osi, 1, osi, 1);
412  Vmath::Sadd(nbq, 0.5, osi, 1, osi, 1);
413 
414  fielddef.clear();
415  fielddata.clear();
416  LibUtilities::Import(infiles[nfiles-1],fielddef,fielddata);
417 
418  // Copy u,v,w, data from last field file for the fields
419  for(j = 0; j < nfields; ++j)
420  {
421  for(int i = 0; i < fielddata.size(); ++i)
422  {
423  Exp[j]->ExtractDataToCoeffs(fielddef[i],
424  fielddata[i],
425  fielddef[i]->m_fields[j],
426  Exp[j]->UpdateCoeffs());
427  }
428  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
429  }
430 
431  for(j = 0; j < addfields; ++j)
432  {
433  for(int i = 0; i < fielddata.size(); ++i)
434  {
435  Exp[j+nfields]->ExtractDataToCoeffs(fielddef[i],
436  fielddata[i],
437  fielddef[i]->m_fields[j],
438  Exp[j+nfields]->UpdateCoeffs());
439  }
440  Exp[j+nfields]->BwdTrans(Exp[j+nfields]->GetCoeffs(),Exp[j+nfields]->UpdatePhys());
441  }
442 
443  Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
444 
445  //get boundary expansions for each field
446  for(j = 0; j < nfields+addfields; ++j)
447  {
448  BndExp[j] = Exp[j]->GetBndCondExpansions();
449  }
450 
451  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
452  {
453  if(n == surfID)
454  {
455  for(i = 0; i < nelem; ++i, cnt++)
456  {
457  // find element and face of this expansion.
458  elmtid = BoundarytoElmtID[cnt];
459  offset = Exp[0]->GetPhys_Offset(elmtid);
460  bndOffset = nfq*i;
461 
462  if(expdim == 2)
463  {
464  }
465  else
466  {
467  // Get face 2D expansion from element expansion
468  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
469 
470  //identify boundary of element looking at.
471  boundary = BoundarytoTraceID[cnt];
472 
473  //Update outfield coefficients in the elemental boundary expansion
474  for (j = 0; j < nfq; j++)
475  {
476  temp[0][j] = trs[bndOffset + j];
477  temp[1][j] = TAwss[bndOffset + j];
478  temp[2][j] = osi[bndOffset + j];
479  }
480 
481  for (j = 0; j < addfields; j++)
482  {
483  values[j] = BndExp[j+nfields][n]->UpdateCoeffs() + BndExp[j+nfields][n]->GetCoeff_Offset(i);
484  bc->FwdTrans(temp[j], values[j]);
485  }
486  /*
487  for (j = 0; j < nfq; j++)
488  {
489  temp[0][j] = sx[bndOffset + j];
490  temp[1][j] = sy[bndOffset + j];
491  temp[2][j] = sz[bndOffset + j];
492  }
493 
494  for (j = 0; j < expdim; j++)
495  {
496  values[j] = BndExp[j+addfields][n]->UpdateCoeffs() + BndExp[j+addfields][n]->GetCoeff_Offset(i);
497  bc->FwdTrans(temp[j], values[j]);
498  }
499  */
500  }
501  }
502  }
503  else
504  {
505  cnt += BndExp[0][n]->GetExpSize();
506  }
507  }
508 
509  for(j = 0; j < nfields + addfields; ++j)
510  {
511  int ncoeffs = Exp[j]->GetNcoeffs();
512  Array<OneD, NekDouble> output(ncoeffs);
513 
514  output=Exp[j]->UpdateCoeffs();
515 
516  int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
517  Array<OneD, NekDouble> outarray(nGlobal,0.0);
518 
519  int bndcnt=0;
520 
521  const Array<OneD,const int>& map = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsMap();
522  NekDouble sign;
523 
524  for(i = 0; i < BndExp[j].num_elements(); ++i)
525  {
526  if(i==surfID)
527  {
528  const Array<OneD,const NekDouble>& coeffs = BndExp[j][i]->GetCoeffs();
529  for(int k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
530  {
531  sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
532  outarray[map[bndcnt++]] = sign * coeffs[k];
533  }
534  }
535  else
536  {
537  bndcnt += BndExp[j][i]->GetNcoeffs();
538  }
539  }
540  m_locToGlobalMap->GlobalToLocal(outarray,output);
541 
542  }
543 
544  //-----------------------------------------------
545  // Write solution to file
546  // ----------------------------------------------
547  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
548  = Exp[0]->GetFieldDefinitions();
549  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
550 
551  vector<string > outname;
552 
553  outname.push_back("TransWSS");
554  outname.push_back("TAWSS");
555  outname.push_back("OSI");
556 
557  for(j = 0; j < nfields+addfields; ++j)
558  {
559  for(i = 0; i < FieldDef.size(); ++i)
560  {
561  if (j >= nfields)
562  {
563  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
564  }
565  else
566  {
567  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
568  }
569  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
570  }
571  }
572 
573  LibUtilities::Write(outfile, FieldDef, FieldData);
574 
575  return 0;
576 }
577 
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2