63#include <boost/program_options.hpp>
68namespace po = boost::program_options;
70int main(
int argc,
char *argv[])
73 po::options_description desc(
"Available options");
74 desc.add_options()(
"help,h",
"Produce this help message.");
76 po::options_description hidden(
"Hidden options");
77 hidden.add_options()(
"input-file", po::value<vector<string>>(),
80 po::options_description cmdline_options;
81 cmdline_options.add(desc).add(hidden);
83 po::options_description visible(
"Allowed options");
86 po::positional_options_description
p;
87 p.add(
"input-file", -1);
93 po::store(po::command_line_parser(argc, argv)
94 .options(cmdline_options)
100 catch (
const std::exception &e)
102 cerr << e.what() << endl;
107 std::vector<std::string> filenames;
108 if (vm.count(
"input-file"))
110 filenames = vm[
"input-file"].as<std::vector<std::string>>();
113 if (vm.count(
"help") || vm.count(
"input-file") != 1)
115 cerr <<
"Description: Extracts a surface from a 2D .fld file created "
116 "by the CompressibleFlowSolver and places it into a .cfs file"
118 cerr <<
"Usage: ExtractSurface2DCFS [options] meshFile fieldFile"
125 LibUtilities::SessionReader::CreateInstance(argc, argv);
127 SpatialDomains::MeshGraph::Read(vSession);
129 fname = vSession->GetSessionName() +
".cfs";
132 std::string fieldFile;
133 for (
auto &file : filenames)
135 if (file.size() > 4 && (file.substr(file.size() - 4, 4) ==
".fld" ||
136 file.substr(file.size() - 4, 4) ==
".chk"))
149 int nBndEdgePts, nBndEdges, nBndRegions;
151 std::string m_ViscosityType;
163 int nDimensions = m_spacedim;
167 ASSERTL0(vSession->DefinesParameter(
"Gamma"),
168 "Compressible flow sessions must define a Gamma parameter.");
169 vSession->LoadParameter(
"Gamma", m_gamma, 1.4);
172 ASSERTL0(vSession->DefinesParameter(
"pInf"),
173 "Compressible flow sessions must define a pInf parameter.");
174 vSession->LoadParameter(
"pInf", m_pInf, 101325);
177 ASSERTL0(vSession->DefinesParameter(
"rhoInf"),
178 "Compressible flow sessions must define a rhoInf parameter.");
179 vSession->LoadParameter(
"rhoInf",
m_rhoInf, 1.225);
182 ASSERTL0(vSession->DefinesParameter(
"uInf"),
183 "Compressible flow sessions must define a uInf parameter.");
184 vSession->LoadParameter(
"uInf",
m_uInf, 0.1);
187 if (m_spacedim == 2 || m_spacedim == 3)
189 ASSERTL0(vSession->DefinesParameter(
"vInf"),
190 "Compressible flow sessions must define a vInf parameter"
191 "for 2D/3D problems.");
192 vSession->LoadParameter(
"vInf",
m_vInf, 0.0);
198 ASSERTL0(vSession->DefinesParameter(
"wInf"),
199 "Compressible flow sessions must define a wInf parameter"
201 vSession->LoadParameter(
"wInf", m_wInf, 0.0);
204 vSession->LoadParameter(
"GasConstant", m_gasConstant, 287.058);
205 vSession->LoadParameter(
"Twall",
m_Twall, 300.15);
206 vSession->LoadSolverInfo(
"ViscosityType", m_ViscosityType,
"Constant");
207 vSession->LoadParameter(
"mu",
m_mu, 1.78e-05);
211 vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
212 vector<vector<NekDouble>> fieldData;
220 vector<vector<LibUtilities::PointsType>> pointsType;
221 for (i = 0; i < fieldDef.size(); ++i)
223 vector<LibUtilities::PointsType> ptype;
224 for (j = 0; j < 2; ++j)
228 pointsType.push_back(ptype);
230 graphShPt->SetExpansionInfo(fieldDef, pointsType);
236 int nfields = vSession->GetVariables().size();
240 for (i = 0; i < pFields.size(); i++)
244 vSession, graphShPt, vSession->GetVariable(i));
253 for (i = 1; i < nfields; ++i)
259 int nSolutionPts = pFields[0]->GetNpoints();
260 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
261 int nElements = pFields[0]->GetExpSize();
277 pFields[0]->GetCoords(x, y,
z);
279 pFields[0]->ExtractTracePhys(x, traceX);
280 pFields[0]->ExtractTracePhys(y, traceY);
281 pFields[0]->ExtractTracePhys(
z, traceZ);
291 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 = 20;
314 for (j = 0; j < nfieldsAdded; ++j)
328 for (i = 0; i < nDimensions; ++i)
332 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
335 for (i = 0; i < nDimensions; ++i)
341 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &traceFieldsAdded[0][0],
345 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &traceFieldsAdded[1][0],
349 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &m_traceTangents[0][0],
351 Vmath::Neg(nTracePts, &m_traceTangents[0][0], 1);
353 Vmath::Vcopy(nTracePts, &m_traceTangents[0][0], 1, &traceFieldsAdded[2][0],
357 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &m_traceTangents[1][0],
360 Vmath::Vcopy(nTracePts, &m_traceTangents[1][0], 1, &traceFieldsAdded[3][0],
371 for (i = 0; i < m_spacedim; i++)
373 Vmath::Vmul(nSolutionPts, &uFields[i + 1][0], 1, &uFields[i + 1][0], 1,
376 Vmath::Smul(nSolutionPts, 0.5, &tmp[0], 1, &tmp[0], 1);
390 pFields[0]->ExtractTracePhys(
pressure, traceFieldsAdded[4]);
402 NekDouble GasConstantInv = 1.0 / m_gasConstant;
403 Vmath::Smul(nSolutionPts, GasConstantInv, &temperature[0], 1,
407 pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[5]);
416 for (i = 0; i < nDimensions; ++i)
422 for (i = 0; i < nDimensions; ++i)
424 for (n = 0; n < nElements; n++)
426 phys_offset = pFields[0]->GetPhys_Offset(n);
428 pFields[i]->GetExp(n)->PhysDeriv(i, temperature + phys_offset,
430 Dtemperature[i] + phys_offset);
433 pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
436 for (i = 0; i < nDimensions; ++i)
439 &traceDtemperature[i][0], 1, &tmp[0], 1);
441 Vmath::Vadd(nTracePts, &traceFieldsAdded[6][0], 1, &tmp[0], 1,
442 &traceFieldsAdded[6][0], 1);
454 for (i = 0; i < nDimensions; ++i)
460 for (i = 0; i < nDimensions; ++i)
462 for (n = 0; n < nElements; n++)
464 phys_offset = pFields[0]->GetPhys_Offset(n);
466 pFields[i]->GetExp(n)->PhysDeriv(i,
pressure + phys_offset,
468 Dpressure[i] + phys_offset);
471 pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
475 for (i = 0; i < nDimensions; ++i)
477 Vmath::Vmul(nTracePts, &m_traceTangents[i][0], 1, &traceDpressure[i][0],
480 Vmath::Vadd(nTracePts, &traceFieldsAdded[7][0], 1, &tmp[0], 1,
481 &traceFieldsAdded[7][0], 1);
485 Vmath::Vcopy(nTracePts, &traceDpressure[0][0], 1, &traceFieldsAdded[8][0],
489 Vmath::Vcopy(nTracePts, &traceDpressure[1][0], 1, &traceFieldsAdded[9][0],
503 for (i = 0; i < nDimensions; ++i)
512 for (j = 0; j < nDimensions; ++j)
519 for (i = 0; i < nDimensions; ++i)
521 for (j = 0; j < nDimensions; ++j)
523 for (n = 0; n < nElements; n++)
525 phys_offset = pFields[0]->GetPhys_Offset(n);
527 pFields[i]->GetExp(n)->PhysDeriv(j,
velocity[i] + phys_offset,
528 auxArray = Dvelocity[i][j] +
533 pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
538 &traceFieldsAdded[10][0], 1);
540 &traceFieldsAdded[11][0], 1);
542 &traceFieldsAdded[12][0], 1);
544 &traceFieldsAdded[13][0], 1);
561 if (m_ViscosityType ==
"Variable")
567 for (
int i = 0; i < nSolutionPts; ++i)
569 ratio = temperature[i] / T_star;
570 mu[i] = mu_star * ratio *
sqrt(ratio) * (T_star + 110.0) /
571 (temperature[i] + 110.0);
584 Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
587 Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[0][0][0], 1, &divVel[0],
589 Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[1][1][0], 1, &divVel[0],
593 Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
594 Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
598 for (j = 0; j < m_spacedim; ++j)
603 Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
606 Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
613 Vmath::Vadd(nSolutionPts, &Dvelocity[0][1][0], 1, &Dvelocity[1][0][0], 1,
617 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
619 pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[14]);
620 pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[15]);
621 pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[16]);
636 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &cosTeta[0], 1);
639 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &sinTeta[0], 1);
642 Vmath::Vsub(nTracePts, &traceFieldsAdded[14][0], 1,
643 &traceFieldsAdded[15][0], 1, &sigma_diff[0], 1);
646 Vmath::Vmul(nTracePts, &cosTeta[0], 1, &sinTeta[0], 1, &tmpTeta[0], 1);
647 Vmath::Smul(nTracePts, 2.0, &tmpTeta[0], 1, &sin2Teta[0], 1);
650 Vmath::Vmul(nTracePts, &cosTeta[0], 1, &cosTeta[0], 1, &cos2Teta[0], 1);
651 Vmath::Vmul(nTracePts, &sinTeta[0], 1, &sinTeta[0], 1, &tmpTeta[0], 1);
652 Vmath::Vsub(nTracePts, &cos2Teta[0], 1, &tmpTeta[0], 1, &cos2Teta[0], 1);
655 Vmath::Smul(nTracePts, -0.5, &sigma_diff[0], 1, &sigma_diff[0], 1);
656 Vmath::Vmul(nTracePts, &sigma_diff[0], 1, &sin2Teta[0], 1, &tau_t[0], 1);
657 Vmath::Vmul(nTracePts, &traceFieldsAdded[16][0], 1, &cos2Teta[0], 1,
659 Vmath::Vadd(nTracePts, &tau_t[0], 1, &tmpTeta[0], 1, &tau_t[0], 1);
661 Vmath::Vcopy(nTracePts, &tau_t[0], 1, &traceFieldsAdded[17][0], 1);
667 pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[18]);
678 Vmath::Smul(nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
679 Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
684 for (
int i = 0; i < m_spacedim; ++i)
686 Vmath::Vvtvp(nSolutionPts, uFields[i + 1], 1, uFields[i + 1], 1, mach,
690 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
691 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
693 Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
695 pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[19]);
700 if (pFields[0]->GetBndCondExpansions().size())
704 nBndRegions = pFields[0]->GetBndCondExpansions().size();
705 for (b = 0; b < nBndRegions; ++b)
707 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
708 for (e = 0; e < nBndEdges; ++e)
710 nBndEdgePts = pFields[0]
711 ->GetBndCondExpansions()[b]
715 id2 = pFields[0]->GetTrace()->GetPhys_Offset(
716 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
719 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
721 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
723 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
726 Vmath::Vcopy(nBndEdgePts, &traceX[id2], 1, &surfaceX[id1],
729 Vmath::Vcopy(nBndEdgePts, &traceY[id2], 1, &surfaceY[id1],
732 Vmath::Vcopy(nBndEdgePts, &traceZ[id2], 1, &surfaceZ[id1],
742 if (pFields[0]->GetBndCondExpansions().size())
744 for (j = 0; j < nfields; ++j)
748 nBndRegions = pFields[j]->GetBndCondExpansions().size();
749 for (b = 0; b < nBndRegions; ++b)
751 nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
752 for (e = 0; e < nBndEdges; ++e)
754 nBndEdgePts = pFields[j]
755 ->GetBndCondExpansions()[b]
759 id2 = pFields[j]->GetTrace()->GetPhys_Offset(
760 pFields[j]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
763 if (pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
765 pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
767 pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
771 &surfaceFields[j][id1], 1);
781 if (pFields[0]->GetBndCondExpansions().size())
783 for (j = 0; j < nfieldsAdded; ++j)
787 nBndRegions = pFields[0]->GetBndCondExpansions().size();
788 for (b = 0; b < nBndRegions; ++b)
790 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
791 for (e = 0; e < nBndEdges; ++e)
793 nBndEdgePts = pFields[0]
794 ->GetBndCondExpansions()[b]
798 id2 = pFields[0]->GetTrace()->GetPhys_Offset(
799 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
802 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
804 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
806 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
810 &surfaceFieldsAdded[j][id1], 1);
821 std::string vEquation = vSession->GetSolverInfo(
"EQType");
824 BndExp = pFields[0]->GetBndCondExpansions();
835 for (
int i = 0; i < BndExp[0]->GetExpSize(); ++i)
858 for (
int j = 0; j < nbc; ++j)
861 nxOnBnd[j] = surfaceFieldsAdded[0][GlobalIndex];
862 nyOnBnd[j] = surfaceFieldsAdded[1][GlobalIndex];
863 txOnBnd[j] = surfaceFieldsAdded[2][GlobalIndex];
864 tyOnBnd[j] = surfaceFieldsAdded[3][GlobalIndex];
866 PressurOnBnd[j] = surfaceFieldsAdded[4][GlobalIndex];
868 if (vEquation ==
"NavierStokesCFE")
870 ShearStressOnBnd[j] = surfaceFieldsAdded[17][GlobalIndex];
890 Vmath::Vmul(nbc, PressurOnBnd, 1, nxOnBnd, 1, drag_p, 1);
891 Vmath::Vmul(nbc, PressurOnBnd, 1, nyOnBnd, 1, lift_p, 1);
896 Fxp += bc->Integral(drag_p);
897 Fyp += bc->Integral(lift_p);
899 if (vEquation ==
"NavierStokesCFE")
901 Vmath::Vmul(nbc, ShearStressOnBnd, 1, txOnBnd, 1, drag_v, 1);
902 Vmath::Vmul(nbc, ShearStressOnBnd, 1, tyOnBnd, 1, lift_v, 1);
907 Fxv += bc->Integral(drag_v);
908 Fyv += bc->Integral(lift_v);
911 Sref += bc->Integral(Unity);
914 cout <<
"\n Sref = " << Sref << endl;
919 cout <<
" Pressure drag (Fxp) = " << Fxp << endl;
920 cout <<
" Pressure lift (Fyp) = " << Fyp << endl;
921 cout <<
" Viscous drag (Fxv) = " << Fxv << endl;
922 cout <<
" Viscous lift (Fyv) = " << Fyv << endl;
923 cout <<
"\n ==> Total drag = " << Fxp + Fxv << endl;
924 cout <<
" ==> Total lift = " << Fyp + Fyv <<
"\n" << endl;
932 outfile.open(fname.c_str());
949 <<
"rhou[kg/(m^2 s)] "
951 <<
"rhov[kg/(m^2 s)] "
989 for (i = 0; i < id1; ++i)
991 outfile << scientific << setw(17) << setprecision(16) << surfaceX[i]
992 <<
" \t " << surfaceY[i] <<
" \t " << surfaceZ[i] <<
" \t "
993 << surfaceFieldsAdded[0][i] <<
" \t "
994 << surfaceFieldsAdded[1][i] <<
" \t "
995 << surfaceFieldsAdded[2][i] <<
" \t "
996 << surfaceFieldsAdded[3][i] <<
" \t " << surfaceFields[0][i]
997 <<
" \t " << surfaceFields[1][i] <<
" \t "
998 << surfaceFields[2][i] <<
" \t " << surfaceFields[3][i]
999 <<
" \t " << surfaceFieldsAdded[4][i] <<
" \t "
1000 << surfaceFieldsAdded[5][i] <<
" \t "
1001 << surfaceFieldsAdded[6][i] <<
" \t "
1002 << surfaceFieldsAdded[7][i] <<
" \t "
1003 << surfaceFieldsAdded[8][i] <<
" \t "
1004 << surfaceFieldsAdded[9][i] <<
" \t "
1005 << surfaceFieldsAdded[10][i] <<
" \t "
1006 << surfaceFieldsAdded[11][i] <<
" \t "
1007 << surfaceFieldsAdded[12][i] <<
" \t "
1008 << surfaceFieldsAdded[13][i] <<
" \t "
1009 << surfaceFieldsAdded[14][i] <<
" \t "
1010 << surfaceFieldsAdded[15][i] <<
" \t "
1011 << surfaceFieldsAdded[16][i] <<
" \t "
1012 << surfaceFieldsAdded[17][i] <<
" \t "
1013 << surfaceFieldsAdded[18][i] <<
" \t "
1014 << surfaceFieldsAdded[19][i]
1019 outfile << endl << endl;
#define ASSERTL0(condition, msg)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
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...
std::map< std::string, std::string > FieldMetaDataMap
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
const std::vector< NekDouble > velocity
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
std::vector< double > z(NPUPPER)
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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.
void Neg(int n, T *x, const int incx)
Negate x = -x.
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 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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
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.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
scalarT< T > sqrt(scalarT< T > in)