Nektar++
Functions
ExtractSurface3DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <string>
#include <iostream>
#include <iomanip>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/DisContField.h>
#include <LocalRegions/MatrixKey.h>
#include <LocalRegions/Expansion3D.h>
#include <LocalRegions/Expansion.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Communication/Comm.h>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ContField.h>
#include <SpatialDomains/MeshGraph.h>
#include <SolverUtils/SolverUtilsDeclspec.h>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Evaluation of the velocity gradient in the cartesian directions Du_x: traceFieldsAdded[17] Du_y: traceFieldsAdded[18] Du_z: traceFieldsAdded[19] Dv_x: traceFieldsAdded[20] Dv_y: traceFieldsAdded[21] Dv_z: traceFieldsAdded[22] Dw_x: traceFieldsAdded[23] Dw_y: traceFieldsAdded[24] Dw_z: traceFieldsAdded[25]

Definition at line 65 of file ExtractSurface3DCFS.cpp.

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

References ASSERTL0, Nektar::LibUtilities::ePolyEvenlySpaced, Vmath::Fill(), Nektar::LibUtilities::Import(), m_mu, m_rhoInf, m_Twall, m_uInf, m_vInf, CG_Iterations::pressure, Vmath::Smul(), tinysimd::sqrt(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1