70 int main(
int argc,
char *argv[])
72 string fname = std::string(argv[2]);
73 int fdot = fname.find_last_of(
'.');
74 if (fdot != std::string::npos)
76 string ending = fname.substr(fdot);
81 if (ending ==
".chk" || ending ==
".fld")
83 fname = fname.substr(0,fdot);
87 fname = fname +
".txt";
94 int nBndEdgePts, nBndEdges, nBndRegions;
99 "Usage: ExtractSurface3DCFS meshfile fieldFile\n");
101 "Extracts a surface from a 3D fld file"
102 "(only for CompressibleFlowSolver and purely 3D .fld files)\n");
109 std::string m_ViscosityType;
123 int nDimensions = m_spacedim;
127 ASSERTL0(vSession->DefinesParameter(
"Gamma"),
128 "Compressible flow sessions must define a Gamma parameter.");
129 vSession->LoadParameter(
"Gamma", m_gamma, 1.4);
132 ASSERTL0(vSession->DefinesParameter(
"pInf"),
133 "Compressible flow sessions must define a pInf parameter.");
134 vSession->LoadParameter(
"pInf", m_pInf, 101325);
137 ASSERTL0(vSession->DefinesParameter(
"rhoInf"),
138 "Compressible flow sessions must define a rhoInf parameter.");
139 vSession->LoadParameter(
"rhoInf", m_rhoInf, 1.225);
142 ASSERTL0(vSession->DefinesParameter(
"uInf"),
143 "Compressible flow sessions must define a uInf parameter.");
144 vSession->LoadParameter(
"uInf", m_uInf, 0.1);
147 if (m_spacedim == 2 || m_spacedim == 3)
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);
158 ASSERTL0(vSession->DefinesParameter(
"wInf"),
159 "Compressible flow sessions must define a wInf parameter"
161 vSession->LoadParameter(
"wInf", m_wInf, 0.0);
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);
173 string meshfile(argv[1]);
180 string fieldFile(argv[2]);
181 vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
182 vector<vector<NekDouble> > fieldData;
189 vector< vector<LibUtilities::PointsType> > pointsType;
190 for (i = 0; i < fieldDef.size(); ++i)
192 vector<LibUtilities::PointsType> ptype;
193 for (j = 0; j < 3; ++j)
197 pointsType.push_back(ptype);
199 graphShPt->SetExpansions(fieldDef, pointsType);
206 int nfields = fieldDef[0]->m_fields.size();
210 for(i = 0; i < pFields.num_elements(); i++)
214 vSession->GetVariable(i));
223 for (i = 1; i < nfields; ++i)
231 if (pFields[0]->GetBndCondExpansions().num_elements())
235 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
236 for (b = 0; b < nBndRegions; ++b)
238 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
239 for (e = 0; e < nBndEdges; ++e)
241 nBndEdgePts = pFields[0]->
242 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
244 if (pFields[0]->GetBndConditions()[b]->
245 GetUserDefined() ==
"WallViscous" ||
246 pFields[0]->GetBndConditions()[b]->
247 GetUserDefined() ==
"WallAdiabatic" ||
248 pFields[0]->GetBndConditions()[b]->
249 GetUserDefined() ==
"Wall")
251 nSurfacePts += nBndEdgePts;
258 int nSolutionPts = pFields[0]->GetNpoints();
259 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
260 int nElements = pFields[0]->GetExpSize();
276 pFields[0]->GetCoords(x, y, z);
278 pFields[0]->ExtractTracePhys(x, traceX);
279 pFields[0]->ExtractTracePhys(y, traceY);
280 pFields[0]->ExtractTracePhys(z, traceZ);
290 for (j = 0; j < nfields; ++j)
297 for (i = 0; i < fieldData.size(); ++i)
299 Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
300 fieldDef[i]->m_fields[j],
301 Exp[j]->UpdateCoeffs());
303 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
304 Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
305 pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
310 int nfieldsAdded = 34;
314 for (j = 0; j < nfieldsAdded; ++j)
334 for(i = 0; i < nDimensions; ++i)
338 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
348 for(i = 0; i < nDimensions; ++i)
359 &m_traceNormals[0][0], 1,
360 &traceFieldsAdded[0][0], 1);
364 &m_traceNormals[1][0], 1,
365 &traceFieldsAdded[1][0], 1);
369 &m_traceNormals[2][0], 1,
370 &traceFieldsAdded[2][0], 1);
376 &m_traceNormals[0][0], 1,
381 &m_traceNormals[1][0], 1,
385 &m_traceNormals[2][0], 1,
389 for (i = 0; i < m_spacedim; i++)
392 &NormH[0],1, &NormH[0],1);
408 &m_traceBinormals[0][0], 1);
411 &m_traceBinormals[0][0], 1,
412 &traceFieldsAdded[3][0], 1);
433 &m_traceBinormals[1][0], 1);
436 &m_traceBinormals[1][0], 1,
437 &traceFieldsAdded[4][0], 1);
453 &m_traceBinormals[2][0], 1);
456 &m_traceBinormals[2][0], 1,
457 &traceFieldsAdded[5][0], 1);
473 &m_traceTangents[0][0], 1);
476 &m_traceTangents[0][0], 1,
477 &traceFieldsAdded[6][0], 1);
481 &m_traceBinormals[2][0], 1,
482 &m_traceTangents[1][0], 1);
485 &m_traceTangents[1][0], 1,
486 &traceFieldsAdded[7][0], 1);
507 &m_traceTangents[2][0], 1);
510 &m_traceTangents[2][0], 1,
511 &traceFieldsAdded[8][0], 1);
523 for (i = 0; i < m_spacedim; i++)
526 &uFields[i + 1][0], 1,
527 &uFields[i + 1][0], 1,
547 &uFields[nfields - 1][0], 1,
556 pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[9]);
570 NekDouble GasConstantInv = 1.0/m_gasConstant;
576 pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[10]);
585 for (i = 0; i < nDimensions; ++ i)
591 for (i = 0; i < nDimensions; ++ i)
593 for (n = 0; n < nElements; n++)
595 phys_offset = pFields[0]->GetPhys_Offset(n);
597 pFields[i]->GetExp(n)->PhysDeriv(
598 i, temperature + phys_offset,
599 auxArray = Dtemperature[i] + phys_offset);
602 pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
605 for(i = 0; i < nDimensions; ++i)
608 &m_traceNormals[i][0], 1,
609 &traceDtemperature[i][0], 1,
613 &traceFieldsAdded[11][0], 1,
615 &traceFieldsAdded[11][0], 1);
629 for (i = 0; i < nDimensions; ++ i)
635 for (i = 0; i < nDimensions; ++ i)
637 for (n = 0; n < nElements; n++)
639 phys_offset = pFields[0]->GetPhys_Offset(n);
641 pFields[i]->GetExp(n)->PhysDeriv(
642 i, pressure + phys_offset,
643 auxArray = Dpressure[i] + phys_offset);
646 pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
650 for(i = 0; i < nDimensions; ++i)
653 &m_traceTangents[i][0], 1,
654 &traceDpressure[i][0], 1,
658 &traceFieldsAdded[12][0], 1,
660 &traceFieldsAdded[12][0], 1);
664 for(i = 0; i < nDimensions; ++i)
667 &m_traceBinormals[i][0], 1,
668 &traceDpressure[i][0], 1,
672 &traceFieldsAdded[13][0], 1,
674 &traceFieldsAdded[13][0], 1);
680 &traceDpressure[0][0], 1,
681 &traceFieldsAdded[14][0], 1);
685 &traceDpressure[1][0], 1,
686 &traceFieldsAdded[15][0], 1);
690 &traceDpressure[2][0], 1,
691 &traceFieldsAdded[16][0], 1);
710 for (i = 0; i < nDimensions; ++ i)
716 Vmath::Vdiv(nSolutionPts, uFields[i+1], 1, uFields[0], 1,
719 for (j = 0; j < nDimensions; ++j)
726 for (i = 0; i < nDimensions; ++i)
728 for (j = 0; j < nDimensions; ++j)
730 for (n = 0; n < nElements; n++)
732 phys_offset = pFields[0]->GetPhys_Offset(n);
734 pFields[i]->GetExp(n)->PhysDeriv(
735 j, velocity[i] + phys_offset,
736 auxArray = Dvelocity[i][j] + phys_offset);
740 pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
745 &traceDvelocity[0][0][0], 1,
746 &traceFieldsAdded[17][0], 1);
748 &traceDvelocity[0][1][0], 1,
749 &traceFieldsAdded[18][0], 1);
751 &traceDvelocity[0][2][0], 1,
752 &traceFieldsAdded[19][0], 1);
754 &traceDvelocity[1][0][0], 1,
755 &traceFieldsAdded[20][0], 1);
757 &traceDvelocity[1][1][0], 1,
758 &traceFieldsAdded[21][0], 1);
760 &traceDvelocity[1][2][0], 1,
761 &traceFieldsAdded[22][0], 1);
763 &traceDvelocity[2][0][0], 1,
764 &traceFieldsAdded[23][0], 1);
766 &traceDvelocity[2][1][0], 1,
767 &traceFieldsAdded[24][0], 1);
769 &traceDvelocity[2][2][0], 1,
770 &traceFieldsAdded[25][0], 1);
791 if (m_ViscosityType ==
"Variable")
794 NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
797 for (
int i = 0; i < nSolutionPts; ++i)
799 ratio = temperature[i] / T_star;
800 mu[i] = mu_star * ratio * sqrt(ratio) *
801 (T_star + 110.0) / (temperature[i] + 110.0);
814 Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
818 &Dvelocity[0][0][0], 1, &divVel[0], 1);
820 &Dvelocity[1][1][0], 1, &divVel[0], 1);
823 Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
824 Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
828 for (j = 0; j < m_spacedim; ++j)
833 Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
836 Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
846 &Dvelocity[1][0][0], 1, &Sxy[0], 1);
850 &Dvelocity[2][0][0], 1, &Sxz[0], 1);
854 &Dvelocity[2][1][0], 1, &Syz[0], 1);
857 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
860 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
863 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
867 pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[26]);
868 pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[27]);
869 pFields[0]->ExtractTracePhys(Sgg[2], traceFieldsAdded[28]);
870 pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[29]);
871 pFields[0]->ExtractTracePhys(Sxz, traceFieldsAdded[30]);
872 pFields[0]->ExtractTracePhys(Syz, traceFieldsAdded[31]);
878 pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[32]);
888 Vmath::Vdiv (nSolutionPts, pressure, 1, uFields[0], 1, soundspeed, 1);
889 Vmath::Smul (nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
890 Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
895 for (
int i = 0; i < m_spacedim; ++i)
903 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
904 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
906 Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
908 pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[33]);
915 if (pFields[0]->GetBndCondExpansions().num_elements())
919 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
920 for (b = 0; b < nBndRegions; ++b)
922 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
923 for (e = 0; e < nBndEdges; ++e)
925 nBndEdgePts = pFields[0]->
926 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
928 id2 = pFields[0]->GetTrace()->
929 GetPhys_Offset(pFields[0]->GetTraceMap()->
930 GetBndCondTraceToGlobalTraceMap(cnt++));
932 if (pFields[0]->GetBndConditions()[b]->
933 GetUserDefined() ==
"WallViscous" ||
934 pFields[0]->GetBndConditions()[b]->
935 GetUserDefined() ==
"WallAdiabatic" ||
936 pFields[0]->GetBndConditions()[b]->
937 GetUserDefined() ==
"Wall")
956 if (pFields[0]->GetBndCondExpansions().num_elements())
959 for (j = 0; j < nfields; ++j)
961 cout <<
"field " << j << endl;
965 nBndRegions = pFields[j]->GetBndCondExpansions().num_elements();
966 for (b = 0; b < nBndRegions; ++b)
968 nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
969 for (e = 0; e < nBndEdges; ++e)
971 nBndEdgePts = pFields[j]->
972 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
974 id2 = pFields[j]->GetTrace()->
975 GetPhys_Offset(pFields[j]->GetTraceMap()->
976 GetBndCondTraceToGlobalTraceMap(cnt++));
978 if (pFields[j]->GetBndConditions()[b]->
979 GetUserDefined() ==
"WallViscous" ||
980 pFields[j]->GetBndConditions()[b]->
981 GetUserDefined() ==
"WallAdiabatic" ||
982 pFields[j]->GetBndConditions()[b]->
983 GetUserDefined() ==
"Wall")
986 &surfaceFields[j][id1], 1);
996 if (pFields[0]->GetBndCondExpansions().num_elements())
998 for (j = 0; j < nfieldsAdded; ++j)
1000 cout <<
"field added " << j << endl;
1004 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
1005 for (b = 0; b < nBndRegions; ++b)
1007 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
1008 for (e = 0; e < nBndEdges; ++e)
1010 nBndEdgePts = pFields[0]->
1011 GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();
1013 id2 = pFields[0]->GetTrace()->
1014 GetPhys_Offset(pFields[0]->GetTraceMap()->
1015 GetBndCondTraceToGlobalTraceMap(cnt++));
1017 if (pFields[0]->GetBndConditions()[b]->
1018 GetUserDefined() ==
"WallViscous" ||
1019 pFields[0]->GetBndConditions()[b]->
1020 GetUserDefined() ==
"WallAdiabatic" ||
1021 pFields[0]->GetBndConditions()[b]->
1022 GetUserDefined() ==
"Wall")
1024 Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
1025 &surfaceFieldsAdded[j][id1], 1);
1040 outfile.open(fname.c_str());
1041 outfile <<
"% x[m] " <<
" \t"
1053 <<
"rho[kg/m^3] " <<
" \t"
1054 <<
"rhou[kg/(m^2 s)] " <<
" \t"
1055 <<
"rhov[kg/(m^2 s)] " <<
" \t"
1056 <<
"rhow[kg/(m^2 s)] " <<
" \t"
1057 <<
"E[Pa] " <<
" \t"
1058 <<
"p[Pa] " <<
" \t"
1060 <<
"dT/dn[k/m] " <<
" \t"
1061 <<
"dp/dT[Pa/m] " <<
" \t"
1062 <<
"dp/dB[Pa/m] " <<
" \t"
1063 <<
"dp/dx[Pa/m] " <<
" \t"
1064 <<
"dp/dy[Pa/m] " <<
" \t"
1065 <<
"dp/dz[Pa/m] " <<
" \t"
1066 <<
"du/dx[s^-1] " <<
" \t"
1067 <<
"du/dy[s^-1] " <<
" \t"
1068 <<
"du/dz[s^-1] " <<
" \t"
1069 <<
"dv/dx[s^-1] " <<
" \t"
1070 <<
"dv/dy[s^-1] " <<
" \t"
1071 <<
"dv/dz[s^-1] " <<
" \t"
1072 <<
"dw/dx[s^-1] " <<
" \t"
1073 <<
"dw/dy[s^-1] " <<
" \t"
1074 <<
"dw/dz[s^-1] " <<
" \t"
1075 <<
"tau_xx[Pa] " <<
" \t"
1076 <<
"tau_yy[Pa] " <<
" \t"
1077 <<
"tau_zz[Pa] " <<
" \t"
1078 <<
"tau_xy[Pa] " <<
" \t"
1079 <<
"tau_xz[Pa] " <<
" \t"
1080 <<
"tau_yz[Pa] " <<
" \t"
1081 <<
"mu[Pa s] " <<
" \t"
1084 for (i = 0; i < nSurfacePts; ++i)
1086 outfile << scientific
1089 << surfaceX[i] <<
" \t "
1090 << surfaceY[i] <<
" \t "
1091 << surfaceZ[i] <<
" \t "
1092 << surfaceFieldsAdded[0][i] <<
" \t "
1093 << surfaceFieldsAdded[1][i] <<
" \t "
1094 << surfaceFieldsAdded[2][i] <<
" \t "
1095 << surfaceFieldsAdded[3][i] <<
" \t "
1096 << surfaceFieldsAdded[4][i] <<
" \t "
1097 << surfaceFieldsAdded[5][i] <<
" \t "
1098 << surfaceFieldsAdded[6][i] <<
" \t "
1099 << surfaceFieldsAdded[7][i] <<
" \t "
1100 << surfaceFieldsAdded[8][i] <<
" \t "
1101 << surfaceFields[0][i] <<
" \t "
1102 << surfaceFields[1][i] <<
" \t "
1103 << surfaceFields[2][i] <<
" \t "
1104 << surfaceFields[3][i] <<
" \t "
1105 << surfaceFields[4][i] <<
" \t "
1106 << surfaceFieldsAdded[9][i] <<
" \t "
1107 << surfaceFieldsAdded[10][i] <<
" \t "
1108 << surfaceFieldsAdded[11][i] <<
" \t "
1109 << surfaceFieldsAdded[12][i] <<
" \t "
1110 << surfaceFieldsAdded[13][i] <<
" \t "
1111 << surfaceFieldsAdded[14][i] <<
" \t "
1112 << surfaceFieldsAdded[15][i] <<
" \t "
1113 << surfaceFieldsAdded[16][i] <<
" \t "
1114 << surfaceFieldsAdded[17][i] <<
" \t "
1115 << surfaceFieldsAdded[18][i] <<
" \t "
1116 << surfaceFieldsAdded[19][i] <<
" \t "
1117 << surfaceFieldsAdded[20][i] <<
" \t "
1118 << surfaceFieldsAdded[21][i] <<
" \t "
1119 << surfaceFieldsAdded[22][i] <<
" \t "
1120 << surfaceFieldsAdded[23][i] <<
" \t "
1121 << surfaceFieldsAdded[24][i] <<
" \t "
1122 << surfaceFieldsAdded[25][i] <<
" \t "
1123 << surfaceFieldsAdded[26][i] <<
" \t "
1124 << surfaceFieldsAdded[27][i] <<
" \t "
1125 << surfaceFieldsAdded[28][i] <<
" \t "
1126 << surfaceFieldsAdded[29][i] <<
" \t "
1127 << surfaceFieldsAdded[30][i] <<
" \t "
1128 << surfaceFieldsAdded[31][i] <<
" \t "
1129 << surfaceFieldsAdded[32][i] <<
" \t "
1130 << surfaceFieldsAdded[33][i] <<
" \t "
1133 outfile << endl << endl;
#define ASSERTL0(condition, msg)
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
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.
1D Evenly-spaced points using Lagrange polynomial
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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.
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.