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