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