42#include <boost/format.hpp>
55 const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams)
56 :
Filter(pSession, pEquation)
59 auto it = pParams.find(
"OutputFile");
60 if (it == pParams.end())
66 ASSERTL0(it->second.length() > 0,
"Missing parameter 'OutputFile'.");
76 it = pParams.find(
"OutputFrequency");
77 if (it == pParams.end())
88 it = pParams.find(
"OutputOneFile");
89 if (it == pParams.end())
95 std::string sOption = it->second.c_str();
97 (boost::iequals(sOption,
"yes"));
104 it = pParams.find(
"OutputPlane");
105 if (it == pParams.end())
115 it = pParams.find(
"WaveSpace");
116 if (it == pParams.end())
122 std::string sOption = it->second.c_str();
124 (boost::iequals(sOption,
"yes"));
129 if (pParams.end() != (it = pParams.find(
"Points")))
134 else if (pParams.end() != (it = pParams.find(
"line")))
136 vector<NekDouble> values;
138 "Failed to interpret line string");
140 ASSERTL0(values.size() > 2,
"line string should contain 2*Dim+1 values "
141 "N,x0,y0,z0,x1,y1,z1");
144 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
145 ASSERTL0(values[0] > 1,
"N is not a valid number");
147 int dim = (values.size() - 1) / 2;
148 int npts = values[0];
152 for (
int i = 0; i < dim; ++i)
154 delta[i + 0] = values[i + 1];
155 delta[i + 3] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
160 else if (pParams.end() != (it = pParams.find(
"plane")))
162 vector<NekDouble> values;
164 "Failed to interpret plane string");
167 "plane string should contain 4 Dim+2 values "
168 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
171 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
172 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
174 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
175 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
177 int dim = (values.size() - 2) / 4;
181 npts[1] = values[0] * values[1];
185 for (
int i = 0; i < dim; ++i)
187 delta[i + 0] = values[2 + i];
188 delta[i + 3] = values[2 + 3 * dim + i];
189 delta[i + 6] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
191 delta[i + 9] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
197 else if (pParams.end() != (it = pParams.find(
"box")))
199 vector<NekDouble> values;
201 "Failed to interpret box string");
203 ASSERTL0(values.size() == 9,
"box string should contain 9 values "
204 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
209 npts[1] = values[0] * values[1];
210 npts[2] = values[0] * values[1] * values[2];
213 for (
int i = 0; i < dim; ++i)
215 delta[i + 0] = values[3 + 2 * i];
217 (values[4 + 2 * i] - values[3 + 2 * i]) / (values[i] - 1);
224 ASSERTL0(
false,
"Missing parameter 'Points'.");
251 for (
int n = 0; n < 3; ++n)
253 gloCoord[n] = values[n] + I * delta[n];
270 for (
int n = 0; n < 3; ++n)
272 gloCoord[n] = (values[n] + i[0] * delta[n]) * (1.0 - i[1] / n1) +
273 (values[3 + n] + i[0] * delta[3 + n]) * (i[1] / n1);
290 for (
int n = 0; n < 3; ++n)
292 gloCoord[n] = values[n] + i[n] * delta[n];
312 std::vector<Array<OneD, const NekDouble>> GloCoords;
313 const int coordim = pFields[0]->GetGraph()->GetSpaceDimension();
322 int nplanes = pFields[0]->GetHomogeneousBasis()->GetZ().size();
323 NekDouble lhom = pFields[0]->GetHomoLen();
332 plane = round((gloCoord[coordim] * nplanes) / lhom);
338 NekDouble Z = (pFields[0]->GetHomogeneousBasis()->GetZ())[plane];
339 Z = (Z + 1) * lhom / 2;
340 if (fabs(gloCoord[coordim] - Z) >
342 vComm->GetRank() == 0)
344 cout <<
"Resetting History point from z = " << gloCoord[coordim]
345 <<
" to z = " << Z << endl;
347 gloCoord[coordim] = Z;
356 const int vRank = vComm->GetRank();
357 const int vColumnRank = vComm->GetColumnComm()->GetRank();
358 const int vRowRank = vComm->GetRowComm()->GetRank();
364 std::vector<Array<OneD, const NekDouble>> LocCoords;
371 gloCoord = GloCoords[i];
376 idList[i] = pFields[0]->GetPlane(0)->GetExpIndex(
381 idList[i] = pFields[0]->GetExpIndex(gloCoord, locCoord,
385 for (
int j = 0; j < locCoord.size(); ++j)
387 locCoord[j] = std::max(locCoord[j], -1.0);
388 locCoord[j] = std::min(locCoord[j], 1.0);
393 if (idList[i] != -1 && vColumnRank == 0)
417 if (procList[i] == vRowRank)
424 std::make_tuple(GloCoords[i], LocCoords[i], i, idList[i]));
429 if (vComm->TreatAsRankZero())
444 if (vComm->TreatAsRankZero())
453 boost::lexical_cast<std::string>(gloCoord[0]) +
", " +
454 boost::lexical_cast<std::string>(gloCoord[1]) +
", " +
455 boost::lexical_cast<std::string>(gloCoord[2]) +
456 " cannot be found in the mesh.");
477 const int numFields = pFields.size();
478 const int coordim = pFields[0]->GetGraph()->GetSpaceDimension();
486 for (
int j = 0; j < numFields; ++j)
491 int nPlanes = pFields[j]->GetHomogeneousBasis()->GetZ().size();
492 NekDouble lHom = pFields[j]->GetHomoLen();
495 ASSERTL0(pFields[j]->GetWaveSpace(),
496 "HistoryPoints in Homogeneous1D require that solution "
498 gloCoord = std::get<0>(x);
499 locCoord = std::get<1>(x);
500 const int indx = std::get<2>(x);
501 const int expId = std::get<3>(x);
505 2. * M_PI * fmod(gloCoord[coordim], lHom) / lHom;
506 for (
size_t n = 0; n < planes.size(); ++n)
512 physvals = pFields[j]->GetPlane(n)->UpdatePhys() +
513 pFields[j]->GetPhys_Offset(expId);
516 if (pFields[j]->GetPhysState() ==
false)
518 pFields[j]->GetPlane(n)->GetExp(expId)->BwdTrans(
519 pFields[j]->GetPlane(n)->GetCoeffs() +
520 pFields[j]->GetCoeff_Offset(expId),
525 pFields[j]->GetPlane(n)->GetExp(expId)->StdPhysEvaluate(
538 else if (planes[n] == 1)
540 value += cos(0.5 * nPlanes * BetaT) * coeff;
542 else if (planes[n] % 2 == 0)
544 NekDouble phase = (planes[n] >> 1) * BetaT;
545 value += cos(phase) * coeff;
549 NekDouble phase = (planes[n] >> 1) * BetaT;
550 value += -sin(phase) * coeff;
555 data[indx * numFields + j] = value;
562 locCoord = std::get<1>(x);
563 const int indx = std::get<2>(x);
564 const int expId = std::get<3>(x);
566 physvals = pFields[j]->UpdatePhys() +
567 pFields[j]->GetPhys_Offset(expId);
570 if (pFields[j]->GetPhysState() ==
false)
572 pFields[j]->GetExp(expId)->BwdTrans(
573 pFields[j]->GetCoeffs() +
574 pFields[j]->GetCoeff_Offset(expId),
579 data[indx * numFields + j] =
580 pFields[j]->GetExp(expId)->StdPhysEvaluate(locCoord,
591 v_WriteData(vComm->GetRank(), data, numFields, time);
596 const int &numFields,
605 std::stringstream vOutputFilename;
627 for (
int i = 0; i < numFields; ++i)
665 for (
int j = 0; j < numFields; ++j)
#define ASSERTL0(condition, msg)
NekDouble Evaluate() const
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
SOLVER_UTILS_EXPORT void v_WriteData(const int &rank, const Array< OneD, NekDouble > &data, const int &numFields, const NekDouble &time)
std::stringstream m_historyPointStream
size_t m_historyPointsSize
SOLVER_UTILS_EXPORT ~FilterHistoryPoints() override
unsigned int m_outputFrequency
Array< OneD, Array< OneD, const NekDouble > > m_historyPoints
bool GetPoint(Array< OneD, NekDouble > gloCoord, int I)
static std::string className
Name of the class.
std::map< LibUtilities::PtsType, Array< OneD, int > > m_pointNumMap
SOLVER_UTILS_EXPORT void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
std::ofstream m_outputStream
int m_outputPlane
plane to take history point from if using a homogeneous1D expansion
unsigned int m_outputIndex
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
SOLVER_UTILS_EXPORT void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
std::list< std::tuple< Array< OneD, const NekDouble >, Array< OneD, const NekDouble >, int, int > > m_historyList
SOLVER_UTILS_EXPORT bool v_IsTimeDependent() override
std::vector< int > m_planeIDs
std::map< LibUtilities::PtsType, Array< OneD, NekDouble > > m_pointDatMap
SOLVER_UTILS_EXPORT FilterHistoryPoints(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
LibUtilities::SessionReaderSharedPtr m_session
std::map< std::string, std::string > ParamMap
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
static const NekDouble kVertexTheSameDouble
static const NekDouble kGeomFactorsTol
FilterFactory & GetFilterFactory()