43 #include <boost/math/special_functions/fpclassify.hpp>
49 ModuleKey ProcessInterpPoints::className =
52 ProcessInterpPoints::create,
53 "Interpolates a set of points to another, requires fromfld and "
54 "fromxml to be defined, a line, plane or block of points can be "
63 "Xml file from which to interpolate field");
66 "Need to specify fromxml=file.xml");
70 "Fld file from which to interpolate field");
73 "Need to specify fromfld=file.fld ");
77 "Lower bound for interpolation value");
80 "Upper bound for interpolation value");
83 "Default value if point is outside domain");
86 "Specify a line of N points using "
87 "line=N,x0,y0,z0,z1,y1,z1");
90 "Specify a plane of N1 x N2 points using "
91 "plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,"
95 "Specify a rectangular box of N1 x N2 x N3 points "
96 "using a box of points limited by box="
97 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
101 "Parameters p0 and q to determine pressure coefficients "
102 "(box only currently)");
113 int rank =
m_f->m_comm->GetRank();
114 int nprocs =
m_f->m_comm->GetSize();
116 if((
m_f->m_verbose)&&(rank == 0))
118 cout <<
"Processing point interpolation" << endl;
126 if(
m_config[
"line"].as<string>().compare(
"NotSet") != 0)
128 string help =
m_config[
"line"].as<
string>();
129 vector<NekDouble> values;
131 m_config[
"line"].as<string>().c_str(),values),
132 "Failed to interpret line string");
135 "line string should contain 2 Dim+1 values "
136 "N,x0,y0,z0,x1,y1,z1");
139 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
140 ASSERTL0(values[0] > 1,
"N is not a valid number");
142 int dim = (values.size()-1)/2;
143 int npts = values[0];
146 for(
int i = 0; i < dim; ++i)
151 for(
int i = 0; i <
npts; ++i)
153 pts[0][i] = values[1]
154 + i/((
NekDouble)(npts-1))*(values[dim+1] - values[1]);
157 pts[1][i] = values[2]
158 + i/((
NekDouble)(npts-1))*(values[dim+2] - values[2]);
162 pts[2][i] = values[3]
163 + i/((
NekDouble)(npts-1))*(values[dim+3] - values[3]);
172 m_f->m_fieldPts->SetPointsPerEdge(ppe);
175 else if(
m_config[
"plane"].as<string>().compare(
"NotSet") != 0)
177 string help =
m_config[
"plane"].as<
string>();
178 vector<NekDouble> values;
180 m_config[
"plane"].as<string>().c_str(),values),
181 "Failed to interpret plane string");
184 "plane string should contain 4 Dim+2 values "
185 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
188 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
189 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
191 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
192 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
194 int dim = (values.size()-2)/4;
196 int npts1 = values[0];
197 int npts2 = values[1];
201 int totpts = npts1*npts2;
202 int nlocpts = totpts/nprocs;
206 for(
int i = 0; i < dim; ++i)
214 for(
int j = 0; j < npts2; ++j)
216 for(
int i = 0; i < npts1; ++i)
219 if((cnt >= rank*nlocpts)&&(cnt < (rank+1)*nlocpts))
222 (values[2] + i/((
NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((
NekDouble)(npts2-1))) +
223 (values[3*dim+2] + i/((
NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((
NekDouble)(npts2-1)));
226 (values[3] + i/((
NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((
NekDouble)(npts2-1))) +
227 (values[3*dim+3] + i/((
NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((
NekDouble)(npts2-1)));
232 (values[4] + i/((
NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((
NekDouble)(npts2-1))) +
233 (values[3*dim+4] + i/((
NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((
NekDouble)(npts2-1)));
243 totpts = totpts - rank*nlocpts;
245 for(
int i = 0; i < dim; ++i)
253 for(
int j = 0; j < npts2; ++j)
255 for(
int i = 0; i < npts1; ++i)
258 if(cnt >= rank*nlocpts)
261 (values[2] + i/((
NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((
NekDouble)(npts2-1))) +
262 (values[3*dim+2] + i/((
NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((
NekDouble)(npts2-1)));
265 (values[3] + i/((
NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((
NekDouble)(npts2-1))) +
266 (values[3*dim+3] + i/((
NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((
NekDouble)(npts2-1)));
271 (values[4] + i/((
NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((
NekDouble)(npts2-1))) +
272 (values[3*dim+4] + i/((
NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((
NekDouble)(npts2-1)));
282 ppe.push_back(npts1);
283 ppe.push_back(npts2);
286 m_f->m_fieldPts->SetPointsPerEdge(ppe);
289 else if(
m_config[
"box"].as<string>().compare(
"NotSet") != 0)
291 string help =
m_config[
"box"].as<
string>();
292 vector<NekDouble> values;
294 m_config[
"box"].as<string>().c_str(),values),
295 "Failed to interpret box string");
298 "box string should contain 9 values "
299 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
303 int npts1 = values[0];
304 int npts2 = values[1];
305 int npts3 = values[2];
309 int totpts = npts1*npts2*npts3;
310 int nlocpts = totpts/nprocs;
316 for(
int i = 0; i < dim; ++i)
324 for(
int k = 0; k < npts3; ++k)
326 for(
int j = 0; j < npts2; ++j)
328 for(
int i = 0; i < npts1; ++i)
330 if((cnt >= rank*nlocpts)&&(cnt < (rank+1)*nlocpts))
332 pts[0][cntloc] = values[3] +
333 i/((
NekDouble)(npts1-1))*(values[4]-values[3]);
334 pts[1][cntloc] = values[5] +
335 j/((
NekDouble)(npts2-1))*(values[6]-values[5]);
336 pts[2][cntloc] = values[7] +
337 k/((
NekDouble)(npts3-1))*(values[8]-values[7]);
347 totpts = totpts - rank*nlocpts;
349 for(
int i = 0; i < dim; ++i)
357 for(
int k = 0; k < npts3; ++k)
359 for(
int j = 0; j < npts2; ++j)
361 for(
int i = 0; i < npts1; ++i)
363 if(cnt >= rank*nlocpts)
365 pts[0][cntloc] = values[3] +
366 i/((
NekDouble)(npts1-1))*(values[4]-values[3]);
367 pts[1][cntloc] = values[5] +
368 j/((
NekDouble)(npts2-1))*(values[6]-values[5]);
369 pts[2][cntloc] = values[7] +
370 k/((
NekDouble)(npts3-1))*(values[8]-values[7]);
380 ppe.push_back(npts1);
381 ppe.push_back(npts2);
382 ppe.push_back(npts3);
385 m_f->m_fieldPts->SetPointsPerEdge(ppe);
386 vector<NekDouble> boxdim;
387 boxdim.assign(&values[3],&values[3]+6);
388 m_f->m_fieldPts->SetBoxSize(boxdim);
395 std::vector<std::string> files;
397 files.push_back(
m_config[
"fromxml"].as<string>());
404 int coordim =
m_f->m_fieldPts->GetDim();
405 int npts =
m_f->m_fieldPts->GetNpoints();
407 m_f->m_fieldPts->GetPts(pts);
409 rng->m_checkShape =
false;
417 rng->m_doZrange =
true;
420 if(rng->m_zmax == rng->m_zmin)
426 rng->m_doYrange =
true;
430 rng->m_doXrange =
true;
435 ASSERTL0(
false,
"too many values specfied in range");
449 SpatialDomains::ExpansionMap::const_iterator expIt;
452 for (expIt = expansions.begin(); expIt != expansions.end();
455 ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
460 ASSERTL0(i > 0,
"No elements are set. Are the interpolated points "
461 "wihtin the domain given by the xml files?");
463 string fromfld =
m_config[
"fromfld"].as<
string>();
464 fromField->m_fld->Import(fromfld,fromField->m_fielddef,
469 int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
473 fromField->m_graph->SetExpansions(fromField->m_fielddef);
475 int nfields = fromField->m_fielddef[0]->m_fields.size();
477 fromField->m_exp.resize(nfields);
478 fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,
true);
480 m_f->m_exp.resize(nfields);
483 for(i = 1; i < nfields; ++i)
485 fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
489 for(
int j = 0; j < nfields; ++j)
491 for (i = 0; i < fromField->m_fielddef.size(); i++)
493 fromField->m_exp[j]->ExtractDataToCoeffs(
494 fromField->m_fielddef[i],
495 fromField->m_data[i],
496 fromField->m_fielddef[0]->m_fields[j],
497 fromField->m_exp[j]->UpdateCoeffs());
499 fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
500 fromField->m_exp[j]->UpdatePhys());
503 m_f->m_fieldPts->AddField(newPts, fromField->m_fielddef[0]->m_fields[j]);
508 cout <<
"Interpolating on proc 0 [" << flush;
516 clamp_low, clamp_up, def_value, !rank);
526 vector<MultiRegions::ExpListSharedPtr> &field0,
533 int dim = pts.num_elements();
536 int nq1 = pts[0].num_elements();
540 int nfields = field0.size();
545 if(!boost::iequals(
m_config[
"cp"].as<string>(),
"NotSet"))
548 vector<NekDouble> values;
550 m_config[
"cp"].as<string>().c_str(),values),
551 "Failed to interpret cp string");
554 "cp string should contain 2 values "
555 "p0 and q (=1/2 rho u^2)");
558 qinv = 1.0/values[1];
562 for(
int i = 0; i < fPts->GetNFields(); ++i)
564 if(boost::iequals(fPts->GetFieldName(i),
"p"))
569 if(boost::iequals(fPts->GetFieldName(i),
"u")||
570 boost::iequals(fPts->GetFieldName(i),
"v")||
571 boost::iequals(fPts->GetFieldName(i),
"w"))
580 m_f->m_fieldPts->AddField(newPts,
"Cp");
586 m_f->m_fieldPts->AddField(newPts,
"Cp0");
591 WARNINGL0(
false,
"Did not find velocity components for Cp0");
596 WARNINGL0(
false,
"Failed to find 'p' field to determine cp0");
602 m_f->m_data.resize(nfields);
604 for (f = 0; f < nfields; ++f)
606 m_f->m_data[f].resize(nq1);
609 for (r = 0; r < nq1; r++)
611 coords[0] = pts[0][r];
612 coords[1] = pts[1][r];
615 coords[2] = pts[2][r];
619 elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3);
623 offset = field0[0]->GetPhys_Offset(field0[0]->
624 GetOffset_Elmt_Id(elmtid));
626 for (f = 0; f < field0.size(); ++f)
629 value = field0[f]->GetExp(elmtid)->
630 StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
632 if ((boost::math::isnan)(value))
634 ASSERTL0(
false,
"new value is not a number");
638 value = (value > clamp_up)? clamp_up :
639 ((value < clamp_low)? clamp_low :
642 m_f->m_data[f][r] = value;
648 for (f = 0; f < field0.size(); ++f)
650 m_f->m_data[f][r] = def_value;
654 if (intpts%1000 == 0&&isRoot)
662 m_f->m_data[nfields-2][r] = qinv*(
m_f->m_data[pfield][r] - p0);
667 for(
int i = 0; i < velid.size(); ++i)
669 q += 0.5*
m_f->m_data[velid[i]][r]*
m_f->m_data[velid[i]][r];
672 m_f->m_data[nfields-1][r] = qinv*(
m_f->m_data[pfield][r]+q - p0);
#define ASSERTL0(condition, msg)
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
void InterpolateFieldToPts(vector< MultiRegions::ExpListSharedPtr > &field0, Array< OneD, Array< OneD, NekDouble > > &pts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value, bool isRoot)
map< string, ConfigOption > m_config
List of configuration values.
FieldSharedPtr m_f
Field object.
boost::shared_ptr< PtsField > PtsFieldSharedPtr
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
#define WARNINGL0(condition, msg)
boost::shared_ptr< DomainRange > DomainRangeShPtr
boost::shared_ptr< Field > FieldSharedPtr
Represents a command-line configuration option.
static PtsFieldSharedPtr NullPtsField
virtual ~ProcessInterpPoints()
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
static FieldMetaDataMap NullFieldMetaDataMap
std::map< int, ExpansionShPtr > ExpansionMap
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.