Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
FldAddWSS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <sstream>
#include <iostream>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/ContField2D.h>
Include dependency graph for FldAddWSS.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 16 of file FldAddWSS.cpp.

References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::LibUtilities::Import(), sign, Vmath::Smul(), Vmath::Svtsvtp(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vsqrt(), Vmath::Vvtvp(), Vmath::Vvtvvtp(), and Nektar::LibUtilities::Write().

17 {
18  int i,j,k;
19  int surfID, nfiles, nStart;
20 
21  if(argc != 5)
22  {
23  fprintf(stderr,"Usage: FldAddWSS meshfile nfiles infld BoundaryID\n");
24  exit(1);
25  }
26 
27  surfID = boost::lexical_cast<int>(argv[argc - 1]);
28  nfiles = boost::lexical_cast<int>(argv[2]);
29 
30  vector<string> infiles(nfiles), outfiles(nfiles);
31 
32  /* format of files if 1 file: name.fld. Output: name_wss.fld
33  format of files if n checkpoint files: name_tn.fld. Output: name_tn_wss.fld, n=0,1,2...
34  */
35 
36 
37  if (nfiles == 1)
38  {
39  infiles[0] = argv[3];
40  string basename = argv[3];
41  basename = basename.substr(0, basename.find_last_of("."));
42  stringstream filename2;
43  filename2 << basename << "_wss.fld";
44  filename2 >> outfiles[0];
45  }
46  else {
47 
48  string basename = argv[3];
49  basename = basename.substr(basename.find_last_of("t")+1, basename.find_last_of(".")-basename.find_last_of("t"));
50  stringstream filename3;
51  filename3 << basename;
52  filename3 >> nStart;
53 
54  for (i = 0; i< nfiles; ++i)
55  {
56  basename = argv[3];
57  string extension = ".fld";
58  basename = basename.substr(0, basename.find_first_of("_"));
59  stringstream filename, filename2;
60  filename << basename << "_t" << i + nStart << extension;
61  filename >> infiles[i];
62  filename2 << basename << "_t" << i +nStart << "_wss.fld";
63  filename2 >> outfiles[i];
64  }
65  }
66 
67  argv[argc - 1] = argv[argc - 2];
68  argv[argc - 3] = argv[argc - 2];
69 
70  LibUtilities::SessionReaderSharedPtr vSession
71  = LibUtilities::SessionReader::CreateInstance(argc, argv);
72 
73  //----------------------------------------------
74  // Read in mesh from input file
75  string meshfile(argv[argc-4]);
76  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);
77  //----------------------------------------------
78 
79  //----------------------------------------------
80  //Get viscosity from file
81  NekDouble m_kinvis;
82  m_kinvis = vSession->GetParameter("Kinvis");
83  //----------------------------------------------
84 
85  //----------------------------------------------
86  // Import first field file.
87  string fieldfile(argv[3]);
88  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
89  vector<vector<NekDouble> > fielddata;
90  LibUtilities::Import(fieldfile,fielddef,fielddata);
91  //----------------------------------------------
92 
93  //----------------------------------------------
94  // Define Expansion
95  int expdim = graphShPt->GetMeshDimension();
96  int nfields = fielddef[0]->m_fields.size()-1;
97  int addfields = (nfields == 3)? 4:3;
98  int nstress = (nfields == 3)? 6:3;
99  Array<OneD, MultiRegions::ExpListSharedPtr> vel(nfields), shear(addfields);
100  MultiRegions::AssemblyMapCGSharedPtr m_locToGlobalMap;
101 
102  switch(expdim)
103  {
104  case 1:
105  {
106  ASSERTL0(false,"Expansion dimension not recognised");
107  }
108  break;
109  case 2:
110  {
111 
112  ASSERTL0(false,"Not implemented in 2D");
113 
114  i=0;
115  MultiRegions::ContField2DSharedPtr cont2D;
116  cont2D = MemoryManager<MultiRegions::ContField2D>
117  ::AllocateSharedPtr(vSession,graphShPt,vSession->GetVariable(i));
118  vel[0] = cont2D;
119 
120  for(i = 1; i < nfields; ++i)
121  {
122  vel[i] = MemoryManager<MultiRegions::ContField2D>
123  ::AllocateSharedPtr(*cont2D);
124  }
125  }
126  break;
127  case 3:
128  {
129  MultiRegions::ContField3DSharedPtr firstfield =
130  MemoryManager<MultiRegions::ContField3D>
131  ::AllocateSharedPtr(vSession, graphShPt,
132  vSession->GetVariable(0));
133 
134  m_locToGlobalMap = firstfield->GetLocalToGlobalMap();
135 
136  vel[0] = firstfield;
137  for(i = 1; i < nfields; ++i)
138  {
139  vel[i] = MemoryManager<MultiRegions::ContField3D>
140  ::AllocateSharedPtr(*firstfield, graphShPt,
141  vSession->GetVariable(i));
142  }
143 
144  for(i = 0; i < addfields; ++i)
145  {
146  shear[i] = MemoryManager<MultiRegions::ContField3D>
147  ::AllocateSharedPtr(*firstfield, graphShPt,
148  vSession->GetVariable(0));
149  }
150  }
151  break;
152  default:
153  ASSERTL0(false,"Expansion dimension not recognised");
154  break;
155  }
156  //----------------------------------------------
157 
158 
159  // Define arrays for WSS (global array), stress components and gradients (local arrays)
160  int n, cnt, elmtid, nq, offset, nt, boundary, nfq;
161  Array<OneD, Array<OneD, NekDouble> > grad(nfields*nfields);
162  Array<OneD, Array<OneD, NekDouble> > stress(nstress), fstress(nstress);
163  Array<OneD, Array<OneD, NekDouble> > values(addfields);
164 
165  // Set up mapping from Boundary condition to element details.
166  StdRegions::StdExpansionSharedPtr elmt;
167  StdRegions::StdExpansion2DSharedPtr bc;
168  Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
169  Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> > BndExp(addfields);
170 
171 
172  for (int fileNo = 0; fileNo < nfiles ; ++fileNo)
173  {
174  //----------------------------------------------
175  // Import field file.
176  fielddef.clear();
177  fielddata.clear();
178  LibUtilities::Import(infiles[fileNo],fielddef,fielddata);
179  //----------------------------------------------
180 
181  //----------------------------------------------
182  // Copy data from field file
183  for(j = 0; j < nfields; ++j)
184  {
185  for(i = 0; i < fielddata.size(); ++i)
186  {
187  vel[j]->ExtractDataToCoeffs(fielddef[i],
188  fielddata[i],
189  fielddef[i]->m_fields[j],
190  vel[j]->UpdateCoeffs());
191  }
192 
193  vel[j]->BwdTrans(vel[j]->GetCoeffs(),vel[j]->UpdatePhys());
194  }
195 
196  for(j = 0; j < addfields; ++j)
197  {
198  for(i = 0; i < fielddata.size(); ++i)
199  {
200  shear[j]->ExtractDataToCoeffs(fielddef[i],
201  fielddata[i],
202  fielddef[i]->m_fields[0],
203  shear[j]->UpdateCoeffs());
204 
205  }
206 
207  shear[j]->BwdTrans(shear[j]->GetCoeffs(),shear[j]->UpdatePhys());
208  }
209 
210  //----------------------------------------------
211 
212  //----------------------------------------------
213  // Compute WSS for each element on boundary
214  // Define total quadrature points, and velocity fields.
215  nt = shear[0]->GetNpoints();
216  Array<OneD, const NekDouble> U(nt),V(nt),W(nt);
217 
218  for(j = 0; j < addfields; ++j)
219  {
220  values[j] = Array<OneD, NekDouble>(nt);
221  }
222 
223 
224  shear[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
225 
226  //get boundary expansions for each field
227  for(j = 0; j < addfields; ++j)
228  {
229  BndExp[j] = shear[j]->GetBndCondExpansions();
230  }
231 
232  // loop over the types of boundary conditions
233  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
234  {
235  // identify boundary which the user wanted
236  if(n == surfID)
237  {
238  for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
239  {
240  // find element and face of this expansion.
241  elmtid = BoundarytoElmtID[cnt];
242  elmt = shear[0]->GetExp(elmtid);
243  nq = elmt->GetTotPoints();
244  offset = shear[0]->GetPhys_Offset(elmtid);
245 
246  // Initialise local arrays for the velocity gradients, and stress components
247  // size of total number of quadrature points for each element (hence local).
248  for(j = 0; j < nfields*nfields; ++j)
249  {
250  grad[j] = Array<OneD, NekDouble>(nq);
251  }
252 
253  for(j = 0; j < nstress; ++j)
254  {
255  stress[j] = Array<OneD, NekDouble>(nq);
256  }
257 
258 
259  if(nfields == 2)
260  {
261  //Not implemented in 2D.
262  }
263  else
264  {
265  // Get face 2D expansion from element expansion
266  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
267  nfq= bc->GetTotPoints();
268 
269  //identify boundary of element looking at.
270  boundary = BoundarytoTraceID[cnt];
271 
272  //Get face normals
273  const SpatialDomains::GeomFactorsSharedPtr m_metricinfo=bc->GetMetricInfo();
274 
275  const Array<OneD, const Array<OneD, NekDouble> > normals
276  = elmt->GetFaceNormal(boundary);
277 
278  // initialise arrays
279  for(j = 0; j < nstress; ++j)
280  {
281  fstress[j] = Array<OneD, NekDouble>(nfq);
282  }
283  Array<OneD, NekDouble> values2(nfq), S(nfq), Sx(nfq), Sy(nfq), Sz(nfq);
284 
285  //Extract Velocities
286  U = vel[0]->GetPhys() + offset;
287  V = vel[1]->GetPhys() + offset;
288  W = vel[2]->GetPhys() + offset;
289 
290  //Compute gradients (velocity correction scheme method)
291  elmt->PhysDeriv(U,grad[0],grad[1],grad[2]);
292  elmt->PhysDeriv(V,grad[3],grad[4],grad[5]);
293  elmt->PhysDeriv(W,grad[6],grad[7],grad[8]);
294 
295  //Compute stress component terms
296  // t_xx = 2.mu.Ux
297  Vmath::Smul (nq,(2*m_kinvis),grad[0],1,stress[0],1);
298  // tyy = 2.mu.Vy
299  Vmath::Smul (nq,(2*m_kinvis),grad[4],1,stress[1],1);
300  // tzz = 2.mu.Wz
301  Vmath::Smul (nq,(2*m_kinvis),grad[8],1,stress[2],1);
302  // txy = mu.(Uy+Vx)
303  Vmath::Vadd (nq,grad[1],1,grad[3],1,stress[3],1);
304  Vmath::Smul (nq,m_kinvis,stress[3],1,stress[3],1);
305  // txz = mu.(Uz+Wx)
306  Vmath::Vadd (nq,grad[2],1,grad[6],1,stress[4],1);
307  Vmath::Smul (nq,m_kinvis,stress[4],1,stress[4],1);
308  // tyz = mu.(Vz+Wy)
309  Vmath::Vadd (nq,grad[5],1,grad[7],1,stress[5],1);
310  Vmath::Smul (nq,m_kinvis,stress[5],1,stress[5],1);
311 
312  // Get face stress values.
313  for(j = 0; j < nstress; ++j)
314  {
315  elmt->GetFacePhysVals(boundary,bc,stress[j],fstress[j]);
316  }
317 
318  //calcuate wss, and update velocity coefficients in the elemental boundary expansion
319  for (j = 0; j< addfields; j++)
320  {
321  values[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[j][n]->GetCoeff_Offset(i);
322  }
323 
324  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
325  {
326  // Sx
327  Vmath::Vvtvvtp(nfq,normals[0],1,fstress[0],1,
328  normals[1],1,fstress[3],1,Sx,1);
329  Vmath::Vvtvp (nfq,normals[2],1,fstress[4],1,Sx,1,Sx,1);
330 
331  // Sy
332  Vmath::Vvtvvtp(nfq,normals[0],1,fstress[3],1,
333  normals[1],1,fstress[1],1,Sy,1);
334  Vmath::Vvtvp (nfq,normals[2],1,fstress[5],1,Sy,1,Sy,1);
335 
336  // Sz
337  Vmath::Vvtvvtp(nfq,normals[0],1,fstress[4],1,
338  normals[1],1,fstress[5],1,Sz,1);
339  Vmath::Vvtvp (nfq,normals[2],1,fstress[2],1,Sz,1,Sz,1);
340  }
341  else
342  {
343  // Sx
344  Vmath::Svtsvtp(nfq,normals[0][0],fstress[0],1,
345  normals[1][0],fstress[3],1,Sx,1);
346  Vmath::Svtvp(nfq,normals[2][0],fstress[4],1,Sx,1,Sx,1);
347 
348  // Sy
349  Vmath::Svtsvtp(nfq,normals[0][0],fstress[3],1,
350  normals[1][0],fstress[1],1,Sy,1);
351  Vmath::Svtvp(nfq,normals[2][0],fstress[5],1,Sy,1,Sy,1);
352 
353  // Sz
354  Vmath::Svtsvtp(nfq,normals[0][0],fstress[4],1,
355  normals[1][0],fstress[5],1,Sz,1);
356  Vmath::Svtvp(nfq,normals[2][0],fstress[2],1,Sz,1,Sz,1);
357  }
358 
359  // T = T - (T.n)n
360  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
361  {
362  Vmath::Vvtvvtp(nfq,normals[0],1,Sx,1,
363  normals[1],1, Sy,1,values2,1);
364  Vmath::Vvtvp (nfq,normals[2],1, Sz,1,values2,1,values2,1);
365 
366  // values2 = (T.n)n
367  Vmath::Smul(nfq, -1.0, values2, 1, values2, 1);
368 
369  //Sx
370  Vmath::Vvtvp(nfq,normals[0],1, values2,1,Sx,1,Sx,1);
371  bc->FwdTrans(Sx, values[0]);
372 
373  //Sy
374  Vmath::Vvtvp(nfq,normals[1],1, values2,1,Sy,1,Sy,1);
375  bc->FwdTrans(Sy, values[1]);
376 
377  //Sz
378  Vmath::Vvtvp(nfq,normals[2],1, values2,1,Sz,1,Sz,1);
379  bc->FwdTrans(Sz, values[2]);
380 
381  }
382  else
383  {
384  Vmath::Svtsvtp(nfq,normals[0][0],Sx,1,
385  normals[1][0],Sy,1,values2,1);
386  Vmath::Svtvp(nfq,normals[2][0],Sz,1,values2,1,values2,1);
387 
388  Vmath::Smul(nfq, -1.0, values2, 1, values2, 1);
389 
390  //Sx
391  Vmath::Svtvp(nfq,normals[0][0],values2,1,Sx,1,Sx,1);
392  bc->FwdTrans(Sx, values[0]);
393 
394  //Sy
395  Vmath::Svtvp(nfq,normals[1][0],values2,1,Sy,1,Sy,1);
396  bc->FwdTrans(Sy, values[1]);
397 
398  //Sz
399  Vmath::Svtvp(nfq,normals[2][0],values2,1,Sz,1,Sz,1);
400  bc->FwdTrans(Sz, values[2]);
401  }
402 
403 
404  //Tw
405  Vmath::Vvtvvtp(nfq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, S, 1);
406  Vmath::Vvtvp(nfq, Sz, 1, Sz, 1, S, 1, S, 1);
407  Vmath::Vsqrt(nfq, S, 1, S, 1);
408  bc->FwdTrans(S, values[3]);
409  }
410  }
411  }
412  else
413  {
414  cnt += BndExp[0][n]->GetExpSize();
415  }
416  }
417 
418 
419  for(j = 0; j < addfields; ++j)
420  {
421  int ncoeffs = shear[j]->GetNcoeffs();
422  Array<OneD, NekDouble> output(ncoeffs);
423 
424  output=shear[j]->UpdateCoeffs();
425 
426  int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
427  Array<OneD, NekDouble> outarray(nGlobal,0.0);
428 
429  int bndcnt=0;
430 
431  const Array<OneD,const int>& map = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsMap();
432  NekDouble sign;
433 
434  for(i = 0; i < BndExp[j].num_elements(); ++i)
435  {
436  if(i==surfID)
437  {
438  const Array<OneD,const NekDouble>& coeffs = BndExp[j][i]->GetCoeffs();
439  for(k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
440  {
441  sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
442  outarray[map[bndcnt++]] = sign * coeffs[k];
443  }
444  }
445  else
446  {
447  bndcnt += BndExp[j][i]->GetNcoeffs();
448  }
449  }
450  m_locToGlobalMap->GlobalToLocal(outarray,output);
451  }
452 
453 
454  //-----------------------------------------------
455  // Write solution to file with additional computed fields
456  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
457  = shear[0]->GetFieldDefinitions();
458  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
459 
460  vector<string > outname;
461 
462  if(addfields == 3)
463  {
464  // Not implemented in 2D.
465  }
466  else
467  {
468  outname.push_back("Tx");
469  outname.push_back("Ty");
470  outname.push_back("Tz");
471  outname.push_back("Tw");
472  }
473 
474  for(j = 0; j < nfields + addfields; ++j)
475  {
476  for(i = 0; i < FieldDef.size(); ++i)
477  {
478  if (j < nfields)
479  {
480  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
481  vel[j]->AppendFieldData(FieldDef[i], FieldData[i]);
482  }
483  else
484  {
485  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
486  shear[j-nfields]->AppendFieldData(FieldDef[i], FieldData[i]);
487  }
488  }
489  }
490  LibUtilities::Write(outfiles[fileNo], FieldDef, FieldData);
491 
492  FieldDef.clear();
493  FieldData.clear();
494  //-----------------------------------------------
495  }
496 
497  return 0;
498 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Vmath::Vvtvp
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
Nektar::MultiRegions::ContField2DSharedPtr
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:293
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
Definition: FieldIO.cpp:115
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
Nektar::StdRegions::StdExpansion2DSharedPtr
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
Definition: StdExpansion2D.h:203
Vmath::Smul
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
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SpatialDomains::GeomFactorsSharedPtr
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:523
Vmath::Svtsvtp
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:577
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:81
Nektar::MultiRegions::ContField3DSharedPtr
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:191
Nektar::MultiRegions::AssemblyMapCGSharedPtr
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:52
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Vmath::Vadd
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
Generated on Sat Jun 4 2016 00:02:27 for Nektar++ by   doxygen 1.8.8