Nektar++
ProcessInterpPtsToPts.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: ProcessInterpPtsToPts.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Interpolate field to a series of specified points.
32//
33////////////////////////////////////////////////////////////////////////////////
34
35#include <iostream>
36#include <string>
37using namespace std;
38
45
47
48namespace Nektar::FieldUtils
49{
50
53 ModuleKey(eProcessModule, "interpptstopts"),
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 "
57 "defined");
58
60 : ProcessModule(f)
61{
62 m_config["topts"] =
63 ConfigOption(false, "NotSet", "Pts file to which interpolate field");
64 m_config["line"] = ConfigOption(false, "NotSet",
65 "Specify a line of N points using "
66 "line=N,x0,y0,z0,z1,y1,z1");
67 m_config["plane"] =
68 ConfigOption(false, "NotSet",
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");
71 m_config["box"] =
72 ConfigOption(false, "NotSet",
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");
76
77 m_config["clamptolowervalue"] =
78 ConfigOption(false, "-10000000", "Lower bound for interpolation value");
79 m_config["clamptouppervalue"] =
80 ConfigOption(false, "10000000", "Upper bound for interpolation value");
81 m_config["defaultvalue"] =
82 ConfigOption(false, "0", "Default value if point is outside domain");
83
84 m_config["cp"] =
85 ConfigOption(false, "NotSet",
86 "Parameters p0 and q to determine pressure coefficients");
87}
88
90{
91}
92
93void ProcessInterpPtsToPts::v_Process(po::variables_map &vm)
94{
96 "Should have a PtsField for ProcessInterpPtsToPts.");
97 ASSERTL0(m_f->m_comm->GetSpaceComm()->GetSize() == 1,
98 "ProcessInterpPtsToPts not implemented in parallel.");
99
100 // Move m_f->m_fieldPts
101 LibUtilities::PtsFieldSharedPtr oldPts = m_f->m_fieldPts;
102 m_f->m_fieldPts = LibUtilities::NullPtsField;
103
104 // Create new fieldPts
105 CreateFieldPts(vm);
106
107 int nfields = m_f->m_variables.size();
108 for (int j = 0; j < nfields; ++j)
109 {
110 Array<OneD, NekDouble> newPts(m_f->m_fieldPts->GetNpoints());
111 m_f->m_fieldPts->AddField(newPts, m_f->m_variables[j]);
112 }
113
114 NekDouble clamp_low = m_config["clamptolowervalue"].as<NekDouble>();
115 NekDouble clamp_up = m_config["clamptouppervalue"].as<NekDouble>();
116 NekDouble def_value = m_config["defaultvalue"].as<NekDouble>();
117
118 InterpolatePtsToPts(oldPts, m_f->m_fieldPts, clamp_low, clamp_up,
119 def_value);
120
121 if (!boost::iequals(m_config["cp"].as<string>(), "NotSet"))
122 {
123 calcCp0();
124 }
125}
126
128 [[maybe_unused]] po::variables_map &vm)
129{
130 int rank = m_f->m_comm->GetSpaceComm()->GetRank();
131 int nprocs = m_f->m_comm->GetSpaceComm()->GetSize();
132 // Check for command line point specification
133 if (m_config["topts"].as<string>().compare("NotSet") != 0)
134 {
135 string inFile = m_config["topts"].as<string>();
136
137 if (fs::path(inFile).extension() == ".pts")
138 {
141 m_f->m_comm);
142
143 ptsIO->Import(inFile, m_f->m_fieldPts);
144 }
145 else if (fs::path(inFile).extension() == ".csv")
146 {
149 m_f->m_comm);
150
151 csvIO->Import(inFile, m_f->m_fieldPts);
152 }
153 else
154 {
155 ASSERTL0(false, "unknown topts file type");
156 }
157 }
158 else if (m_config["line"].as<string>().compare("NotSet") != 0)
159 {
160 vector<NekDouble> values;
161 ASSERTL0(
162 ParseUtils::GenerateVector(m_config["line"].as<string>(), values),
163 "Failed to interpret line string");
164
165 ASSERTL0(values.size() > 2, "line string should contain 2*Dim+1 values "
166 "N,x0,y0,z0,x1,y1,z1");
167
168 double tmp;
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");
171
172 int dim = (values.size() - 1) / 2;
173 int npts = values[0];
174
175 // Information for partitioning
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;
181
183 Array<OneD, NekDouble> delta(dim);
184 for (int i = 0; i < dim; ++i)
185 {
186 pts[i] = Array<OneD, NekDouble>(locPts);
187 delta[i] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
188 }
189
190 for (int i = 0, cntLoc = 0; i < npts; ++i)
191 {
192 if (i >= start && i < end)
193 {
194 for (int n = 0; n < dim; ++n)
195 {
196 pts[n][cntLoc] = values[n + 1] + i * delta[n];
197 }
198 ++cntLoc;
199 }
200 }
201
202 vector<size_t> ppe;
203 ppe.push_back(npts);
204 m_f->m_fieldPts =
206 m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsLine);
207 m_f->m_fieldPts->SetPointsPerEdge(ppe);
208 }
209 else if (m_config["plane"].as<string>().compare("NotSet") != 0)
210 {
211 vector<NekDouble> values;
212 ASSERTL0(
213 ParseUtils::GenerateVector(m_config["plane"].as<string>(), values),
214 "Failed to interpret plane string");
215
216 ASSERTL0(values.size() > 9,
217 "plane string should contain 4 Dim+2 values "
218 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
219
220 double tmp;
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");
223
224 ASSERTL0(values[0] > 1, "N1 is not a valid number");
225 ASSERTL0(values[1] > 1, "N2 is not a valid number");
226
227 int dim = (values.size() - 2) / 4;
228
229 Array<OneD, int> npts(2);
230 npts[0] = values[0];
231 npts[1] = values[1];
232
233 int totpts = npts[0] * npts[1];
234
235 // Information for partitioning
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;
241
243 Array<OneD, NekDouble> delta1(dim);
244 Array<OneD, NekDouble> delta2(dim);
245 for (int i = 0; i < dim; ++i)
246 {
247 pts[i] = Array<OneD, NekDouble>(locPts);
248 delta1[i] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
249 (npts[0] - 1);
250 delta2[i] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
251 (npts[0] - 1);
252 }
253
254 for (int j = 0, cnt = 0, cntLoc = 0; j < npts[1]; ++j)
255 {
256 for (int i = 0; i < npts[0]; ++i, ++cnt)
257 {
258 if (cnt >= start && cnt < end)
259 {
260 for (int n = 0; n < dim; ++n)
261 {
262 pts[n][cntLoc] =
263 (values[2 + n] + i * delta1[n]) *
264 (1.0 - j / ((NekDouble)(npts[1] - 1))) +
265 (values[2 + 3 * dim + n] + i * delta2[n]) *
266 (j / ((NekDouble)(npts[1] - 1)));
267 }
268 ++cntLoc;
269 }
270 }
271 }
272
273 vector<size_t> ppe;
274 ppe.push_back(npts[0]);
275 ppe.push_back(npts[1]);
276 m_f->m_fieldPts =
278 m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsPlane);
279 m_f->m_fieldPts->SetPointsPerEdge(ppe);
280 }
281 else if (m_config["box"].as<string>().compare("NotSet") != 0)
282 {
283 vector<NekDouble> values;
284 ASSERTL0(
285 ParseUtils::GenerateVector(m_config["box"].as<string>(), values),
286 "Failed to interpret box string");
287
288 ASSERTL0(values.size() == 9, "box string should contain 9 values "
289 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
290
291 int dim = 3;
292 Array<OneD, int> npts(3);
293 npts[0] = values[0];
294 npts[1] = values[1];
295 npts[2] = values[2];
296
297 int totpts = npts[0] * npts[1] * npts[2];
298
300 Array<OneD, NekDouble> delta(dim);
301
302 // Information for partitioning
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;
308
309 for (int i = 0; i < dim; ++i)
310 {
311 pts[i] = Array<OneD, NekDouble>(locPts);
312 delta[i] = (values[4 + 2 * i] - values[3 + 2 * i]) / (npts[i] - 1);
313 }
314
315 for (int k = 0, cnt = 0, cntLoc = 0; k < npts[2]; ++k)
316 {
317 for (int j = 0; j < npts[1]; ++j)
318 {
319 for (int i = 0; i < npts[0]; ++i, ++cnt)
320 {
321 if (cnt >= start && cnt < end)
322 {
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];
326 ++cntLoc;
327 }
328 }
329 }
330 }
331
332 vector<size_t> ppe;
333 ppe.push_back(npts[0]);
334 ppe.push_back(npts[1]);
335 ppe.push_back(npts[2]);
336 m_f->m_fieldPts =
338 m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsBox);
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);
343 }
344 else
345 {
346 ASSERTL0(false,
347 "ProcessInterpPtsToPts requires line, plane or box option.");
348 }
349}
350
354 NekDouble clamp_up, [[maybe_unused]] NekDouble def_value)
355{
356 ASSERTL0(toPts->GetNFields() >= fromPts->GetNFields(),
357 "ptField has too few fields");
358
359 int nfields = fromPts->GetNFields();
360
362 if (m_f->m_comm->GetRank() == 0)
363 {
365 this);
366 }
367 interp.Interpolate(fromPts, toPts);
368 if (m_f->m_comm->GetRank() == 0)
369 {
370 cout << endl;
371 }
372
373 for (int f = 0; f < nfields; ++f)
374 {
375 for (int i = 0; i < toPts->GetNpoints(); ++i)
376 {
377 if (toPts->GetPointVal(f, i) > clamp_up)
378 {
379 toPts->SetPointVal(f, i, clamp_up);
380 }
381 else if (toPts->GetPointVal(f, i) < clamp_low)
382 {
383 toPts->SetPointVal(f, i, clamp_low);
384 }
385 }
386 }
387}
388
390{
391 LibUtilities::PtsFieldSharedPtr pts = m_f->m_fieldPts;
392 int dim = pts->GetDim();
393 int nq1 = pts->GetNpoints();
394 int r, f;
395 int pfield = -1;
396 NekDouble p0, qinv;
397 vector<int> velid;
398
399 vector<NekDouble> values;
400 ASSERTL0(ParseUtils::GenerateVector(m_config["cp"].as<string>(), values),
401 "Failed to interpret cp string");
402
403 ASSERTL0(values.size() == 2, "cp string should contain 2 values "
404 "p0 and q (=1/2 rho u^2)");
405
406 p0 = values[0];
407 qinv = 1.0 / values[1];
408
409 for (int i = 0; i < pts->GetNFields(); ++i)
410 {
411 if (boost::iequals(pts->GetFieldName(i), "p"))
412 {
413 pfield = i;
414 }
415
416 if (boost::iequals(pts->GetFieldName(i), "u") ||
417 boost::iequals(pts->GetFieldName(i), "v") ||
418 boost::iequals(pts->GetFieldName(i), "w"))
419 {
420 velid.push_back(i);
421 }
422 }
423
424 if (pfield != -1)
425 {
426 if (!velid.size())
427 {
428 WARNINGL0(false, "Did not find velocity components for Cp0");
429 }
430 }
431 else
432 {
433 WARNINGL0(false, "Failed to find 'p' field to determine cp0");
434 }
435
436 // Allocate data storage
438
439 for (f = 0; f < 2; ++f)
440 {
441 data[f] = Array<OneD, NekDouble>(nq1, 0.0);
442 }
443
444 for (r = 0; r < nq1; r++)
445 {
446 if (pfield != -1) // calculate cp
447 {
448 data[0][r] = qinv * (pts->GetPointVal(dim + pfield, r) - p0);
449
450 if (velid.size()) // calculate cp0
451 {
452 NekDouble q = 0;
453 for (int i = 0; i < velid.size(); ++i)
454 {
455 q += 0.5 * pts->GetPointVal(dim + velid[i], r) *
456 pts->GetPointVal(dim + velid[i], r);
457 }
458 data[1][r] =
459 qinv * (pts->GetPointVal(dim + pfield, r) + q - p0);
460 }
461 }
462 }
463
464 if (pfield != -1)
465 {
466 pts->AddField(data[0], "Cp");
467 m_f->m_variables.push_back("Cp");
468 if (velid.size())
469 {
470 pts->AddField(data[1], "Cp0");
471 m_f->m_variables.push_back("Cp0");
472 }
473 }
474}
475
477 const int goal) const
478{
479 LibUtilities::PrintProgressbar(position, goal, "Interpolating");
480}
481} // namespace Nektar::FieldUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
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.
Definition: Module.h:239
std::map< std::string, ConfigOption > m_config
List of configuration values.
Definition: Module.h:272
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.
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.
Definition: Module.h:301
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.
Definition: ParseUtils.cpp:130
std::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:1026
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:180
ModuleFactory & GetModuleFactory()
Definition: Module.cpp:47
int PrintProgressbar(const int position, const int goal, const std::string message, int lastprogress=-1)
Prints a progressbar.
Definition: Progressbar.hpp:65
std::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:184
std::shared_ptr< CsvIO > CsvIOSharedPtr
Definition: CsvIO.h:74
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:185
std::shared_ptr< PtsIO > PtsIOSharedPtr
Definition: PtsIO.h:90
std::vector< double > q(NPUPPER *NPUPPER)
double NekDouble
STL namespace.
Represents a command-line configuration option.
Definition: Module.h:129