38#include <boost/geometry.hpp>
39#include <boost/lexical_cast.hpp>
52namespace bg = boost::geometry;
53namespace bgi = boost::geometry::index;
61 "Interpolates a field to a set of points. Requires fromfld, fromxml "
62 "to be defined, and a topts, line, plane or block of target points ");
67 false,
"NotSet",
"Xml file from which to interpolate field");
70 false,
"NotSet",
"Fld file from which to interpolate field");
73 ConfigOption(
false,
"NotSet",
"Pts file to which interpolate field");
75 "Specify a line of N points using "
76 "line=N,x0,y0,z0,z1,y1,z1");
79 "Specify a plane of N1 x N2 points using "
80 "plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,y3,z3");
83 "Specify a rectangular box of N1 x N2 x N3 points "
84 "using a box of points limited by box="
85 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
88 ConfigOption(
false,
"-10000000",
"Lower bound for interpolation value");
90 ConfigOption(
false,
"10000000",
"Upper bound for interpolation value");
92 ConfigOption(
false,
"0",
"Default value if point is outside domain");
96 "Parameters p0 and q to determine pressure coefficients");
98 ConfigOption(
false,
"NotSet",
"Take fields as sin mode");
112 std::vector<std::string> files;
116 char *argv[] = {
const_cast<char *
>(
"FieldConvert"),
nullptr};
125 int coordim =
m_f->m_fieldPts->GetDim();
126 int npts =
m_f->m_fieldPts->GetNpoints();
127 std::vector<std::string> fieldNames =
m_f->m_fieldPts->GetFieldNames();
128 for (
auto &it : fieldNames)
130 m_f->m_fieldPts->RemoveField(it);
134 m_f->m_fieldPts->GetPts(pts);
136 rng->m_checkShape =
false;
144 rng->m_doZrange =
true;
147 if (rng->m_zmax == rng->m_zmin)
154 rng->m_doYrange =
true;
159 rng->m_doXrange =
true;
173 fromField->m_graph->GetExpansionInfo();
177 for (
auto &expIt : expansions)
179 ElementGIDs[i++] = expIt.second->m_geomShPtr->GetGlobalID();
183 ASSERTL0(i > 0,
"No elements are set. Are the interpolated points "
184 "within the domain given by the xml files?");
185 string fromfld =
m_config[
"fromfld"].as<
string>();
186 m_f->FieldIOForFile(fromfld)->Import(
187 fromfld, fromField->m_fielddef, fromField->m_data,
189 int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
190 for (i = 0; i < fromField->m_fielddef.size(); ++i)
192 int d1 = fromField->m_fielddef[i]->m_basis.size();
194 if (d1 >= 0 && (fromField->m_fielddef[i]->m_basis[d1] ==
196 fromField->m_fielddef[i]->m_basis[d1] ==
199 fromField->m_fielddef[i]->m_homogeneousZIDs[0] += 2;
200 fromField->m_fielddef[i]->m_numModes[d1] = 4;
207 fromField->m_graph->SetExpansionInfo(fromField->m_fielddef);
208 int nfields = fromField->m_fielddef[0]->m_fields.size();
209 fromField->m_exp.resize(nfields);
210 fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,
true);
211 m_f->m_exp.resize(nfields);
214 for (i = 1; i < nfields; ++i)
216 fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
221 if (
m_config[
"realmodetoimag"].as<string>().compare(
"NotSet"))
224 m_f->m_variables, sinmode);
226 for (
int j = 0; j < nfields; ++j)
228 for (i = 0; i < fromField->m_fielddef.size(); i++)
230 fromField->m_exp[j]->ExtractDataToCoeffs(
231 fromField->m_fielddef[i], fromField->m_data[i],
232 fromField->m_fielddef[0]->m_fields[j],
233 fromField->m_exp[j]->UpdateCoeffs());
235 if (NumHomogeneousDir == 1)
237 fromField->m_exp[j]->SetWaveSpace(
true);
238 if (sinmode.count(j))
240 int Ncoeff = fromField->m_exp[j]->GetPlane(2)->GetNcoeffs();
242 Ncoeff, -1., fromField->m_exp[j]->GetPlane(2)->GetCoeffs(),
243 1, fromField->m_exp[j]->GetPlane(3)->UpdateCoeffs(), 1);
245 fromField->m_exp[j]->GetPlane(2)->UpdateCoeffs(),
249 fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
250 fromField->m_exp[j]->UpdatePhys());
253 m_f->m_fieldPts->AddField(newPts,
254 fromField->m_fielddef[0]->m_fields[j]);
255 m_f->m_variables.push_back(fromField->m_fielddef[0]->m_fields[j]);
263 clamp_up, def_value);
265 if (!boost::iequals(
m_config[
"cp"].as<string>(),
"NotSet"))
273 int rank =
m_f->m_comm->GetSpaceComm()->GetRank();
274 int nprocs =
m_f->m_comm->GetSpaceComm()->GetSize();
277 if (
m_config[
"topts"].as<string>().compare(
"NotSet") != 0)
279 string inFile =
m_config[
"topts"].as<
string>();
281 if (fs::path(inFile).extension() ==
".pts")
287 ptsIO->Import(inFile,
m_f->m_fieldPts);
289 else if (fs::path(inFile).extension() ==
".csv")
295 csvIO->Import(inFile,
m_f->m_fieldPts);
299 ASSERTL0(
false,
"unknown topts file type");
302 else if (
m_config[
"line"].as<string>().compare(
"NotSet") != 0)
304 vector<NekDouble> values;
307 "Failed to interpret line string");
309 ASSERTL0(values.size() > 2,
"line string should contain 2*Dim+1 values "
310 "N,x0,y0,z0,x1,y1,z1");
313 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
314 ASSERTL0(values[0] > 1,
"N is not a valid number");
316 int dim = (values.size() - 1) / 2;
317 int npts = values[0];
320 int ptsPerProc = npts / nprocs;
321 int extraPts = (rank < nprocs - 1) ? 0 : npts % nprocs;
322 int locPts = ptsPerProc + extraPts;
323 int start = rank * ptsPerProc;
324 int end = start + locPts;
328 for (
int i = 0; i < dim; ++i)
331 delta[i] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
334 for (
int i = 0, cntLoc = 0; i < npts; ++i)
336 if (i >= start && i < end)
338 for (
int n = 0; n < dim; ++n)
340 pts[n][cntLoc] = values[n + 1] + i * delta[n];
351 m_f->m_fieldPts->SetPointsPerEdge(ppe);
353 else if (
m_config[
"plane"].as<string>().compare(
"NotSet") != 0)
355 vector<NekDouble> values;
358 "Failed to interpret plane string");
361 "plane string should contain 4 Dim+2 values "
362 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
365 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
366 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
368 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
369 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
371 int dim = (values.size() - 2) / 4;
377 int totpts = npts[0] * npts[1];
380 int ptsPerProc = totpts / nprocs;
381 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
382 int locPts = ptsPerProc + extraPts;
383 int start = rank * ptsPerProc;
384 int end = start + locPts;
389 for (
int i = 0; i < dim; ++i)
392 delta1[i] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
394 delta2[i] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
398 for (
int j = 0, cnt = 0, cntLoc = 0; j < npts[1]; ++j)
400 for (
int i = 0; i < npts[0]; ++i, ++cnt)
402 if (cnt >= start && cnt < end)
404 for (
int n = 0; n < dim; ++n)
407 (values[2 + n] + i * delta1[n]) *
409 (values[2 + 3 * dim + n] + i * delta2[n]) *
418 ppe.push_back(npts[0]);
419 ppe.push_back(npts[1]);
423 m_f->m_fieldPts->SetPointsPerEdge(ppe);
425 else if (
m_config[
"box"].as<string>().compare(
"NotSet") != 0)
427 vector<NekDouble> values;
430 "Failed to interpret box string");
432 ASSERTL0(values.size() == 9,
"box string should contain 9 values "
433 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
441 int totpts = npts[0] * npts[1] * npts[2];
447 int ptsPerProc = totpts / nprocs;
448 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
449 int locPts = ptsPerProc + extraPts;
450 int start = rank * ptsPerProc;
451 int end = start + locPts;
453 for (
int i = 0; i < dim; ++i)
456 delta[i] = (values[4 + 2 * i] - values[3 + 2 * i]) / (npts[i] - 1);
459 for (
int k = 0, cnt = 0, cntLoc = 0; k < npts[2]; ++k)
461 for (
int j = 0; j < npts[1]; ++j)
463 for (
int i = 0; i < npts[0]; ++i, ++cnt)
465 if (cnt >= start && cnt < end)
467 pts[0][cntLoc] = values[3] + i * delta[0];
468 pts[1][cntLoc] = values[5] + j * delta[1];
469 pts[2][cntLoc] = values[7] + k * delta[2];
477 ppe.push_back(npts[0]);
478 ppe.push_back(npts[1]);
479 ppe.push_back(npts[2]);
483 m_f->m_fieldPts->SetPointsPerEdge(ppe);
484 vector<NekDouble> boxdim;
485 boxdim.assign(&values[3], &values[3] + 6);
486 m_f->m_fieldPts->SetBoxSize(boxdim);
490 ASSERTL0(
false,
"Missing target points for ProcessInterpPoints.");
495 vector<MultiRegions::ExpListSharedPtr> &field0,
499 ASSERTL0(pts->GetNFields() == field0.size(),
"ptField has too few fields");
501 int nfields = field0.size();
504 if (
m_f->m_comm->GetRank() == 0)
511 if (
m_f->m_comm->GetRank() == 0)
516 for (
int f = 0; f < nfields; ++f)
518 for (
int i = 0; i < pts->GetNpoints(); ++i)
520 if (pts->GetPointVal(f, i) > clamp_up)
522 pts->SetPointVal(f, i, clamp_up);
524 else if (pts->GetPointVal(f, i) < clamp_low)
526 pts->SetPointVal(f, i, clamp_low);
535 int dim = pts->GetDim();
536 int nq1 = pts->GetNpoints();
542 vector<NekDouble> values;
544 "Failed to interpret cp string");
546 ASSERTL0(values.size() == 2,
"cp string should contain 2 values "
547 "p0 and q (=1/2 rho u^2)");
550 qinv = 1.0 / values[1];
552 for (
int i = 0; i < pts->GetNFields(); ++i)
554 if (boost::iequals(pts->GetFieldName(i),
"p"))
559 if (boost::iequals(pts->GetFieldName(i),
"u") ||
560 boost::iequals(pts->GetFieldName(i),
"v") ||
561 boost::iequals(pts->GetFieldName(i),
"w"))
571 WARNINGL0(
false,
"Did not find velocity components for Cp0");
576 WARNINGL0(
false,
"Failed to find 'p' field to determine cp0");
582 for (f = 0; f < 2; ++f)
587 for (r = 0; r < nq1; r++)
591 data[0][r] = qinv * (pts->GetPointVal(dim + pfield, r) - p0);
596 for (
int i = 0; i < velid.size(); ++i)
598 q += 0.5 * pts->GetPointVal(dim + velid[i], r) *
599 pts->GetPointVal(dim + velid[i], r);
602 qinv * (pts->GetPointVal(dim + pfield, r) +
q - p0);
609 pts->AddField(data[0],
"Cp");
610 m_f->m_variables.push_back(
"Cp");
613 pts->AddField(data[1],
"Cp0");
614 m_f->m_variables.push_back(
"Cp0");
620 const int goal)
const
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define WARNINGL0(condition, msg)
A class that contains algorithms for interpolation between pts fields, expansions and different meshe...
FIELD_UTILS_EXPORT void Interpolate(const T expInField, T &expOutField, NekDouble def_value=0.0)
Interpolate from an expansion to an expansion.
FieldSharedPtr m_f
Field object.
std::map< std::string, ConfigOption > m_config
List of configuration values.
void InterpolateFieldToPts(std::vector< MultiRegions::ExpListSharedPtr > &field0, LibUtilities::PtsFieldSharedPtr &pts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value)
~ProcessInterpPoints() override
void CreateFieldPts(po::variables_map &vm)
void PrintProgressbar(const int position, const int goal) const
static ModuleKey className
ProcessInterpPoints(FieldSharedPtr f)
static std::shared_ptr< Module > create(FieldSharedPtr f)
Creates an instance of this class.
void v_Process(po::variables_map &vm) override
Write mesh to output file.
Abstract base class for processing modules.
void SetProgressCallback(FuncPointerT func, ObjectPointerT obj)
sets a callback funtion which gets called every time the interpolation progresses
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static bool GenerateVariableSet(const std::string &str, const std::vector< std::string > &variables, std::set< int > &out)
Generate a set of variable locations.
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
std::shared_ptr< Field > FieldSharedPtr
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
int PrintProgressbar(const int position, const int goal, const std::string message, int lastprogress=-1)
Prints a progressbar.
static FieldMetaDataMap NullFieldMetaDataMap
std::shared_ptr< PtsField > PtsFieldSharedPtr
std::shared_ptr< DomainRange > DomainRangeShPtr
std::shared_ptr< CsvIO > CsvIOSharedPtr
std::shared_ptr< PtsIO > PtsIOSharedPtr
CommFactory & GetCommFactory()
@ eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
@ eFourier
Fourier Expansion .
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
std::vector< double > q(NPUPPER *NPUPPER)
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Represents a command-line configuration option.