Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ExtractSurface3DCFS.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExtractSurface3DCFS.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: Extract 2D surfaces from 3D file and output relevant quantities
33 // for compressible flow solver.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <cstdio>
38 #include <cstdlib>
39 #include <string>
40 #include <iostream>
41 #include <iomanip>
42 
43 #include <MultiRegions/ExpList.h>
44 #include <MultiRegions/ExpList1D.h>
45 #include <MultiRegions/ExpList2D.h>
46 #include <MultiRegions/ExpList3D.h>
48 
50 #include <LocalRegions/MatrixKey.h>
52 #include <LocalRegions/Expansion.h>
53 
59 
63 
65 
66 using namespace std;
67 using namespace Nektar;
68 
69 int main(int argc, char *argv[])
70 {
71  string fname = std::string(argv[2]);
72  int fdot = fname.find_last_of('.');
73  if (fdot != std::string::npos)
74  {
75  string ending = fname.substr(fdot);
76 
77  // If .chk or .fld we exchange the extension in the output file.
78  // For all other files (e.g. .bse) we append the extension to avoid
79  // conflicts.
80  if (ending == ".chk" || ending == ".fld")
81  {
82  fname = fname.substr(0,fdot);
83  }
84  }
85 
86  fname = fname + ".txt";
87 
88  int cnt;
89  int id1, id2;
90  int i, j, n, e, b;
91  Array<OneD, NekDouble> auxArray;
92 
93  int nBndEdgePts, nBndEdges, nBndRegions;
94 
95  if (argc < 3)
96  {
97  fprintf(stderr,
98  "Usage: ExtractSurface3DCFS meshfile fieldFile\n");
99  fprintf(stderr,
100  "Extracts a surface from a 3D fld file"
101  "(only for CompressibleFlowSolver and purely 3D .fld files)\n");
102  exit(1);
103  }
104 
106  = LibUtilities::SessionReader::CreateInstance(3, argv);
107 
108  std::string m_ViscosityType;
109 
110  NekDouble m_gamma;
111  NekDouble m_pInf;
115  NekDouble m_wInf;
116  NekDouble m_gasConstant;
118  NekDouble m_mu;
119 
120  int m_spacedim = 3;
121  int nDimensions = m_spacedim;
122  int phys_offset;
123 
124  // Get gamma parameter from session file.
125  ASSERTL0(vSession->DefinesParameter("Gamma"),
126  "Compressible flow sessions must define a Gamma parameter.");
127  vSession->LoadParameter("Gamma", m_gamma, 1.4);
128 
129  // Get E0 parameter from session file.
130  ASSERTL0(vSession->DefinesParameter("pInf"),
131  "Compressible flow sessions must define a pInf parameter.");
132  vSession->LoadParameter("pInf", m_pInf, 101325);
133 
134  // Get rhoInf parameter from session file.
135  ASSERTL0(vSession->DefinesParameter("rhoInf"),
136  "Compressible flow sessions must define a rhoInf parameter.");
137  vSession->LoadParameter("rhoInf", m_rhoInf, 1.225);
138 
139  // Get uInf parameter from session file.
140  ASSERTL0(vSession->DefinesParameter("uInf"),
141  "Compressible flow sessions must define a uInf parameter.");
142  vSession->LoadParameter("uInf", m_uInf, 0.1);
143 
144  // Get vInf parameter from session file.
145  if (m_spacedim == 2 || m_spacedim == 3)
146  {
147  ASSERTL0(vSession->DefinesParameter("vInf"),
148  "Compressible flow sessions must define a vInf parameter"
149  "for 2D/3D problems.");
150  vSession->LoadParameter("vInf", m_vInf, 0.0);
151  }
152 
153  // Get wInf parameter from session file.
154  if (m_spacedim == 3)
155  {
156  ASSERTL0(vSession->DefinesParameter("wInf"),
157  "Compressible flow sessions must define a wInf parameter"
158  "for 3D problems.");
159  vSession->LoadParameter("wInf", m_wInf, 0.0);
160  }
161 
162  vSession->LoadParameter ("GasConstant", m_gasConstant, 287.058);
163  vSession->LoadParameter ("Twall", m_Twall, 300.15);
164  vSession->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
165  vSession->LoadParameter ("mu", m_mu, 1.78e-05);
166 
167  //--------------------------------------------------------------------------
168  // Read in mesh from input file
169  string meshfile(argv[1]);
171  SpatialDomains::MeshGraph::Read(vSession);
172  //--------------------------------------------------------------------------
173 
174  //--------------------------------------------------------------------------
175  // Import field file
176  string fieldFile(argv[2]);
177  vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
178  vector<vector<NekDouble> > fieldData;
179 
180  LibUtilities::Import(fieldFile, fieldDef, fieldData);
181  //--------------------------------------------------------------------------
182 
183  //--------------------------------------------------------------------------
184  // Set up Expansion information
185  vector< vector<LibUtilities::PointsType> > pointsType;
186  for (i = 0; i < fieldDef.size(); ++i)
187  {
188  vector<LibUtilities::PointsType> ptype;
189  for (j = 0; j < 3; ++j)
190  {
191  ptype.push_back(LibUtilities::ePolyEvenlySpaced);
192  }
193  pointsType.push_back(ptype);
194  }
195  graphShPt->SetExpansions(fieldDef, pointsType);
196 
197  //--------------------------------------------------------------------------
198 
199 
200  //--------------------------------------------------------------------------
201  // Define Expansion
202  int nfields = fieldDef[0]->m_fields.size();
205 
206  for(i = 0; i < pFields.num_elements(); i++)
207  {
208  pFields[i] = MemoryManager<MultiRegions
209  ::DisContField3D>::AllocateSharedPtr(vSession, graphShPt,
210  vSession->GetVariable(i));
211  }
212 
215  ::AllocateSharedPtr(vSession, graphShPt);
216 
217  Exp[0] = Exp3D;
218 
219  for (i = 1; i < nfields; ++i)
220  {
222  ::AllocateSharedPtr(*Exp3D);
223  }
224 
225  // Count of the point on the surface
226  int nSurfacePts = 0;
227  if (pFields[0]->GetBndCondExpansions().num_elements())
228  {
229  nSurfacePts = 0;
230  cnt = 0;
231  nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
232  for (b = 0; b < nBndRegions; ++b)
233  {
234  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
235  for (e = 0; e < nBndEdges; ++e)
236  {
237  nBndEdgePts = pFields[0]->
238  GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
239 
240  if (pFields[0]->GetBndConditions()[b]->
241  GetUserDefined() == "WallViscous" ||
242  pFields[0]->GetBndConditions()[b]->
243  GetUserDefined() == "WallAdiabatic" ||
244  pFields[0]->GetBndConditions()[b]->
245  GetUserDefined() == "Wall")
246  {
247  nSurfacePts += nBndEdgePts;
248  }
249  }
250  }
251  }
252 
253 
254  int nSolutionPts = pFields[0]->GetNpoints();
255  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
256  int nElements = pFields[0]->GetExpSize();
257 
258  Array<OneD, NekDouble> tmp(nSolutionPts, 0.0);
259 
260  Array<OneD, NekDouble> x(nSolutionPts);
261  Array<OneD, NekDouble> y(nSolutionPts);
262  Array<OneD, NekDouble> z(nSolutionPts);
263 
264  Array<OneD, NekDouble> traceX(nTracePts);
265  Array<OneD, NekDouble> traceY(nTracePts);
266  Array<OneD, NekDouble> traceZ(nTracePts);
267 
268  Array<OneD, NekDouble> surfaceX(nSurfacePts);
269  Array<OneD, NekDouble> surfaceY(nSurfacePts);
270  Array<OneD, NekDouble> surfaceZ(nSurfacePts);
271 
272  pFields[0]->GetCoords(x, y, z);
273 
274  pFields[0]->ExtractTracePhys(x, traceX);
275  pFields[0]->ExtractTracePhys(y, traceY);
276  pFields[0]->ExtractTracePhys(z, traceZ);
277  //--------------------------------------------------------------------------
278 
279  //--------------------------------------------------------------------------
280  // Copy data from field file
281  Array<OneD, Array<OneD, NekDouble> > uFields(nfields);
282  Array<OneD, Array<OneD, NekDouble> > traceFields(nfields);
283  Array<OneD, Array<OneD, NekDouble> > surfaceFields(nfields);
284 
285  // Extract the physical values of the solution at the boundaries
286  for (j = 0; j < nfields; ++j)
287  {
288  uFields[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
289  traceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
290  surfaceFields[j] = Array<OneD, NekDouble>(nSurfacePts, 0.0);
291 
292 
293  for (i = 0; i < fieldData.size(); ++i)
294  {
295  Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
296  fieldDef[i]->m_fields[j],
297  Exp[j]->UpdateCoeffs());
298  }
299  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
300  Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
301  pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
302  }
303 
304 
305  //Fields to add in the output file
306  int nfieldsAdded = 34;
307  Array<OneD, Array<OneD, NekDouble> > traceFieldsAdded(nfieldsAdded);
308  Array<OneD, Array<OneD, NekDouble> > surfaceFieldsAdded(nfieldsAdded);
309 
310  for (j = 0; j < nfieldsAdded; ++j)
311  {
312  traceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
313  surfaceFieldsAdded[j] = Array<OneD, NekDouble>(nSurfacePts, 0.0);
314  }
315 
316  /******** Evaluation of normals and tangents on the trace *****************
317  * nx -> traceFieldsAdded[0];
318  * ny -> traceFieldsAdded[1];
319  * nz -> traceFieldsAdded[2];
320  * bx -> traceFieldsAdded[3];
321  * by -> traceFieldsAdded[4];
322  * bz -> traceFieldsAdded[5];
323  * tx -> traceFieldsAdded[6];
324  * ty -> traceFieldsAdded[7];
325  * tz -> traceFieldsAdded[8];
326 
327  ***************************************************************************/
328 
329  Array<OneD, Array<OneD, NekDouble> > m_traceNormals (nDimensions);
330  for(i = 0; i < nDimensions; ++i)
331  {
332  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
333  }
334  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
335 
336  Array<OneD, Array<OneD, NekDouble> > m_traceTangents (nDimensions);
337  Array<OneD, Array<OneD, NekDouble> > m_traceBinormals (nDimensions);
338  Array<OneD, Array<OneD, NekDouble> > h (nDimensions);
339  Array<OneD, NekDouble > tmpNorm (nTracePts, 1.0);
340  Array<OneD, NekDouble > NormH (nTracePts, 0.0);
341  Array<OneD, NekDouble > tmpTrace (nTracePts, 0.0);
342 
343 
344  for(i = 0; i < nDimensions; ++i)
345  {
346  m_traceTangents[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
347  m_traceBinormals[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
348  h[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
349  }
350 
351  // Normals
352 
353  // nx
354  Vmath::Vcopy(nTracePts,
355  &m_traceNormals[0][0], 1,
356  &traceFieldsAdded[0][0], 1);
357 
358  // ny
359  Vmath::Vcopy(nTracePts,
360  &m_traceNormals[1][0], 1,
361  &traceFieldsAdded[1][0], 1);
362 
363  // nz
364  Vmath::Vcopy(nTracePts,
365  &m_traceNormals[2][0], 1,
366  &traceFieldsAdded[2][0], 1);
367 
368 
369  // Tangents and Binormals
370  // h1
371  Vmath::Vadd(nTracePts,
372  &m_traceNormals[0][0], 1,
373  &tmpNorm[0], 1,
374  &h[0][0], 1);
375  // h2
376  Vmath::Vcopy(nTracePts,
377  &m_traceNormals[1][0], 1,
378  &h[1][0], 1);
379  // h3
380  Vmath::Vcopy(nTracePts,
381  &m_traceNormals[2][0], 1,
382  &h[2][0], 1);
383 
384  // Norm of h
385  for (i = 0; i < m_spacedim; i++)
386  {
387  Vmath::Vvtvp (nTracePts, &h[i][0], 1, &h[i][0], 1,
388  &NormH[0],1, &NormH[0],1);
389  }
390 
391  //b1
392  Vmath::Vmul(nTracePts,
393  &h[0][0], 1,
394  &h[1][0], 1,
395  &tmpTrace[0],1);
396 
397  Vmath::Vdiv(nTracePts,
398  &tmpTrace[0],1,
399  &NormH[0], 1,
400  &tmpTrace[0],1);
401 
402  Vmath::Smul(nTracePts, -2.0,
403  &tmpTrace[0], 1,
404  &m_traceBinormals[0][0], 1);
405 
406  Vmath::Vcopy(nTracePts,
407  &m_traceBinormals[0][0], 1,
408  &traceFieldsAdded[3][0], 1);
409 
410 
411  //b2
412  Vmath::Vmul(nTracePts,
413  &h[1][0], 1,
414  &h[1][0], 1,
415  &tmpTrace[0],1);
416 
417  Vmath::Vdiv(nTracePts,
418  &tmpTrace[0],1,
419  &NormH[0], 1,
420  &tmpTrace[0],1);
421 
422  Vmath::Smul(nTracePts, -2.0,
423  &tmpTrace[0], 1,
424  &tmpTrace[0], 1);
425 
426  Vmath::Vadd(nTracePts,
427  &tmpTrace[0], 1,
428  &tmpNorm[0], 1,
429  &m_traceBinormals[1][0], 1);
430 
431  Vmath::Vcopy(nTracePts,
432  &m_traceBinormals[1][0], 1,
433  &traceFieldsAdded[4][0], 1);
434 
435 
436  //b3
437  Vmath::Vmul(nTracePts,
438  &h[1][0], 1,
439  &h[2][0], 1,
440  &tmpTrace[0],1);
441 
442  Vmath::Vdiv(nTracePts,
443  &tmpTrace[0],1,
444  &NormH[0], 1,
445  &tmpTrace[0],1);
446 
447  Vmath::Smul(nTracePts, -2.0,
448  &tmpTrace[0], 1,
449  &m_traceBinormals[2][0], 1);
450 
451  Vmath::Vcopy(nTracePts,
452  &m_traceBinormals[2][0], 1,
453  &traceFieldsAdded[5][0], 1);
454 
455 
456  //t1
457  Vmath::Vmul(nTracePts,
458  &h[0][0], 1,
459  &h[2][0], 1,
460  &tmpTrace[0],1);
461 
462  Vmath::Vdiv(nTracePts,
463  &tmpTrace[0],1,
464  &NormH[0], 1,
465  &tmpTrace[0],1);
466 
467  Vmath::Smul(nTracePts, -2.0,
468  &tmpTrace[0], 1,
469  &m_traceTangents[0][0], 1);
470 
471  Vmath::Vcopy(nTracePts,
472  &m_traceTangents[0][0], 1,
473  &traceFieldsAdded[6][0], 1);
474 
475  //t2
476  Vmath::Vcopy(nTracePts,
477  &m_traceBinormals[2][0], 1,
478  &m_traceTangents[1][0], 1);
479 
480  Vmath::Vcopy(nTracePts,
481  &m_traceTangents[1][0], 1,
482  &traceFieldsAdded[7][0], 1);
483 
484 
485  //t3
486  Vmath::Vmul(nTracePts,
487  &h[2][0], 1,
488  &h[2][0], 1,
489  &tmpTrace[0],1);
490 
491  Vmath::Vdiv(nTracePts,
492  &tmpTrace[0],1,
493  &NormH[0], 1,
494  &tmpTrace[0],1);
495 
496  Vmath::Smul(nTracePts, -2.0,
497  &tmpTrace[0], 1,
498  &tmpTrace[0], 1);
499 
500  Vmath::Vadd(nTracePts,
501  &tmpTrace[0], 1,
502  &tmpNorm[0], 1,
503  &m_traceTangents[2][0], 1);
504 
505  Vmath::Vcopy(nTracePts,
506  &m_traceTangents[2][0], 1,
507  &traceFieldsAdded[8][0], 1);
508 
509 
510  /******** Evaluation of the pressure ***************************************
511  * P = (E-1/2.*rho.*((rhou./rho).^2+(rhov./rho).^2))*(gamma - 1);
512  * P -> traceFieldsAdded[9];
513  ***************************************************************************/
514 
515  Array<OneD, NekDouble> pressure(nSolutionPts, 0.0);
516  NekDouble gammaMinusOne = m_gamma - 1.0;
517 
518 
519  for (i = 0; i < m_spacedim; i++)
520  {
521  Vmath::Vmul(nSolutionPts,
522  &uFields[i + 1][0], 1,
523  &uFields[i + 1][0], 1,
524  &tmp[0],1);
525 
526 
527  Vmath::Smul(nSolutionPts, 0.5,
528  &tmp[0], 1,
529  &tmp[0], 1);
530 
531  Vmath::Vadd(nSolutionPts,
532  &pressure[0], 1,
533  &tmp[0], 1,
534  &pressure[0], 1);
535  }
536 
537  Vmath::Vdiv(nSolutionPts,
538  &pressure[0], 1,
539  &uFields[0][0], 1,
540  &pressure[0],1);
541 
542  Vmath::Vsub(nSolutionPts,
543  &uFields[nfields - 1][0], 1,
544  &pressure[0], 1,
545  &pressure[0],1);
546 
547  Vmath::Smul(nSolutionPts, gammaMinusOne,
548  &pressure[0], 1,
549  &pressure[0], 1);
550 
551  // Extract trace
552  pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[9]);
553 
554  /******** Evaluation of the temperature ************************************
555  * T = P/(R*rho);
556  * T -> traceFieldsAdded[10];
557  ***************************************************************************/
558 
559  Array<OneD, NekDouble> temperature(nSolutionPts, 0.0);
560 
561  Vmath::Vdiv(nSolutionPts,
562  &pressure[0], 1,
563  &uFields[0][0], 1,
564  &temperature[0],1);
565 
566  NekDouble GasConstantInv = 1.0/m_gasConstant;
567  Vmath::Smul(nSolutionPts, GasConstantInv,
568  &temperature[0], 1,
569  &temperature[0], 1);
570 
571  // Extract trace
572  pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[10]);
573 
574  /*** Evaluation of the temperature gradient in the normal direction ********
575  * DT_n -> traceFieldsAdded[11]
576  ***************************************************************************/
577 
578  Array<OneD, Array<OneD, NekDouble> > Dtemperature(nDimensions);
579  Array<OneD, Array<OneD, NekDouble> > traceDtemperature(nDimensions);
580 
581  for (i = 0; i < nDimensions; ++ i)
582  {
583  Dtemperature[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
584  traceDtemperature[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
585  }
586 
587  for (i = 0; i < nDimensions; ++ i)
588  {
589  for (n = 0; n < nElements; n++)
590  {
591  phys_offset = pFields[0]->GetPhys_Offset(n);
592 
593  pFields[i]->GetExp(n)->PhysDeriv(
594  i, temperature + phys_offset,
595  auxArray = Dtemperature[i] + phys_offset);
596  }
597  // Extract trace
598  pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
599  }
600 
601  for(i = 0; i < nDimensions; ++i)
602  {
603  Vmath::Vmul(nTracePts,
604  &m_traceNormals[i][0], 1,
605  &traceDtemperature[i][0], 1,
606  &tmp[0],1);
607 
608  Vmath::Vadd(nTracePts,
609  &traceFieldsAdded[11][0], 1,
610  &tmp[0], 1,
611  &traceFieldsAdded[11][0], 1);
612  }
613 
614  /*** Evaluation of the pressure gradient ***********************************
615  * DP_t -> traceFieldsAdded[12] tangent direction
616  * DP_b -> traceFieldsAdded[13] binormal direction
617  * DP_x -> traceFieldsAdded[14]
618  * DP_y -> traceFieldsAdded[15]
619  * DP_z -> traceFieldsAdded[16]
620  ***************************************************************************/
621 
622  Array<OneD, Array<OneD, NekDouble> > Dpressure(nDimensions);
623  Array<OneD, Array<OneD, NekDouble> > traceDpressure(nDimensions);
624 
625  for (i = 0; i < nDimensions; ++ i)
626  {
627  Dpressure[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
628  traceDpressure[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
629  }
630 
631  for (i = 0; i < nDimensions; ++ i)
632  {
633  for (n = 0; n < nElements; n++)
634  {
635  phys_offset = pFields[0]->GetPhys_Offset(n);
636 
637  pFields[i]->GetExp(n)->PhysDeriv(
638  i, pressure + phys_offset,
639  auxArray = Dpressure[i] + phys_offset);
640  }
641  // Extract trace
642  pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
643  }
644 
645  // Dp_t
646  for(i = 0; i < nDimensions; ++i)
647  {
648  Vmath::Vmul(nTracePts,
649  &m_traceTangents[i][0], 1,
650  &traceDpressure[i][0], 1,
651  &tmp[0],1);
652 
653  Vmath::Vadd(nTracePts,
654  &traceFieldsAdded[12][0], 1,
655  &tmp[0], 1,
656  &traceFieldsAdded[12][0], 1);
657  }
658 
659  // Dp_b
660  for(i = 0; i < nDimensions; ++i)
661  {
662  Vmath::Vmul(nTracePts,
663  &m_traceBinormals[i][0], 1,
664  &traceDpressure[i][0], 1,
665  &tmp[0],1);
666 
667  Vmath::Vadd(nTracePts,
668  &traceFieldsAdded[13][0], 1,
669  &tmp[0], 1,
670  &traceFieldsAdded[13][0], 1);
671  }
672 
673 
674  // Dp_x
675  Vmath::Vcopy(nTracePts,
676  &traceDpressure[0][0], 1,
677  &traceFieldsAdded[14][0], 1);
678 
679  // Dp_y
680  Vmath::Vcopy(nTracePts,
681  &traceDpressure[1][0], 1,
682  &traceFieldsAdded[15][0], 1);
683 
684  // Dp_z
685  Vmath::Vcopy(nTracePts,
686  &traceDpressure[2][0], 1,
687  &traceFieldsAdded[16][0], 1);
688 
689 
690  /** Evaluation of the velocity gradient in the cartesian directions
691  * Du_x: traceFieldsAdded[17]
692  * Du_y: traceFieldsAdded[18]
693  * Du_z: traceFieldsAdded[19]
694  * Dv_x: traceFieldsAdded[20]
695  * Dv_y: traceFieldsAdded[21]
696  * Dv_z: traceFieldsAdded[22]
697  * Dw_x: traceFieldsAdded[23]
698  * Dw_y: traceFieldsAdded[24]
699  * Dw_z: traceFieldsAdded[25]
700  **/
701 
702  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Dvelocity(nDimensions);
703  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > traceDvelocity(nDimensions);
704  Array<OneD, Array<OneD, NekDouble> > velocity(nDimensions);
705 
706  for (i = 0; i < nDimensions; ++ i)
707  {
708  Dvelocity[i] = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
709  traceDvelocity[i] = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
710  velocity[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
711 
712  Vmath::Vdiv(nSolutionPts, uFields[i+1], 1, uFields[0], 1,
713  velocity[i], 1);
714 
715  for (j = 0; j < nDimensions; ++j)
716  {
717  Dvelocity[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
718  traceDvelocity[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
719  }
720  }
721 
722  for (i = 0; i < nDimensions; ++i)
723  {
724  for (j = 0; j < nDimensions; ++j)
725  {
726  for (n = 0; n < nElements; n++)
727  {
728  phys_offset = pFields[0]->GetPhys_Offset(n);
729 
730  pFields[i]->GetExp(n)->PhysDeriv(
731  j, velocity[i] + phys_offset,
732  auxArray = Dvelocity[i][j] + phys_offset);
733  }
734 
735  // Extract trace
736  pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
737  }
738  }
739 
740  Vmath::Vcopy(nTracePts,
741  &traceDvelocity[0][0][0], 1,
742  &traceFieldsAdded[17][0], 1);
743  Vmath::Vcopy(nTracePts,
744  &traceDvelocity[0][1][0], 1,
745  &traceFieldsAdded[18][0], 1);
746  Vmath::Vcopy(nTracePts,
747  &traceDvelocity[0][2][0], 1,
748  &traceFieldsAdded[19][0], 1);
749  Vmath::Vcopy(nTracePts,
750  &traceDvelocity[1][0][0], 1,
751  &traceFieldsAdded[20][0], 1);
752  Vmath::Vcopy(nTracePts,
753  &traceDvelocity[1][1][0], 1,
754  &traceFieldsAdded[21][0], 1);
755  Vmath::Vcopy(nTracePts,
756  &traceDvelocity[1][2][0], 1,
757  &traceFieldsAdded[22][0], 1);
758  Vmath::Vcopy(nTracePts,
759  &traceDvelocity[2][0][0], 1,
760  &traceFieldsAdded[23][0], 1);
761  Vmath::Vcopy(nTracePts,
762  &traceDvelocity[2][1][0], 1,
763  &traceFieldsAdded[24][0], 1);
764  Vmath::Vcopy(nTracePts,
765  &traceDvelocity[2][2][0], 1,
766  &traceFieldsAdded[25][0], 1);
767 
768 
769  /*** Evaluation of shear stresses ******************************************
770  * tau_xx -> traceFieldsAdded[26]
771  * tau_yy -> traceFieldsAdded[27]
772  * tau_zz -> traceFieldsAdded[28]
773  * tau_xy -> traceFieldsAdded[29]
774  * tau_xz -> traceFieldsAdded[30]
775  * tau_yz -> traceFieldsAdded[31]
776  ***************************************************************************/
777 
778  // Stokes hypotesis
779  const NekDouble lambda = -2.0/3.0;
780 
781  // Auxiliary variables
782  Array<OneD, NekDouble > mu (nSolutionPts, 0.0);
783  Array<OneD, NekDouble > mu2 (nSolutionPts, 0.0);
784  Array<OneD, NekDouble > divVel(nSolutionPts, 0.0);
785 
786  // Variable viscosity through the Sutherland's law
787  if (m_ViscosityType == "Variable")
788  {
789  NekDouble mu_star = m_mu;
790  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
791  NekDouble ratio;
792 
793  for (int i = 0; i < nSolutionPts; ++i)
794  {
795  ratio = temperature[i] / T_star;
796  mu[i] = mu_star * ratio * sqrt(ratio) *
797  (T_star + 110.0) / (temperature[i] + 110.0);
798  }
799  }
800  else
801  {
802  Vmath::Fill(nSolutionPts, m_mu, &mu[0], 1);
803  }
804 
805  // Computing diagonal terms of viscous stress tensor
806  Array<OneD, Array<OneD, NekDouble> > temp(m_spacedim);
807  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
808 
809  // mu2 = 2 * mu
810  Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
811 
812  // Velocity divergence
813  Vmath::Vadd(nSolutionPts, &divVel[0], 1,
814  &Dvelocity[0][0][0], 1, &divVel[0], 1);
815  Vmath::Vadd(nSolutionPts, &divVel[0], 1,
816  &Dvelocity[1][1][0], 1, &divVel[0], 1);
817 
818  // Velocity divergence scaled by lambda * mu
819  Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
820  Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
821 
822  // Diagonal terms of viscous stress tensor (Sxx, Syy)
823  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
824  for (j = 0; j < m_spacedim; ++j)
825  {
826  temp[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
827  Sgg[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
828 
829  Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
830  &temp[j][0], 1);
831 
832  Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
833  }
834 
835  // Extra diagonal terms of viscous stress tensor
836  Array<OneD, NekDouble > Sxy(nSolutionPts, 0.0);
837  Array<OneD, NekDouble > Sxz(nSolutionPts, 0.0);
838  Array<OneD, NekDouble > Syz(nSolutionPts, 0.0);
839 
840  // Sxy = (du/dy + dv/dx)
841  Vmath::Vadd(nSolutionPts, &Dvelocity[0][1][0], 1,
842  &Dvelocity[1][0][0], 1, &Sxy[0], 1);
843 
844  // Sxz = (du/dz + dw/dx)
845  Vmath::Vadd(nSolutionPts, &Dvelocity[0][2][0], 1,
846  &Dvelocity[2][0][0], 1, &Sxz[0], 1);
847 
848  // Syz = (dv/dz + dw/dy)
849  Vmath::Vadd(nSolutionPts, &Dvelocity[1][2][0], 1,
850  &Dvelocity[2][1][0], 1, &Syz[0], 1);
851 
852  // Sxy = mu * (du/dy + dv/dx)
853  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
854 
855  // Sxz = mu * (du/dy + dv/dx)
856  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
857 
858  // Syz = mu * (du/dy + dv/dx)
859  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
860 
861 
862 
863  pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[26]);
864  pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[27]);
865  pFields[0]->ExtractTracePhys(Sgg[2], traceFieldsAdded[28]);
866  pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[29]);
867  pFields[0]->ExtractTracePhys(Sxz, traceFieldsAdded[30]);
868  pFields[0]->ExtractTracePhys(Syz, traceFieldsAdded[31]);
869 
870  /*** Evaluation of dinamic viscosity ***************************************
871  * mu -> traceFieldsAdded[32]
872  ***************************************************************************/
873 
874  pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[32]);
875 
876  /*** Evaluation of Mach number *********************************************
877  * M -> traceFieldsAdded[33]
878  ***************************************************************************/
879  NekDouble gamma = m_gamma;
880 
881  // Speed of sound
882  Array<OneD, NekDouble> soundspeed(nSolutionPts, 0.0);
883 
884  Vmath::Vdiv (nSolutionPts, pressure, 1, uFields[0], 1, soundspeed, 1);
885  Vmath::Smul (nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
886  Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
887 
888  // Mach
889  Array<OneD, NekDouble> mach(nSolutionPts, 0.0);
890 
891  for (int i = 0; i < m_spacedim; ++i)
892  {
893  Vmath::Vvtvp(nSolutionPts,
894  uFields[i + 1], 1,
895  uFields[i + 1], 1,
896  mach, 1, mach, 1);
897  }
898 
899  Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
900  Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
901  Vmath::Vsqrt(nSolutionPts, mach, 1, mach, 1);
902  Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
903 
904  pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[33]);
905 
906 
907 
908  /**************************************************************************/
909  // Extract coordinates
910 
911  if (pFields[0]->GetBndCondExpansions().num_elements())
912  {
913  id1 = 0;
914  cnt = 0;
915  nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
916  for (b = 0; b < nBndRegions; ++b)
917  {
918  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
919  for (e = 0; e < nBndEdges; ++e)
920  {
921  nBndEdgePts = pFields[0]->
922  GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
923 
924  id2 = pFields[0]->GetTrace()->
925  GetPhys_Offset(pFields[0]->GetTraceMap()->
926  GetBndCondTraceToGlobalTraceMap(cnt++));
927 
928  if (pFields[0]->GetBndConditions()[b]->
929  GetUserDefined() == "WallViscous" ||
930  pFields[0]->GetBndConditions()[b]->
931  GetUserDefined() == "WallAdiabatic" ||
932  pFields[0]->GetBndConditions()[b]->
933  GetUserDefined() == "Wall")
934  {
935 
936  Vmath::Vcopy(nBndEdgePts, &traceX[id2], 1,
937  &surfaceX[id1], 1);
938 
939  Vmath::Vcopy(nBndEdgePts, &traceY[id2], 1,
940  &surfaceY[id1], 1);
941 
942  Vmath::Vcopy(nBndEdgePts, &traceZ[id2], 1,
943  &surfaceZ[id1], 1);
944 
945  id1 += nBndEdgePts;
946  }
947  }
948  }
949  }
950 
951  // Extract fields
952  if (pFields[0]->GetBndCondExpansions().num_elements())
953  {
954 
955  for (j = 0; j < nfields; ++j)
956  {
957  cout << "field " << j << endl;
958 
959  id1 = 0;
960  cnt = 0;
961  nBndRegions = pFields[j]->GetBndCondExpansions().num_elements();
962  for (b = 0; b < nBndRegions; ++b)
963  {
964  nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
965  for (e = 0; e < nBndEdges; ++e)
966  {
967  nBndEdgePts = pFields[j]->
968  GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
969 
970  id2 = pFields[j]->GetTrace()->
971  GetPhys_Offset(pFields[j]->GetTraceMap()->
972  GetBndCondTraceToGlobalTraceMap(cnt++));
973 
974  if (pFields[j]->GetBndConditions()[b]->
975  GetUserDefined() == "WallViscous" ||
976  pFields[j]->GetBndConditions()[b]->
977  GetUserDefined() == "WallAdiabatic" ||
978  pFields[j]->GetBndConditions()[b]->
979  GetUserDefined() == "Wall")
980  {
981  Vmath::Vcopy(nBndEdgePts, &traceFields[j][id2], 1,
982  &surfaceFields[j][id1], 1);
983 
984  id1 += nBndEdgePts;
985  }
986  }
987  }
988  }
989  }
990 
991  // Extract fields added
992  if (pFields[0]->GetBndCondExpansions().num_elements())
993  {
994  for (j = 0; j < nfieldsAdded; ++j)
995  {
996  cout << "field added " << j << endl;
997 
998  id1 = 0;
999  cnt = 0;
1000  nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
1001  for (b = 0; b < nBndRegions; ++b)
1002  {
1003  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
1004  for (e = 0; e < nBndEdges; ++e)
1005  {
1006  nBndEdgePts = pFields[0]->
1007  GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
1008 
1009  id2 = pFields[0]->GetTrace()->
1010  GetPhys_Offset(pFields[0]->GetTraceMap()->
1011  GetBndCondTraceToGlobalTraceMap(cnt++));
1012 
1013  if (pFields[0]->GetBndConditions()[b]->
1014  GetUserDefined() == "WallViscous" ||
1015  pFields[0]->GetBndConditions()[b]->
1016  GetUserDefined() == "WallAdiabatic" ||
1017  pFields[0]->GetBndConditions()[b]->
1018  GetUserDefined() == "Wall")
1019  {
1020  Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
1021  &surfaceFieldsAdded[j][id1], 1);
1022 
1023  id1 += nBndEdgePts;
1024  }
1025  }
1026  }
1027  }
1028  }
1029 
1030  //==========================================================================
1031  //==========================================================================
1032  //==========================================================================
1033 
1034  // Print the surface coordinates and the surface solution in a .txt file
1035  ofstream outfile;
1036  outfile.open(fname.c_str());
1037  outfile << "% x[m] " << " \t"
1038  << "y[m] " << " \t"
1039  << "z[m] " << " \t"
1040  << "nx[] " << " \t"
1041  << "ny[] " << " \t"
1042  << "nz[] " << " \t"
1043  << "bx[] " << " \t"
1044  << "by[] " << " \t"
1045  << "bz[] " << " \t"
1046  << "tx[] " << " \t"
1047  << "ty[] " << " \t"
1048  << "tz[] " << " \t"
1049  << "rho[kg/m^3] " << " \t"
1050  << "rhou[kg/(m^2 s)] " << " \t"
1051  << "rhov[kg/(m^2 s)] " << " \t"
1052  << "rhow[kg/(m^2 s)] " << " \t"
1053  << "E[Pa] " << " \t"
1054  << "p[Pa] " << " \t"
1055  << "T[k] " << " \t"
1056  << "dT/dn[k/m] " << " \t"
1057  << "dp/dT[Pa/m] " << " \t"
1058  << "dp/dB[Pa/m] " << " \t"
1059  << "dp/dx[Pa/m] " << " \t"
1060  << "dp/dy[Pa/m] " << " \t"
1061  << "dp/dz[Pa/m] " << " \t"
1062  << "du/dx[s^-1] " << " \t"
1063  << "du/dy[s^-1] " << " \t"
1064  << "du/dz[s^-1] " << " \t"
1065  << "dv/dx[s^-1] " << " \t"
1066  << "dv/dy[s^-1] " << " \t"
1067  << "dv/dz[s^-1] " << " \t"
1068  << "dw/dx[s^-1] " << " \t"
1069  << "dw/dy[s^-1] " << " \t"
1070  << "dw/dz[s^-1] " << " \t"
1071  << "tau_xx[Pa] " << " \t"
1072  << "tau_yy[Pa] " << " \t"
1073  << "tau_zz[Pa] " << " \t"
1074  << "tau_xy[Pa] " << " \t"
1075  << "tau_xz[Pa] " << " \t"
1076  << "tau_yz[Pa] " << " \t"
1077  << "mu[Pa s] " << " \t"
1078  << "M[] " << " \t"
1079  << endl;
1080  for (i = 0; i < nSurfacePts; ++i)
1081  {
1082  outfile << scientific
1083  << setw (17)
1084  << setprecision(16)
1085  << surfaceX[i] << " \t "
1086  << surfaceY[i] << " \t "
1087  << surfaceZ[i] << " \t "
1088  << surfaceFieldsAdded[0][i] << " \t "
1089  << surfaceFieldsAdded[1][i] << " \t "
1090  << surfaceFieldsAdded[2][i] << " \t "
1091  << surfaceFieldsAdded[3][i] << " \t "
1092  << surfaceFieldsAdded[4][i] << " \t "
1093  << surfaceFieldsAdded[5][i] << " \t "
1094  << surfaceFieldsAdded[6][i] << " \t "
1095  << surfaceFieldsAdded[7][i] << " \t "
1096  << surfaceFieldsAdded[8][i] << " \t "
1097  << surfaceFields[0][i] << " \t "
1098  << surfaceFields[1][i] << " \t "
1099  << surfaceFields[2][i] << " \t "
1100  << surfaceFields[3][i] << " \t "
1101  << surfaceFields[4][i] << " \t "
1102  << surfaceFieldsAdded[9][i] << " \t "
1103  << surfaceFieldsAdded[10][i] << " \t "
1104  << surfaceFieldsAdded[11][i] << " \t "
1105  << surfaceFieldsAdded[12][i] << " \t "
1106  << surfaceFieldsAdded[13][i] << " \t "
1107  << surfaceFieldsAdded[14][i] << " \t "
1108  << surfaceFieldsAdded[15][i] << " \t "
1109  << surfaceFieldsAdded[16][i] << " \t "
1110  << surfaceFieldsAdded[17][i] << " \t "
1111  << surfaceFieldsAdded[18][i] << " \t "
1112  << surfaceFieldsAdded[19][i] << " \t "
1113  << surfaceFieldsAdded[20][i] << " \t "
1114  << surfaceFieldsAdded[21][i] << " \t "
1115  << surfaceFieldsAdded[22][i] << " \t "
1116  << surfaceFieldsAdded[23][i] << " \t "
1117  << surfaceFieldsAdded[24][i] << " \t "
1118  << surfaceFieldsAdded[25][i] << " \t "
1119  << surfaceFieldsAdded[26][i] << " \t "
1120  << surfaceFieldsAdded[27][i] << " \t "
1121  << surfaceFieldsAdded[28][i] << " \t "
1122  << surfaceFieldsAdded[29][i] << " \t "
1123  << surfaceFieldsAdded[30][i] << " \t "
1124  << surfaceFieldsAdded[31][i] << " \t "
1125  << surfaceFieldsAdded[32][i] << " \t "
1126  << surfaceFieldsAdded[33][i] << " \t "
1127  << endl;
1128  }
1129  outfile << endl << endl;
1130  outfile.close();
1131 
1132  return 0;
1133 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
int main(int argc, char *argv[])
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:279
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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:442
NekDouble m_vInf
NekDouble m_Twall
NekDouble m_uInf
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:241
NekDouble m_mu
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
NekDouble m_rhoInf
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
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:213
double NekDouble
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:112
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:299
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:183