55 "Interpolates a set of points to another, requires fromfld and "
56 "fromxml to be defined, a line, plane or block of points can be "
63 ConfigOption(
false,
"NotSet",
"Pts file to which interpolate field");
65 "Specify a line of N points using "
66 "line=N,x0,y0,z0,z1,y1,z1");
69 "Specify a plane of N1 x N2 points using "
70 "plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,y3,z3");
73 "Specify a rectangular box of N1 x N2 x N3 points "
74 "using a box of points limited by box="
75 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
78 ConfigOption(
false,
"-10000000",
"Lower bound for interpolation value");
80 ConfigOption(
false,
"10000000",
"Upper bound for interpolation value");
82 ConfigOption(
false,
"0",
"Default value if point is outside domain");
86 "Parameters p0 and q to determine pressure coefficients");
96 "Should have a PtsField for ProcessInterpPtsToPts.");
97 ASSERTL0(
m_f->m_comm->GetSpaceComm()->GetSize() == 1,
98 "ProcessInterpPtsToPts not implemented in parallel.");
107 int nfields =
m_f->m_variables.size();
108 for (
int j = 0; j < nfields; ++j)
111 m_f->m_fieldPts->AddField(newPts,
m_f->m_variables[j]);
121 if (!boost::iequals(
m_config[
"cp"].as<string>(),
"NotSet"))
128 [[maybe_unused]] po::variables_map &vm)
130 int rank =
m_f->m_comm->GetSpaceComm()->GetRank();
131 int nprocs =
m_f->m_comm->GetSpaceComm()->GetSize();
133 if (
m_config[
"topts"].as<string>().compare(
"NotSet") != 0)
135 string inFile =
m_config[
"topts"].as<
string>();
137 if (fs::path(inFile).extension() ==
".pts")
143 ptsIO->Import(inFile,
m_f->m_fieldPts);
145 else if (fs::path(inFile).extension() ==
".csv")
151 csvIO->Import(inFile,
m_f->m_fieldPts);
155 ASSERTL0(
false,
"unknown topts file type");
158 else if (
m_config[
"line"].as<string>().compare(
"NotSet") != 0)
160 vector<NekDouble> values;
163 "Failed to interpret line string");
165 ASSERTL0(values.size() > 2,
"line string should contain 2*Dim+1 values "
166 "N,x0,y0,z0,x1,y1,z1");
169 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
170 ASSERTL0(values[0] > 1,
"N is not a valid number");
172 int dim = (values.size() - 1) / 2;
173 int npts = values[0];
176 int ptsPerProc = npts / nprocs;
177 int extraPts = (rank < nprocs - 1) ? 0 : npts % nprocs;
178 int locPts = ptsPerProc + extraPts;
179 int start = rank * ptsPerProc;
180 int end = start + locPts;
184 for (
int i = 0; i < dim; ++i)
187 delta[i] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
190 for (
int i = 0, cntLoc = 0; i < npts; ++i)
192 if (i >= start && i < end)
194 for (
int n = 0; n < dim; ++n)
196 pts[n][cntLoc] = values[n + 1] + i * delta[n];
207 m_f->m_fieldPts->SetPointsPerEdge(ppe);
209 else if (
m_config[
"plane"].as<string>().compare(
"NotSet") != 0)
211 vector<NekDouble> values;
214 "Failed to interpret plane string");
217 "plane string should contain 4 Dim+2 values "
218 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
221 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
222 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
224 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
225 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
227 int dim = (values.size() - 2) / 4;
233 int totpts = npts[0] * npts[1];
236 int ptsPerProc = totpts / nprocs;
237 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
238 int locPts = ptsPerProc + extraPts;
239 int start = rank * ptsPerProc;
240 int end = start + locPts;
245 for (
int i = 0; i < dim; ++i)
248 delta1[i] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
250 delta2[i] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
254 for (
int j = 0, cnt = 0, cntLoc = 0; j < npts[1]; ++j)
256 for (
int i = 0; i < npts[0]; ++i, ++cnt)
258 if (cnt >= start && cnt < end)
260 for (
int n = 0; n < dim; ++n)
263 (values[2 + n] + i * delta1[n]) *
265 (values[2 + 3 * dim + n] + i * delta2[n]) *
274 ppe.push_back(npts[0]);
275 ppe.push_back(npts[1]);
279 m_f->m_fieldPts->SetPointsPerEdge(ppe);
281 else if (
m_config[
"box"].as<string>().compare(
"NotSet") != 0)
283 vector<NekDouble> values;
286 "Failed to interpret box string");
288 ASSERTL0(values.size() == 9,
"box string should contain 9 values "
289 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
297 int totpts = npts[0] * npts[1] * npts[2];
303 int ptsPerProc = totpts / nprocs;
304 int extraPts = (rank < nprocs - 1) ? 0 : totpts % nprocs;
305 int locPts = ptsPerProc + extraPts;
306 int start = rank * ptsPerProc;
307 int end = start + locPts;
309 for (
int i = 0; i < dim; ++i)
312 delta[i] = (values[4 + 2 * i] - values[3 + 2 * i]) / (npts[i] - 1);
315 for (
int k = 0, cnt = 0, cntLoc = 0; k < npts[2]; ++k)
317 for (
int j = 0; j < npts[1]; ++j)
319 for (
int i = 0; i < npts[0]; ++i, ++cnt)
321 if (cnt >= start && cnt < end)
323 pts[0][cntLoc] = values[3] + i * delta[0];
324 pts[1][cntLoc] = values[5] + j * delta[1];
325 pts[2][cntLoc] = values[7] + k * delta[2];
333 ppe.push_back(npts[0]);
334 ppe.push_back(npts[1]);
335 ppe.push_back(npts[2]);
339 m_f->m_fieldPts->SetPointsPerEdge(ppe);
340 vector<NekDouble> boxdim;
341 boxdim.assign(&values[3], &values[3] + 6);
342 m_f->m_fieldPts->SetBoxSize(boxdim);
347 "ProcessInterpPtsToPts requires line, plane or box option.");
356 ASSERTL0(toPts->GetNFields() >= fromPts->GetNFields(),
357 "ptField has too few fields");
359 int nfields = fromPts->GetNFields();
362 if (
m_f->m_comm->GetRank() == 0)
368 if (
m_f->m_comm->GetRank() == 0)
373 for (
int f = 0; f < nfields; ++f)
375 for (
int i = 0; i < toPts->GetNpoints(); ++i)
377 if (toPts->GetPointVal(f, i) > clamp_up)
379 toPts->SetPointVal(f, i, clamp_up);
381 else if (toPts->GetPointVal(f, i) < clamp_low)
383 toPts->SetPointVal(f, i, clamp_low);
392 int dim = pts->GetDim();
393 int nq1 = pts->GetNpoints();
399 vector<NekDouble> values;
401 "Failed to interpret cp string");
403 ASSERTL0(values.size() == 2,
"cp string should contain 2 values "
404 "p0 and q (=1/2 rho u^2)");
407 qinv = 1.0 / values[1];
409 for (
int i = 0; i < pts->GetNFields(); ++i)
411 if (boost::iequals(pts->GetFieldName(i),
"p"))
416 if (boost::iequals(pts->GetFieldName(i),
"u") ||
417 boost::iequals(pts->GetFieldName(i),
"v") ||
418 boost::iequals(pts->GetFieldName(i),
"w"))
428 WARNINGL0(
false,
"Did not find velocity components for Cp0");
433 WARNINGL0(
false,
"Failed to find 'p' field to determine cp0");
439 for (f = 0; f < 2; ++f)
444 for (r = 0; r < nq1; r++)
448 data[0][r] = qinv * (pts->GetPointVal(dim + pfield, r) - p0);
453 for (
int i = 0; i < velid.size(); ++i)
455 q += 0.5 * pts->GetPointVal(dim + velid[i], r) *
456 pts->GetPointVal(dim + velid[i], r);
459 qinv * (pts->GetPointVal(dim + pfield, r) +
q - p0);
466 pts->AddField(data[0],
"Cp");
467 m_f->m_variables.push_back(
"Cp");
470 pts->AddField(data[1],
"Cp0");
471 m_f->m_variables.push_back(
"Cp0");
477 const int goal)
const
#define ASSERTL0(condition, msg)
#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.
static ModuleKey className
void CreateFieldPts(po::variables_map &vm)
void InterpolatePtsToPts(LibUtilities::PtsFieldSharedPtr &fromPts, LibUtilities::PtsFieldSharedPtr &toPts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value)
static std::shared_ptr< Module > create(FieldSharedPtr f)
Creates an instance of this class.
~ProcessInterpPtsToPts() override
ProcessInterpPtsToPts(FieldSharedPtr f)
void PrintProgressbar(const int position, const int goal) const
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 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.
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.
std::shared_ptr< PtsField > PtsFieldSharedPtr
std::shared_ptr< CsvIO > CsvIOSharedPtr
static PtsFieldSharedPtr NullPtsField
std::shared_ptr< PtsIO > PtsIOSharedPtr
std::vector< double > q(NPUPPER *NPUPPER)
Represents a command-line configuration option.