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]
70 string fname = std::string(argv[2]);
71 int fdot = fname.find_last_of(
'.');
72 if (fdot != std::string::npos)
74 string ending = fname.substr(fdot);
79 if (ending ==
".chk" || ending ==
".fld")
81 fname = fname.substr(0,fdot);
85 fname = fname +
".txt";
92 int nBndEdgePts, nBndEdges, nBndRegions;
97 "Usage: ExtractSurface3DCFS meshfile fieldFile\n");
99 "Extracts a surface from a 3D fld file" 100 "(only for CompressibleFlowSolver and purely 3D .fld files)\n");
105 = LibUtilities::SessionReader::CreateInstance(3, argv);
107 SpatialDomains::MeshGraph::Read(vSession);
109 std::string m_ViscosityType;
122 int nDimensions = m_spacedim;
126 ASSERTL0(vSession->DefinesParameter(
"Gamma"),
127 "Compressible flow sessions must define a Gamma parameter.");
128 vSession->LoadParameter(
"Gamma", m_gamma, 1.4);
131 ASSERTL0(vSession->DefinesParameter(
"pInf"),
132 "Compressible flow sessions must define a pInf parameter.");
133 vSession->LoadParameter(
"pInf", m_pInf, 101325);
136 ASSERTL0(vSession->DefinesParameter(
"rhoInf"),
137 "Compressible flow sessions must define a rhoInf parameter.");
138 vSession->LoadParameter(
"rhoInf", m_rhoInf, 1.225);
141 ASSERTL0(vSession->DefinesParameter(
"uInf"),
142 "Compressible flow sessions must define a uInf parameter.");
143 vSession->LoadParameter(
"uInf", m_uInf, 0.1);
146 if (m_spacedim == 2 || m_spacedim == 3)
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);
157 ASSERTL0(vSession->DefinesParameter(
"wInf"),
158 "Compressible flow sessions must define a wInf parameter" 160 vSession->LoadParameter(
"wInf", m_wInf, 0.0);
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);
170 string fieldFile(argv[2]);
171 vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
172 vector<vector<NekDouble> > fieldData;
179 vector< vector<LibUtilities::PointsType> > pointsType;
180 for (i = 0; i < fieldDef.size(); ++i)
182 vector<LibUtilities::PointsType> ptype;
183 for (j = 0; j < 3; ++j)
187 pointsType.push_back(ptype);
189 graphShPt->SetExpansions(fieldDef, pointsType);
196 int nfields = fieldDef[0]->m_fields.size();
200 for(i = 0; i < pFields.num_elements(); i++)
203 ::DisContField3D>::AllocateSharedPtr(vSession, graphShPt,
204 vSession->GetVariable(i));
213 for (i = 1; i < nfields; ++i)
221 if (pFields[0]->GetBndCondExpansions().num_elements())
225 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
226 for (b = 0; b < nBndRegions; ++b)
228 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
229 for (e = 0; e < nBndEdges; ++e)
231 nBndEdgePts = pFields[0]->
232 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
234 if (pFields[0]->GetBndConditions()[b]->
235 GetUserDefined() ==
"WallViscous" ||
236 pFields[0]->GetBndConditions()[b]->
237 GetUserDefined() ==
"WallAdiabatic" ||
238 pFields[0]->GetBndConditions()[b]->
239 GetUserDefined() ==
"Wall")
241 nSurfacePts += nBndEdgePts;
248 int nSolutionPts = pFields[0]->GetNpoints();
249 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
250 int nElements = pFields[0]->GetExpSize();
266 pFields[0]->GetCoords(x, y, z);
268 pFields[0]->ExtractTracePhys(x, traceX);
269 pFields[0]->ExtractTracePhys(y, traceY);
270 pFields[0]->ExtractTracePhys(z, traceZ);
280 for (j = 0; j < nfields; ++j)
287 for (i = 0; i < fieldData.size(); ++i)
289 Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
290 fieldDef[i]->m_fields[j],
291 Exp[j]->UpdateCoeffs());
293 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
294 Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
295 pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
300 int nfieldsAdded = 34;
304 for (j = 0; j < nfieldsAdded; ++j)
324 for(i = 0; i < nDimensions; ++i)
328 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
338 for(i = 0; i < nDimensions; ++i)
349 &m_traceNormals[0][0], 1,
350 &traceFieldsAdded[0][0], 1);
354 &m_traceNormals[1][0], 1,
355 &traceFieldsAdded[1][0], 1);
359 &m_traceNormals[2][0], 1,
360 &traceFieldsAdded[2][0], 1);
366 &m_traceNormals[0][0], 1,
371 &m_traceNormals[1][0], 1,
375 &m_traceNormals[2][0], 1,
379 for (i = 0; i < m_spacedim; i++)
382 &NormH[0],1, &NormH[0],1);
398 &m_traceBinormals[0][0], 1);
401 &m_traceBinormals[0][0], 1,
402 &traceFieldsAdded[3][0], 1);
423 &m_traceBinormals[1][0], 1);
426 &m_traceBinormals[1][0], 1,
427 &traceFieldsAdded[4][0], 1);
443 &m_traceBinormals[2][0], 1);
446 &m_traceBinormals[2][0], 1,
447 &traceFieldsAdded[5][0], 1);
463 &m_traceTangents[0][0], 1);
466 &m_traceTangents[0][0], 1,
467 &traceFieldsAdded[6][0], 1);
471 &m_traceBinormals[2][0], 1,
472 &m_traceTangents[1][0], 1);
475 &m_traceTangents[1][0], 1,
476 &traceFieldsAdded[7][0], 1);
497 &m_traceTangents[2][0], 1);
500 &m_traceTangents[2][0], 1,
501 &traceFieldsAdded[8][0], 1);
513 for (i = 0; i < m_spacedim; i++)
516 &uFields[i + 1][0], 1,
517 &uFields[i + 1][0], 1,
537 &uFields[nfields - 1][0], 1,
546 pFields[0]->ExtractTracePhys(
pressure, traceFieldsAdded[9]);
560 NekDouble GasConstantInv = 1.0/m_gasConstant;
566 pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[10]);
575 for (i = 0; i < nDimensions; ++ i)
581 for (i = 0; i < nDimensions; ++ i)
583 for (n = 0; n < nElements; n++)
585 phys_offset = pFields[0]->GetPhys_Offset(n);
587 pFields[i]->GetExp(n)->PhysDeriv(
588 i, temperature + phys_offset,
589 auxArray = Dtemperature[i] + phys_offset);
592 pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
595 for(i = 0; i < nDimensions; ++i)
598 &m_traceNormals[i][0], 1,
599 &traceDtemperature[i][0], 1,
603 &traceFieldsAdded[11][0], 1,
605 &traceFieldsAdded[11][0], 1);
619 for (i = 0; i < nDimensions; ++ i)
625 for (i = 0; i < nDimensions; ++ i)
627 for (n = 0; n < nElements; n++)
629 phys_offset = pFields[0]->GetPhys_Offset(n);
631 pFields[i]->GetExp(n)->PhysDeriv(
633 auxArray = Dpressure[i] + phys_offset);
636 pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
640 for(i = 0; i < nDimensions; ++i)
643 &m_traceTangents[i][0], 1,
644 &traceDpressure[i][0], 1,
648 &traceFieldsAdded[12][0], 1,
650 &traceFieldsAdded[12][0], 1);
654 for(i = 0; i < nDimensions; ++i)
657 &m_traceBinormals[i][0], 1,
658 &traceDpressure[i][0], 1,
662 &traceFieldsAdded[13][0], 1,
664 &traceFieldsAdded[13][0], 1);
670 &traceDpressure[0][0], 1,
671 &traceFieldsAdded[14][0], 1);
675 &traceDpressure[1][0], 1,
676 &traceFieldsAdded[15][0], 1);
680 &traceDpressure[2][0], 1,
681 &traceFieldsAdded[16][0], 1);
700 for (i = 0; i < nDimensions; ++ i)
706 Vmath::Vdiv(nSolutionPts, uFields[i+1], 1, uFields[0], 1,
709 for (j = 0; j < nDimensions; ++j)
716 for (i = 0; i < nDimensions; ++i)
718 for (j = 0; j < nDimensions; ++j)
720 for (n = 0; n < nElements; n++)
722 phys_offset = pFields[0]->GetPhys_Offset(n);
724 pFields[i]->GetExp(n)->PhysDeriv(
725 j, velocity[i] + phys_offset,
726 auxArray = Dvelocity[i][j] + phys_offset);
730 pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
735 &traceDvelocity[0][0][0], 1,
736 &traceFieldsAdded[17][0], 1);
738 &traceDvelocity[0][1][0], 1,
739 &traceFieldsAdded[18][0], 1);
741 &traceDvelocity[0][2][0], 1,
742 &traceFieldsAdded[19][0], 1);
744 &traceDvelocity[1][0][0], 1,
745 &traceFieldsAdded[20][0], 1);
747 &traceDvelocity[1][1][0], 1,
748 &traceFieldsAdded[21][0], 1);
750 &traceDvelocity[1][2][0], 1,
751 &traceFieldsAdded[22][0], 1);
753 &traceDvelocity[2][0][0], 1,
754 &traceFieldsAdded[23][0], 1);
756 &traceDvelocity[2][1][0], 1,
757 &traceFieldsAdded[24][0], 1);
759 &traceDvelocity[2][2][0], 1,
760 &traceFieldsAdded[25][0], 1);
781 if (m_ViscosityType ==
"Variable")
784 NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
787 for (
int i = 0; i < nSolutionPts; ++i)
789 ratio = temperature[i] / T_star;
790 mu[i] = mu_star * ratio * sqrt(ratio) *
791 (T_star + 110.0) / (temperature[i] + 110.0);
804 Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
808 &Dvelocity[0][0][0], 1, &divVel[0], 1);
810 &Dvelocity[1][1][0], 1, &divVel[0], 1);
813 Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
814 Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
818 for (j = 0; j < m_spacedim; ++j)
823 Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
826 Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
836 &Dvelocity[1][0][0], 1, &Sxy[0], 1);
840 &Dvelocity[2][0][0], 1, &Sxz[0], 1);
844 &Dvelocity[2][1][0], 1, &Syz[0], 1);
847 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
850 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
853 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
857 pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[26]);
858 pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[27]);
859 pFields[0]->ExtractTracePhys(Sgg[2], traceFieldsAdded[28]);
860 pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[29]);
861 pFields[0]->ExtractTracePhys(Sxz, traceFieldsAdded[30]);
862 pFields[0]->ExtractTracePhys(Syz, traceFieldsAdded[31]);
868 pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[32]);
879 Vmath::Smul (nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
880 Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
885 for (
int i = 0; i < m_spacedim; ++i)
893 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
894 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
896 Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
898 pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[33]);
905 if (pFields[0]->GetBndCondExpansions().num_elements())
909 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
910 for (b = 0; b < nBndRegions; ++b)
912 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
913 for (e = 0; e < nBndEdges; ++e)
915 nBndEdgePts = pFields[0]->
916 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
918 id2 = pFields[0]->GetTrace()->
919 GetPhys_Offset(pFields[0]->GetTraceMap()->
920 GetBndCondTraceToGlobalTraceMap(cnt++));
922 if (pFields[0]->GetBndConditions()[b]->
923 GetUserDefined() ==
"WallViscous" ||
924 pFields[0]->GetBndConditions()[b]->
925 GetUserDefined() ==
"WallAdiabatic" ||
926 pFields[0]->GetBndConditions()[b]->
927 GetUserDefined() ==
"Wall")
946 if (pFields[0]->GetBndCondExpansions().num_elements())
949 for (j = 0; j < nfields; ++j)
951 cout <<
"field " << j << endl;
955 nBndRegions = pFields[j]->GetBndCondExpansions().num_elements();
956 for (b = 0; b < nBndRegions; ++b)
958 nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
959 for (e = 0; e < nBndEdges; ++e)
961 nBndEdgePts = pFields[j]->
962 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
964 id2 = pFields[j]->GetTrace()->
965 GetPhys_Offset(pFields[j]->GetTraceMap()->
966 GetBndCondTraceToGlobalTraceMap(cnt++));
968 if (pFields[j]->GetBndConditions()[b]->
969 GetUserDefined() ==
"WallViscous" ||
970 pFields[j]->GetBndConditions()[b]->
971 GetUserDefined() ==
"WallAdiabatic" ||
972 pFields[j]->GetBndConditions()[b]->
973 GetUserDefined() ==
"Wall")
976 &surfaceFields[j][id1], 1);
986 if (pFields[0]->GetBndCondExpansions().num_elements())
988 for (j = 0; j < nfieldsAdded; ++j)
990 cout <<
"field added " << j << endl;
994 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
995 for (b = 0; b < nBndRegions; ++b)
997 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
998 for (e = 0; e < nBndEdges; ++e)
1000 nBndEdgePts = pFields[0]->
1001 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
1003 id2 = pFields[0]->GetTrace()->
1004 GetPhys_Offset(pFields[0]->GetTraceMap()->
1005 GetBndCondTraceToGlobalTraceMap(cnt++));
1007 if (pFields[0]->GetBndConditions()[b]->
1008 GetUserDefined() ==
"WallViscous" ||
1009 pFields[0]->GetBndConditions()[b]->
1010 GetUserDefined() ==
"WallAdiabatic" ||
1011 pFields[0]->GetBndConditions()[b]->
1012 GetUserDefined() ==
"Wall")
1014 Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
1015 &surfaceFieldsAdded[j][id1], 1);
1030 outfile.open(fname.c_str());
1031 outfile <<
"% x[m] " <<
" \t" 1043 <<
"rho[kg/m^3] " <<
" \t" 1044 <<
"rhou[kg/(m^2 s)] " <<
" \t" 1045 <<
"rhov[kg/(m^2 s)] " <<
" \t" 1046 <<
"rhow[kg/(m^2 s)] " <<
" \t" 1047 <<
"E[Pa] " <<
" \t" 1048 <<
"p[Pa] " <<
" \t" 1050 <<
"dT/dn[k/m] " <<
" \t" 1051 <<
"dp/dT[Pa/m] " <<
" \t" 1052 <<
"dp/dB[Pa/m] " <<
" \t" 1053 <<
"dp/dx[Pa/m] " <<
" \t" 1054 <<
"dp/dy[Pa/m] " <<
" \t" 1055 <<
"dp/dz[Pa/m] " <<
" \t" 1056 <<
"du/dx[s^-1] " <<
" \t" 1057 <<
"du/dy[s^-1] " <<
" \t" 1058 <<
"du/dz[s^-1] " <<
" \t" 1059 <<
"dv/dx[s^-1] " <<
" \t" 1060 <<
"dv/dy[s^-1] " <<
" \t" 1061 <<
"dv/dz[s^-1] " <<
" \t" 1062 <<
"dw/dx[s^-1] " <<
" \t" 1063 <<
"dw/dy[s^-1] " <<
" \t" 1064 <<
"dw/dz[s^-1] " <<
" \t" 1065 <<
"tau_xx[Pa] " <<
" \t" 1066 <<
"tau_yy[Pa] " <<
" \t" 1067 <<
"tau_zz[Pa] " <<
" \t" 1068 <<
"tau_xy[Pa] " <<
" \t" 1069 <<
"tau_xz[Pa] " <<
" \t" 1070 <<
"tau_yz[Pa] " <<
" \t" 1071 <<
"mu[Pa s] " <<
" \t" 1074 for (i = 0; i < nSurfacePts; ++i)
1076 outfile << scientific
1079 << surfaceX[i] <<
" \t " 1080 << surfaceY[i] <<
" \t " 1081 << surfaceZ[i] <<
" \t " 1082 << surfaceFieldsAdded[0][i] <<
" \t " 1083 << surfaceFieldsAdded[1][i] <<
" \t " 1084 << surfaceFieldsAdded[2][i] <<
" \t " 1085 << surfaceFieldsAdded[3][i] <<
" \t " 1086 << surfaceFieldsAdded[4][i] <<
" \t " 1087 << surfaceFieldsAdded[5][i] <<
" \t " 1088 << surfaceFieldsAdded[6][i] <<
" \t " 1089 << surfaceFieldsAdded[7][i] <<
" \t " 1090 << surfaceFieldsAdded[8][i] <<
" \t " 1091 << surfaceFields[0][i] <<
" \t " 1092 << surfaceFields[1][i] <<
" \t " 1093 << surfaceFields[2][i] <<
" \t " 1094 << surfaceFields[3][i] <<
" \t " 1095 << surfaceFields[4][i] <<
" \t " 1096 << surfaceFieldsAdded[9][i] <<
" \t " 1097 << surfaceFieldsAdded[10][i] <<
" \t " 1098 << surfaceFieldsAdded[11][i] <<
" \t " 1099 << surfaceFieldsAdded[12][i] <<
" \t " 1100 << surfaceFieldsAdded[13][i] <<
" \t " 1101 << surfaceFieldsAdded[14][i] <<
" \t " 1102 << surfaceFieldsAdded[15][i] <<
" \t " 1103 << surfaceFieldsAdded[16][i] <<
" \t " 1104 << surfaceFieldsAdded[17][i] <<
" \t " 1105 << surfaceFieldsAdded[18][i] <<
" \t " 1106 << surfaceFieldsAdded[19][i] <<
" \t " 1107 << surfaceFieldsAdded[20][i] <<
" \t " 1108 << surfaceFieldsAdded[21][i] <<
" \t " 1109 << surfaceFieldsAdded[22][i] <<
" \t " 1110 << surfaceFieldsAdded[23][i] <<
" \t " 1111 << surfaceFieldsAdded[24][i] <<
" \t " 1112 << surfaceFieldsAdded[25][i] <<
" \t " 1113 << surfaceFieldsAdded[26][i] <<
" \t " 1114 << surfaceFieldsAdded[27][i] <<
" \t " 1115 << surfaceFieldsAdded[28][i] <<
" \t " 1116 << surfaceFieldsAdded[29][i] <<
" \t " 1117 << surfaceFieldsAdded[30][i] <<
" \t " 1118 << surfaceFieldsAdded[31][i] <<
" \t " 1119 << surfaceFieldsAdded[32][i] <<
" \t " 1120 << surfaceFieldsAdded[33][i] <<
" \t " 1123 outfile << endl << endl;
#define ASSERTL0(condition, msg)
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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...
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.
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
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.
1D Evenly-spaced points using Lagrange polynomial
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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.
std::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
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.
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.