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