39 #include <boost/core/ignore_unused.hpp>
40 #include <boost/geometry.hpp>
41 #include <boost/lexical_cast.hpp>
42 #include <boost/math/special_functions/fpclassify.hpp>
53 namespace bg = boost::geometry;
54 namespace bgi = boost::geometry::index;
61 ModuleKey ProcessInterpPoints::className =
64 "Interpolates a field to a set of points. Requires fromfld, fromxml "
65 "to be defined, and a topts, line, plane or block of target points ");
70 false,
"NotSet",
"Xml file from which to interpolate field");
73 false,
"NotSet",
"Fld file from which to interpolate field");
76 ConfigOption(
false,
"NotSet",
"Pts file to which interpolate field");
78 "Specify a line of N points using "
79 "line=N,x0,y0,z0,z1,y1,z1");
82 "Specify a plane of N1 x N2 points using "
83 "plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,y3,z3");
86 "Specify a rectangular box of N1 x N2 x N3 points "
87 "using a box of points limited by box="
88 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
91 ConfigOption(
false,
"-10000000",
"Lower bound for interpolation value");
93 ConfigOption(
false,
"10000000",
"Upper bound for interpolation value");
95 ConfigOption(
false,
"0",
"Default value if point is outside domain");
99 "Parameters p0 and q to determine pressure coefficients");
101 ConfigOption(
false,
"NotSet",
"Take fields as sin mode");
115 std::vector<std::string> files;
119 char *argv[] = {
const_cast<char *
>(
"FieldConvert"),
nullptr};
128 int coordim =
m_f->m_fieldPts->GetDim();
129 int npts =
m_f->m_fieldPts->GetNpoints();
130 std::vector<std::string> fieldNames =
m_f->m_fieldPts->GetFieldNames();
131 for (
auto &it : fieldNames)
133 m_f->m_fieldPts->RemoveField(it);
137 m_f->m_fieldPts->GetPts(pts);
139 rng->m_checkShape =
false;
147 rng->m_doZrange =
true;
150 if (rng->m_zmax == rng->m_zmin)
157 rng->m_doYrange =
true;
162 rng->m_doXrange =
true;
176 fromField->m_graph->GetExpansionInfo();
180 for (
auto &expIt : expansions)
182 ElementGIDs[i++] = expIt.second->m_geomShPtr->GetGlobalID();
186 ASSERTL0(i > 0,
"No elements are set. Are the interpolated points "
187 "within the domain given by the xml files?");
188 string fromfld =
m_config[
"fromfld"].as<
string>();
189 m_f->FieldIOForFile(fromfld)->Import(
190 fromfld, fromField->m_fielddef, fromField->m_data,
192 int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
193 for (i = 0; i < fromField->m_fielddef.size(); ++i)
195 int d1 = fromField->m_fielddef[i]->m_basis.size();
197 if (d1 >= 0 && (fromField->m_fielddef[i]->m_basis[d1] ==
199 fromField->m_fielddef[i]->m_basis[d1] ==
202 fromField->m_fielddef[i]->m_homogeneousZIDs[0] += 2;
203 fromField->m_fielddef[i]->m_numModes[d1] = 4;
210 fromField->m_graph->SetExpansionInfo(fromField->m_fielddef);
211 int nfields = fromField->m_fielddef[0]->m_fields.size();
212 fromField->m_exp.resize(nfields);
213 fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,
true);
214 m_f->m_exp.resize(nfields);
217 for (i = 1; i < nfields; ++i)
219 fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
223 if (
m_config[
"realmodetoimag"].as<string>().compare(
"NotSet"))
227 m_config[
"realmodetoimag"].as<string>(), value),
228 "Failed to interpret realmodetoimag string");
234 for (
int j = 0; j < nfields; ++j)
236 for (i = 0; i < fromField->m_fielddef.size(); i++)
238 fromField->m_exp[j]->ExtractDataToCoeffs(
239 fromField->m_fielddef[i], fromField->m_data[i],
240 fromField->m_fielddef[0]->m_fields[j],
241 fromField->m_exp[j]->UpdateCoeffs());
243 if (NumHomogeneousDir == 1)
245 fromField->m_exp[j]->SetWaveSpace(
true);
246 if (sinmode.count(j))
248 int Ncoeff = fromField->m_exp[j]->GetPlane(2)->GetNcoeffs();
250 Ncoeff, -1., fromField->m_exp[j]->GetPlane(2)->GetCoeffs(),
251 1, fromField->m_exp[j]->GetPlane(3)->UpdateCoeffs(), 1);
253 fromField->m_exp[j]->GetPlane(2)->UpdateCoeffs(),
257 fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
258 fromField->m_exp[j]->UpdatePhys());
260 m_f->m_fieldPts->AddField(newPts,
261 fromField->m_fielddef[0]->m_fields[j]);
262 m_f->m_variables.push_back(fromField->m_fielddef[0]->m_fields[j]);
270 clamp_up, def_value);
272 if (!boost::iequals(
m_config[
"cp"].as<string>(),
"NotSet"))
280 boost::ignore_unused(vm);
282 int rank =
m_f->m_comm->GetRank();
283 int nprocs =
m_f->m_comm->GetSize();
286 if (
m_config[
"topts"].as<string>().compare(
"NotSet") != 0)
288 string inFile =
m_config[
"topts"].as<
string>();
290 if (boost::filesystem::path(inFile).extension() ==
".pts")
296 ptsIO->Import(inFile,
m_f->m_fieldPts);
298 else if (boost::filesystem::path(inFile).extension() ==
".csv")
304 csvIO->Import(inFile,
m_f->m_fieldPts);
308 ASSERTL0(
false,
"unknown topts file type");
311 else if (
m_config[
"line"].as<string>().compare(
"NotSet") != 0)
313 vector<NekDouble> values;
316 "Failed to interpret line string");
318 ASSERTL0(values.size() > 2,
"line string should contain 2*Dim+1 values "
319 "N,x0,y0,z0,x1,y1,z1");
322 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
323 ASSERTL0(values[0] > 1,
"N is not a valid number");
325 int dim = (values.size() - 1) / 2;
326 int npts = values[0];
329 int ptsPerProc = npts / nprocs;
330 int extraPts = (rank < nprocs - 1) ? 0 : npts % nprocs;
331 int locPts = ptsPerProc + extraPts;
332 int start = rank * ptsPerProc;
333 int end = start + locPts;
337 for (
int i = 0; i < dim; ++i)
340 delta[i] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
343 for (
int i = 0, cntLoc = 0; i < npts; ++i)
345 if (i >= start && i < end)
347 for (
int n = 0; n < dim; ++n)
349 pts[n][cntLoc] = values[n + 1] + i * delta[n];
360 m_f->m_fieldPts->SetPointsPerEdge(ppe);
362 else if (
m_config[
"plane"].as<string>().compare(
"NotSet") != 0)
364 vector<NekDouble> values;
367 "Failed to interpret plane string");
370 "plane string should contain 4 Dim+2 values "
371 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
374 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
375 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
377 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
378 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
380 int dim = (values.size() - 2) / 4;
386 int totpts = npts[0] * npts[1];
389 int ptsPerProc = totpts / nprocs;
390 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
391 int locPts = ptsPerProc + extraPts;
392 int start = rank * ptsPerProc;
393 int end = start + locPts;
398 for (
int i = 0; i < dim; ++i)
401 delta1[i] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
403 delta2[i] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
407 for (
int j = 0, cnt = 0, cntLoc = 0; j < npts[1]; ++j)
409 for (
int i = 0; i < npts[0]; ++i, ++cnt)
411 if (cnt >= start && cnt < end)
413 for (
int n = 0; n < dim; ++n)
416 (values[2 + n] + i * delta1[n]) *
418 (values[2 + 3 * dim + n] + i * delta2[n]) *
427 ppe.push_back(npts[0]);
428 ppe.push_back(npts[1]);
432 m_f->m_fieldPts->SetPointsPerEdge(ppe);
434 else if (
m_config[
"box"].as<string>().compare(
"NotSet") != 0)
436 vector<NekDouble> values;
439 "Failed to interpret box string");
441 ASSERTL0(values.size() == 9,
"box string should contain 9 values "
442 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
450 int totpts = npts[0] * npts[1] * npts[2];
456 int ptsPerProc = totpts / nprocs;
457 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
458 int locPts = ptsPerProc + extraPts;
459 int start = rank * ptsPerProc;
460 int end = start + locPts;
462 for (
int i = 0; i < dim; ++i)
465 delta[i] = (values[4 + 2 * i] - values[3 + 2 * i]) / (npts[i] - 1);
468 for (
int k = 0, cnt = 0, cntLoc = 0; k < npts[2]; ++k)
470 for (
int j = 0; j < npts[1]; ++j)
472 for (
int i = 0; i < npts[0]; ++i, ++cnt)
474 if (cnt >= start && cnt < end)
476 pts[0][cntLoc] = values[3] + i * delta[0];
477 pts[1][cntLoc] = values[5] + j * delta[1];
478 pts[2][cntLoc] = values[7] + k * delta[2];
486 ppe.push_back(npts[0]);
487 ppe.push_back(npts[1]);
488 ppe.push_back(npts[2]);
492 m_f->m_fieldPts->SetPointsPerEdge(ppe);
493 vector<NekDouble> boxdim;
494 boxdim.assign(&values[3], &values[3] + 6);
495 m_f->m_fieldPts->SetBoxSize(boxdim);
499 ASSERTL0(
false,
"Missing target points for ProcessInterpPoints.");
504 vector<MultiRegions::ExpListSharedPtr> &field0,
508 boost::ignore_unused(def_value);
510 ASSERTL0(pts->GetNFields() == field0.size(),
"ptField has too few fields");
512 int nfields = field0.size();
515 if (
m_f->m_comm->GetRank() == 0)
521 if (
m_f->m_comm->GetRank() == 0)
526 for (
int f = 0; f < nfields; ++f)
528 for (
int i = 0; i < pts->GetNpoints(); ++i)
530 if (pts->GetPointVal(f, i) > clamp_up)
532 pts->SetPointVal(f, i, clamp_up);
534 else if (pts->GetPointVal(f, i) < clamp_low)
536 pts->SetPointVal(f, i, clamp_low);
545 int dim = pts->GetDim();
546 int nq1 = pts->GetNpoints();
552 vector<NekDouble> values;
554 "Failed to interpret cp string");
556 ASSERTL0(values.size() == 2,
"cp string should contain 2 values "
557 "p0 and q (=1/2 rho u^2)");
560 qinv = 1.0 / values[1];
562 for (
int i = 0; i < pts->GetNFields(); ++i)
564 if (boost::iequals(pts->GetFieldName(i),
"p"))
569 if (boost::iequals(pts->GetFieldName(i),
"u") ||
570 boost::iequals(pts->GetFieldName(i),
"v") ||
571 boost::iequals(pts->GetFieldName(i),
"w"))
581 WARNINGL0(
false,
"Did not find velocity components for Cp0");
586 WARNINGL0(
false,
"Failed to find 'p' field to determine cp0");
592 for (f = 0; f < 2; ++f)
597 for (r = 0; r < nq1; r++)
601 data[0][r] = qinv * (pts->GetPointVal(dim + pfield, r) - p0);
606 for (
int i = 0; i < velid.size(); ++i)
608 q += 0.5 * pts->GetPointVal(dim + velid[i], r) *
609 pts->GetPointVal(dim + velid[i], r);
612 qinv * (pts->GetPointVal(dim + pfield, r) + q - p0);
619 pts->AddField(data[0],
"Cp");
620 m_f->m_variables.push_back(
"Cp");
623 pts->AddField(data[1],
"Cp0");
624 m_f->m_variables.push_back(
"Cp0");
630 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 std::vector< MultiRegions::ExpListSharedPtr > expInField, std::vector< MultiRegions::ExpListSharedPtr > &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)
virtual void Process(po::variables_map &vm)
Write mesh to output file.
virtual ~ProcessInterpPoints()
void CreateFieldPts(po::variables_map &vm)
void PrintProgressbar(const int position, const int goal) const
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 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)
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
The above copyright notice and this permission notice shall be included.
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.