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