38#include <boost/geometry.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");
100 ConfigOption(
false,
"1e-5",
"The maximum acceptable distance.");
114 std::vector<std::string> files;
118 char *argv[] = {
const_cast<char *
>(
"FieldConvert"),
nullptr};
127 int coordim =
m_f->m_fieldPts->GetDim();
128 int npts =
m_f->m_fieldPts->GetNpoints();
129 std::vector<std::string> fieldNames =
m_f->m_fieldPts->GetFieldNames();
130 for (
auto &it : fieldNames)
132 m_f->m_fieldPts->RemoveField(it);
136 m_f->m_fieldPts->GetPts(pts);
138 rng->m_checkShape =
false;
146 rng->m_doZrange =
true;
149 if (rng->m_zmax == rng->m_zmin)
156 rng->m_doYrange =
true;
161 rng->m_doXrange =
true;
175 fromField->m_graph->GetExpansionInfo();
179 for (
auto &expIt : expansions)
181 ElementGIDs[i++] = expIt.second->m_geomShPtr->GetGlobalID();
185 ASSERTL0(i > 0,
"No elements are set. Are the interpolated points "
186 "within the domain given by the xml files?");
187 string fromfld =
m_config[
"fromfld"].as<
string>();
188 m_f->FieldIOForFile(fromfld)->Import(
189 fromfld, fromField->m_fielddef, fromField->m_data,
191 int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
192 for (i = 0; i < fromField->m_fielddef.size(); ++i)
194 int d1 = fromField->m_fielddef[i]->m_basis.size();
196 if (d1 >= 0 && (fromField->m_fielddef[i]->m_basis[d1] ==
198 fromField->m_fielddef[i]->m_basis[d1] ==
201 fromField->m_fielddef[i]->m_homogeneousZIDs[0] += 2;
202 fromField->m_fielddef[i]->m_numModes[d1] = 4;
209 fromField->m_graph->SetExpansionInfo(fromField->m_fielddef);
210 int nfields = fromField->m_fielddef[0]->m_fields.size();
211 fromField->m_exp.resize(nfields);
212 fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,
true);
213 m_f->m_exp.resize(nfields);
216 for (i = 1; i < nfields; ++i)
218 fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
223 if (
m_config[
"realmodetoimag"].as<string>().compare(
"NotSet"))
226 m_f->m_variables, sinmode);
228 for (
int j = 0; j < nfields; ++j)
230 for (i = 0; i < fromField->m_fielddef.size(); i++)
232 fromField->m_exp[j]->ExtractDataToCoeffs(
233 fromField->m_fielddef[i], fromField->m_data[i],
234 fromField->m_fielddef[0]->m_fields[j],
235 fromField->m_exp[j]->UpdateCoeffs());
237 if (NumHomogeneousDir == 1)
239 fromField->m_exp[j]->SetWaveSpace(
true);
240 if (sinmode.count(j))
242 int Ncoeff = fromField->m_exp[j]->GetPlane(2)->GetNcoeffs();
244 Ncoeff, -1., fromField->m_exp[j]->GetPlane(2)->GetCoeffs(),
245 1, fromField->m_exp[j]->GetPlane(3)->UpdateCoeffs(), 1);
247 fromField->m_exp[j]->GetPlane(2)->UpdateCoeffs(),
251 fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
252 fromField->m_exp[j]->UpdatePhys());
255 m_f->m_fieldPts->AddField(newPts,
256 fromField->m_fielddef[0]->m_fields[j]);
257 m_f->m_variables.push_back(fromField->m_fielddef[0]->m_fields[j]);
266 clamp_up, def_value, tolerance);
268 if (!boost::iequals(
m_config[
"cp"].as<string>(),
"NotSet"))
276 int rank =
m_f->m_comm->GetSpaceComm()->GetRank();
277 int nprocs =
m_f->m_comm->GetSpaceComm()->GetSize();
280 if (
m_config[
"topts"].as<string>().compare(
"NotSet") != 0)
282 string inFile =
m_config[
"topts"].as<
string>();
284 if (fs::path(inFile).extension() ==
".pts")
290 ptsIO->Import(inFile,
m_f->m_fieldPts);
292 else if (fs::path(inFile).extension() ==
".csv")
298 csvIO->Import(inFile,
m_f->m_fieldPts);
302 ASSERTL0(
false,
"unknown topts file type");
305 else if (
m_config[
"line"].as<string>().compare(
"NotSet") != 0)
307 vector<NekDouble> values;
310 "Failed to interpret line string");
312 ASSERTL0(values.size() > 2,
"line string should contain 2*Dim+1 values "
313 "N,x0,y0,z0,x1,y1,z1");
316 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
317 ASSERTL0(values[0] > 1,
"N is not a valid number");
319 int dim = (values.size() - 1) / 2;
320 int npts = values[0];
323 int ptsPerProc = npts / nprocs;
324 int extraPts = (rank < nprocs - 1) ? 0 : npts % nprocs;
325 int locPts = ptsPerProc + extraPts;
326 int start = rank * ptsPerProc;
327 int end = start + locPts;
331 for (
int i = 0; i < dim; ++i)
334 delta[i] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
337 for (
int i = 0, cntLoc = 0; i < npts; ++i)
339 if (i >= start && i < end)
341 for (
int n = 0; n < dim; ++n)
343 pts[n][cntLoc] = values[n + 1] + i * delta[n];
354 m_f->m_fieldPts->SetPointsPerEdge(ppe);
356 else if (
m_config[
"plane"].as<string>().compare(
"NotSet") != 0)
358 vector<NekDouble> values;
361 "Failed to interpret plane string");
364 "plane string should contain 4 Dim+2 values "
365 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
368 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
369 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
371 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
372 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
374 int dim = (values.size() - 2) / 4;
380 int totpts = npts[0] * npts[1];
383 int ptsPerProc = totpts / nprocs;
384 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
385 int locPts = ptsPerProc + extraPts;
386 int start = rank * ptsPerProc;
387 int end = start + locPts;
392 for (
int i = 0; i < dim; ++i)
395 delta1[i] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
397 delta2[i] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
401 for (
int j = 0, cnt = 0, cntLoc = 0; j < npts[1]; ++j)
403 for (
int i = 0; i < npts[0]; ++i, ++cnt)
405 if (cnt >= start && cnt < end)
407 for (
int n = 0; n < dim; ++n)
410 (values[2 + n] + i * delta1[n]) *
412 (values[2 + 3 * dim + n] + i * delta2[n]) *
421 ppe.push_back(npts[0]);
422 ppe.push_back(npts[1]);
426 m_f->m_fieldPts->SetPointsPerEdge(ppe);
428 else if (
m_config[
"box"].as<string>().compare(
"NotSet") != 0)
430 vector<NekDouble> values;
433 "Failed to interpret box string");
435 ASSERTL0(values.size() == 9,
"box string should contain 9 values "
436 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
444 int totpts = npts[0] * npts[1] * npts[2];
450 int ptsPerProc = totpts / nprocs;
451 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
452 int locPts = ptsPerProc + extraPts;
453 int start = rank * ptsPerProc;
454 int end = start + locPts;
456 for (
int i = 0; i < dim; ++i)
459 delta[i] = (values[4 + 2 * i] - values[3 + 2 * i]) / (npts[i] - 1);
462 for (
int k = 0, cnt = 0, cntLoc = 0; k < npts[2]; ++k)
464 for (
int j = 0; j < npts[1]; ++j)
466 for (
int i = 0; i < npts[0]; ++i, ++cnt)
468 if (cnt >= start && cnt < end)
470 pts[0][cntLoc] = values[3] + i * delta[0];
471 pts[1][cntLoc] = values[5] + j * delta[1];
472 pts[2][cntLoc] = values[7] + k * delta[2];
480 ppe.push_back(npts[0]);
481 ppe.push_back(npts[1]);
482 ppe.push_back(npts[2]);
486 m_f->m_fieldPts->SetPointsPerEdge(ppe);
487 vector<NekDouble> boxdim;
488 boxdim.assign(&values[3], &values[3] + 6);
489 m_f->m_fieldPts->SetBoxSize(boxdim);
493 ASSERTL0(
false,
"Missing target points for ProcessInterpPoints.");
498 vector<MultiRegions::ExpListSharedPtr> &field0,
502 ASSERTL0(pts->GetNFields() == field0.size(),
"ptField has too few fields");
504 int nfields = field0.size();
507 if (
m_f->m_comm->GetRank() == 0)
512 interp.
Interpolate(field0, pts, def_value, tolerance);
514 if (
m_f->m_comm->GetRank() == 0)
519 for (
int f = 0; f < nfields; ++f)
521 for (
int i = 0; i < pts->GetNpoints(); ++i)
523 if (pts->GetPointVal(f, i) > clamp_up)
525 pts->SetPointVal(f, i, clamp_up);
527 else if (pts->GetPointVal(f, i) < clamp_low)
529 pts->SetPointVal(f, i, clamp_low);
538 int dim = pts->GetDim();
539 int nq1 = pts->GetNpoints();
545 vector<NekDouble> values;
547 "Failed to interpret cp string");
549 ASSERTL0(values.size() == 2,
"cp string should contain 2 values "
550 "p0 and q (=1/2 rho u^2)");
553 qinv = 1.0 / values[1];
555 for (
int i = 0; i < pts->GetNFields(); ++i)
557 if (boost::iequals(pts->GetFieldName(i),
"p"))
562 if (boost::iequals(pts->GetFieldName(i),
"u") ||
563 boost::iequals(pts->GetFieldName(i),
"v") ||
564 boost::iequals(pts->GetFieldName(i),
"w"))
574 WARNINGL0(
false,
"Did not find velocity components for Cp0");
579 WARNINGL0(
false,
"Failed to find 'p' field to determine cp0");
585 for (f = 0; f < 2; ++f)
590 for (r = 0; r < nq1; r++)
594 data[0][r] = qinv * (pts->GetPointVal(dim + pfield, r) - p0);
599 for (
int i = 0; i < velid.size(); ++i)
601 q += 0.5 * pts->GetPointVal(dim + velid[i], r) *
602 pts->GetPointVal(dim + velid[i], r);
605 qinv * (pts->GetPointVal(dim + pfield, r) +
q - p0);
612 pts->AddField(data[0],
"Cp");
613 m_f->m_variables.push_back(
"Cp");
616 pts->AddField(data[1],
"Cp0");
617 m_f->m_variables.push_back(
"Cp0");
623 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., NekDouble tolerance=NekConstants::kFindDistanceMin)
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, NekDouble tolerance)
~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.