Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | List of all members
Nektar::Utilities::ProcessInterpPoints Class Reference

This processing module interpolates one field to another. More...

#include <ProcessInterpPoints.h>

Inheritance diagram for Nektar::Utilities::ProcessInterpPoints:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessInterpPoints:
Collaboration graph
[legend]

Public Member Functions

 ProcessInterpPoints (FieldSharedPtr f)
 
virtual ~ProcessInterpPoints ()
 
virtual void Process (po::variables_map &vm)
 Write mesh to output file. More...
 
virtual std::string GetModuleName ()
 
- Public Member Functions inherited from Nektar::Utilities::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
 ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
 
void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
void PrintConfig ()
 Print out all configuration options for a module. More...
 
void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
bool GetRequireEquiSpaced (void)
 
void SetRequireEquiSpaced (bool pVal)
 
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 
 Module (MeshSharedPtr p_m)
 
virtual void Process ()=0
 
void RegisterConfig (std::string key, std::string value)
 
void PrintConfig ()
 
void SetDefaults ()
 
MeshSharedPtr GetMesh ()
 
virtual void ProcessVertices ()
 Extract element vertices. More...
 
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual void ProcessElements ()
 Generate element IDs. More...
 
virtual void ProcessComposites ()
 Generate composites. More...
 
virtual void ClearElementLinks ()
 

Static Public Member Functions

static boost::shared_ptr< Modulecreate (FieldSharedPtr f)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Private Member Functions

void InterpolateFieldToPts (vector< MultiRegions::ExpListSharedPtr > &field0, Array< OneD, Array< OneD, NekDouble > > &pts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value, bool isRoot)
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::Utilities::Module
 Module ()
 
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
void PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
MeshSharedPtr m_mesh
 Mesh object. More...
 
std::map< std::string,
ConfigOption
m_config
 List of configuration values. More...
 

Detailed Description

This processing module interpolates one field to another.

Definition at line 49 of file ProcessInterpPoints.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessInterpPoints::ProcessInterpPoints ( FieldSharedPtr  f)

Definition at line 58 of file ProcessInterpPoints.cpp.

References ASSERTL0, and Nektar::Utilities::Module::m_config.

58  : ProcessModule(f)
59 {
60 
61  m_config["fromxml"] =
62  ConfigOption(false, "NotSet",
63  "Xml file from which to interpolate field");
64 
65  ASSERTL0(m_config["fromxml"].as<string>().compare("NotSet") != 0,
66  "Need to specify fromxml=file.xml");
67 
68  m_config["fromfld"] =
69  ConfigOption(false,"NotSet",
70  "Fld file from which to interpolate field");
71 
72  ASSERTL0(m_config["fromfld"].as<string>().compare("NotSet") != 0,
73  "Need to specify fromfld=file.fld ");
74 
75  m_config["clamptolowervalue"] =
76  ConfigOption(false, "-10000000",
77  "Lower bound for interpolation value");
78  m_config["clamptouppervalue"] =
79  ConfigOption(false, "10000000",
80  "Upper bound for interpolation value");
81  m_config["defaultvalue"] =
82  ConfigOption(false, "0",
83  "Default value if point is outside domain");
84  m_config["line"] =
85  ConfigOption(false, "NotSet",
86  "Specify a line of N points using "
87  "line=N,x0,y0,z0,z1,y1,z1");
88  m_config["plane"] =
89  ConfigOption(false, "NotSet",
90  "Specify a plane of N1 x N2 points using "
91  "plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,"
92  "y3,z3");
93  m_config["box"] =
94  ConfigOption(false, "NotSet",
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");
98 
99  m_config["cp"] =
100  ConfigOption(false,"NotSet",
101  "Parameters p0 and q to determine pressure coefficients "
102  "(box only currently)");
103 
104 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
map< string, ConfigOption > m_config
List of configuration values.
Nektar::Utilities::ProcessInterpPoints::~ProcessInterpPoints ( )
virtual

Definition at line 106 of file ProcessInterpPoints.cpp.

107 {
108 }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessInterpPoints::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 53 of file ProcessInterpPoints.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

53  {
55  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual std::string Nektar::Utilities::ProcessInterpPoints::GetModuleName ( )
inlinevirtual

Implements Nektar::Utilities::Module.

Definition at line 64 of file ProcessInterpPoints.h.

65  {
66  return "ProcessInterpPoints";
67  }
void Nektar::Utilities::ProcessInterpPoints::InterpolateFieldToPts ( vector< MultiRegions::ExpListSharedPtr > &  field0,
Array< OneD, Array< OneD, NekDouble > > &  pts,
NekDouble  clamp_low,
NekDouble  clamp_up,
NekDouble  def_value,
bool  isRoot 
)
private

Definition at line 525 of file ProcessInterpPoints.cpp.

References ASSERTL0, Nektar::ParseUtils::GenerateUnOrderedVector(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, and WARNINGL0.

Referenced by Process().

532 {
533  int dim = pts.num_elements();
534 
535  Array<OneD, NekDouble> coords(dim), Lcoords(dim);
536  int nq1 = pts[0].num_elements();
537  int elmtid, offset;
538  int r, f;
539  int intpts = 0;
540  int nfields = field0.size();
541  int pfield = -1;
542  NekDouble p0,qinv;
543  vector<int> velid;
544 
545  if(!boost::iequals(m_config["cp"].as<string>(),"NotSet"))
546  {
547 
548  vector<NekDouble> values;
550  m_config["cp"].as<string>().c_str(),values),
551  "Failed to interpret cp string");
552 
553  ASSERTL0(values.size() == 2,
554  "cp string should contain 2 values "
555  "p0 and q (=1/2 rho u^2)");
556 
557  p0 = values[0];
558  qinv = 1.0/values[1];
559 
560  LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
561 
562  for(int i = 0; i < fPts->GetNFields(); ++i)
563  {
564  if(boost::iequals(fPts->GetFieldName(i),"p"))
565  {
566  pfield = i;
567  }
568 
569  if(boost::iequals(fPts->GetFieldName(i),"u")||
570  boost::iequals(fPts->GetFieldName(i),"v")||
571  boost::iequals(fPts->GetFieldName(i),"w"))
572  {
573  velid.push_back(i);
574  }
575  }
576 
577  if(pfield != -1)
578  {
579  Array< OneD, NekDouble > newPts(m_f->m_fieldPts->GetNpoints());
580  m_f->m_fieldPts->AddField(newPts, "Cp");
581  nfields += 1;
582 
583  if(velid.size())
584  {
585  Array< OneD, NekDouble > newPts(m_f->m_fieldPts->GetNpoints());
586  m_f->m_fieldPts->AddField(newPts, "Cp0");
587  nfields += 1;
588  }
589  else
590  {
591  WARNINGL0(false,"Did not find velocity components for Cp0");
592  }
593  }
594  else
595  {
596  WARNINGL0(false,"Failed to find 'p' field to determine cp0");
597  }
598 
599  }
600 
601  // resize data field
602  m_f->m_data.resize(nfields);
603 
604  for (f = 0; f < nfields; ++f)
605  {
606  m_f->m_data[f].resize(nq1);
607  }
608 
609  for (r = 0; r < nq1; r++)
610  {
611  coords[0] = pts[0][r];
612  coords[1] = pts[1][r];
613  if (dim == 3)
614  {
615  coords[2] = pts[2][r];
616  }
617 
618  // Obtain Element and LocalCoordinate to interpolate
619  elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3);
620 
621  if(elmtid >= 0)
622  {
623  offset = field0[0]->GetPhys_Offset(field0[0]->
624  GetOffset_Elmt_Id(elmtid));
625 
626  for (f = 0; f < field0.size(); ++f)
627  {
628  NekDouble value;
629  value = field0[f]->GetExp(elmtid)->
630  StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
631 
632  if ((boost::math::isnan)(value))
633  {
634  ASSERTL0(false, "new value is not a number");
635  }
636  else
637  {
638  value = (value > clamp_up)? clamp_up :
639  ((value < clamp_low)? clamp_low :
640  value);
641 
642  m_f->m_data[f][r] = value;
643  }
644  }
645  }
646  else
647  {
648  for (f = 0; f < field0.size(); ++f)
649  {
650  m_f->m_data[f][r] = def_value;
651  }
652  }
653 
654  if (intpts%1000 == 0&&isRoot)
655  {
656  cout <<"." << flush;
657  }
658  intpts ++;
659 
660  if(pfield != -1) // calculate cp
661  {
662  m_f->m_data[nfields-2][r] = qinv*(m_f->m_data[pfield][r] - p0);
663 
664  if(velid.size()) // calculate cp0
665  {
666  NekDouble q = 0;
667  for(int i = 0; i < velid.size(); ++i)
668  {
669  q += 0.5*m_f->m_data[velid[i]][r]*m_f->m_data[velid[i]][r];
670  }
671 
672  m_f->m_data[nfields-1][r] = qinv*(m_f->m_data[pfield][r]+q - p0);
673  }
674  }
675 
676  }
677 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
map< string, ConfigOption > m_config
List of configuration values.
FieldSharedPtr m_f
Field object.
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:262
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:194
double NekDouble
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
Definition: ParseUtils.hpp:128
void Nektar::Utilities::ProcessInterpPoints::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 110 of file ProcessInterpPoints.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::ePtsBox, Nektar::LibUtilities::ePtsLine, Nektar::LibUtilities::ePtsPlane, Nektar::ParseUtils::GenerateUnOrderedVector(), InterpolateFieldToPts(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, npts, Nektar::LibUtilities::NullFieldMetaDataMap, Nektar::LibUtilities::NullPtsField, Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Vmax(), and Vmath::Vmin().

111 {
112  if(m_f->m_verbose)
113  {
114  if(m_f->m_comm->GetRank() == 0)
115  {
116  cout << "ProcessInterpPoints: interpolating to points..." << endl;
117  }
118  }
119 
120  int rank = m_f->m_comm->GetRank();
121  int nprocs = m_f->m_comm->GetSize();
122 
123  // Check for command line point specification if no .pts file
124  // specified
125  if(m_f->m_fieldPts == LibUtilities::NullPtsField)
126  {
127  if(m_config["line"].as<string>().compare("NotSet") != 0)
128  {
129  string help = m_config["line"].as<string>();
130  vector<NekDouble> values;
132  m_config["line"].as<string>().c_str(),values),
133  "Failed to interpret line string");
134 
135  ASSERTL0(values.size() > 2,
136  "line string should contain 2 Dim+1 values "
137  "N,x0,y0,z0,x1,y1,z1");
138 
139  double tmp;
140  ASSERTL0(std::modf(values[0], &tmp) == 0.0, "N is not an integer");
141  ASSERTL0(values[0] > 1, "N is not a valid number");
142 
143  int dim = (values.size()-1)/2;
144  int npts = values[0];
145  Array<OneD, Array<OneD, NekDouble> > pts(dim);
146 
147  for(int i = 0; i < dim; ++i)
148  {
149  pts[i] = Array<OneD,NekDouble>(npts);
150  }
151 
152  for(int i = 0; i < npts; ++i)
153  {
154  pts[0][i] = values[1]
155  + i/((NekDouble)(npts-1))*(values[dim+1] - values[1]);
156  if(dim > 1)
157  {
158  pts[1][i] = values[2]
159  + i/((NekDouble)(npts-1))*(values[dim+2] - values[2]);
160 
161  if(dim > 2)
162  {
163  pts[2][i] = values[3]
164  + i/((NekDouble)(npts-1))*(values[dim+3] - values[3]);
165  }
166  }
167  }
168 
169  vector<int> ppe;
170  ppe.push_back(npts);
172  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsLine);
173  m_f->m_fieldPts->SetPointsPerEdge(ppe);
174 
175  }
176  else if(m_config["plane"].as<string>().compare("NotSet") != 0)
177  {
178  string help = m_config["plane"].as<string>();
179  vector<NekDouble> values;
181  m_config["plane"].as<string>().c_str(),values),
182  "Failed to interpret plane string");
183 
184  ASSERTL0(values.size() > 9,
185  "plane string should contain 4 Dim+2 values "
186  "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
187 
188  double tmp;
189  ASSERTL0(std::modf(values[0], &tmp) == 0.0, "N1 is not an integer");
190  ASSERTL0(std::modf(values[1], &tmp) == 0.0, "N2 is not an integer");
191 
192  ASSERTL0(values[0] > 1, "N1 is not a valid number");
193  ASSERTL0(values[1] > 1, "N2 is not a valid number");
194 
195  int dim = (values.size()-2)/4;
196 
197  int npts1 = values[0];
198  int npts2 = values[1];
199 
200  Array<OneD, Array<OneD, NekDouble> > pts(dim);
201 
202  int totpts = npts1*npts2;
203  int nlocpts = totpts/nprocs;
204 
205  if(rank < nprocs-1)
206  {
207  for(int i = 0; i < dim; ++i)
208  {
209  pts[i] = Array<OneD,NekDouble>(nlocpts);
210  }
211 
212  int cnt = 0;
213  int cntloc = 0;
214 
215  for(int j = 0; j < npts2; ++j)
216  {
217  for(int i = 0; i < npts1; ++i)
218  {
219 
220  if((cnt >= rank*nlocpts)&&(cnt < (rank+1)*nlocpts))
221  {
222  pts[0][cntloc] =
223  (values[2] + i/((NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((NekDouble)(npts2-1))) +
224  (values[3*dim+2] + i/((NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((NekDouble)(npts2-1)));
225 
226  pts[1][cntloc] =
227  (values[3] + i/((NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((NekDouble)(npts2-1))) +
228  (values[3*dim+3] + i/((NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((NekDouble)(npts2-1)));
229 
230  if(dim > 2)
231  {
232  pts[2][cntloc] =
233  (values[4] + i/((NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((NekDouble)(npts2-1))) +
234  (values[3*dim+4] + i/((NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((NekDouble)(npts2-1)));
235  }
236  cntloc++;
237  }
238  cnt++;
239  }
240  }
241  }
242  else
243  {
244  totpts = totpts - rank*nlocpts;
245 
246  for(int i = 0; i < dim; ++i)
247  {
248  pts[i] = Array<OneD,NekDouble>(totpts);
249  }
250 
251  int cnt = 0;
252  int cntloc = 0;
253 
254  for(int j = 0; j < npts2; ++j)
255  {
256  for(int i = 0; i < npts1; ++i)
257  {
258 
259  if(cnt >= rank*nlocpts)
260  {
261  pts[0][cntloc] =
262  (values[2] + i/((NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((NekDouble)(npts2-1))) +
263  (values[3*dim+2] + i/((NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((NekDouble)(npts2-1)));
264 
265  pts[1][cntloc] =
266  (values[3] + i/((NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((NekDouble)(npts2-1))) +
267  (values[3*dim+3] + i/((NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((NekDouble)(npts2-1)));
268 
269  if(dim > 2)
270  {
271  pts[2][cntloc] =
272  (values[4] + i/((NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((NekDouble)(npts2-1))) +
273  (values[3*dim+4] + i/((NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((NekDouble)(npts2-1)));
274  }
275  cntloc++;
276  }
277  cnt++;
278  }
279  }
280  }
281 
282  vector<int> ppe;
283  ppe.push_back(npts1);
284  ppe.push_back(npts2);
286  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsPlane);
287  m_f->m_fieldPts->SetPointsPerEdge(ppe);
288 
289  }
290  else if(m_config["box"].as<string>().compare("NotSet") != 0)
291  {
292  string help = m_config["box"].as<string>();
293  vector<NekDouble> values;
295  m_config["box"].as<string>().c_str(),values),
296  "Failed to interpret box string");
297 
298  ASSERTL0(values.size() == 9,
299  "box string should contain 9 values "
300  "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
301 
302  int dim = 3;
303 
304  int npts1 = values[0];
305  int npts2 = values[1];
306  int npts3 = values[2];
307 
308  Array<OneD, Array<OneD, NekDouble> > pts(dim);
309 
310  int totpts = npts1*npts2*npts3;
311  int nlocpts = totpts/nprocs;
312 
313  if(rank < nprocs-1) // for rank 0 to nproc-1
314  {
315  totpts = nlocpts;
316 
317  for(int i = 0; i < dim; ++i)
318  {
319  pts[i] = Array<OneD,NekDouble>(totpts);
320  }
321 
322  int cnt = 0;
323  int cntloc = 0;
324 
325  for(int k = 0; k < npts3; ++k)
326  {
327  for(int j = 0; j < npts2; ++j)
328  {
329  for(int i = 0; i < npts1; ++i)
330  {
331  if((cnt >= rank*nlocpts)&&(cnt < (rank+1)*nlocpts))
332  {
333  pts[0][cntloc] = values[3] +
334  i/((NekDouble)(npts1-1))*(values[4]-values[3]);
335  pts[1][cntloc] = values[5] +
336  j/((NekDouble)(npts2-1))*(values[6]-values[5]);
337  pts[2][cntloc] = values[7] +
338  k/((NekDouble)(npts3-1))*(values[8]-values[7]);
339  cntloc++;
340  }
341  cnt++;
342  }
343  }
344  }
345  }
346  else // give last rank all remaining points
347  {
348  totpts = totpts - rank*nlocpts;
349 
350  for(int i = 0; i < dim; ++i)
351  {
352  pts[i] = Array<OneD,NekDouble>(totpts);
353  }
354 
355  int cnt = 0;
356  int cntloc = 0;
357 
358  for(int k = 0; k < npts3; ++k)
359  {
360  for(int j = 0; j < npts2; ++j)
361  {
362  for(int i = 0; i < npts1; ++i)
363  {
364  if(cnt >= rank*nlocpts)
365  {
366  pts[0][cntloc] = values[3] +
367  i/((NekDouble)(npts1-1))*(values[4]-values[3]);
368  pts[1][cntloc] = values[5] +
369  j/((NekDouble)(npts2-1))*(values[6]-values[5]);
370  pts[2][cntloc] = values[7] +
371  k/((NekDouble)(npts3-1))*(values[8]-values[7]);
372  cntloc++;
373  }
374  cnt++;
375  }
376  }
377  }
378  }
379 
380  vector<int> ppe;
381  ppe.push_back(npts1);
382  ppe.push_back(npts2);
383  ppe.push_back(npts3);
385  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsBox);
386  m_f->m_fieldPts->SetPointsPerEdge(ppe);
387  vector<NekDouble> boxdim;
388  boxdim.assign(&values[3],&values[3]+6);
389  m_f->m_fieldPts->SetBoxSize(boxdim);
390  }
391  }
392 
393 
394  FieldSharedPtr fromField = boost::shared_ptr<Field>(new Field());
395 
396  std::vector<std::string> files;
397  // set up session file for from field
398  files.push_back(m_config["fromxml"].as<string>());
399  fromField->m_session = LibUtilities::SessionReader::
400  CreateInstance(0, 0, files);
401 
402  // Set up range based on min and max of local parallel partition
404 
405  int coordim = m_f->m_fieldPts->GetDim();
406  int npts = m_f->m_fieldPts->GetNpoints();
407  Array<OneD, Array<OneD, NekDouble> > pts;
408  m_f->m_fieldPts->GetPts(pts);
409 
410  rng->m_checkShape = false;
411  rng->m_zmin = -1;
412  rng->m_zmax = 1;
413  rng->m_ymin = -1;
414  rng->m_ymax = 1;
415  switch(coordim)
416  {
417  case 3:
418  rng->m_doZrange = true;
419  rng->m_zmin = Vmath::Vmin(npts, pts[2],1);
420  rng->m_zmax = Vmath::Vmax(npts, pts[2],1);
421  if(rng->m_zmax == rng->m_zmin)
422  {
423  rng->m_zmin -=1;
424  rng->m_zmax +=1;
425  }
426  case 2:
427  rng->m_doYrange = true;
428  rng->m_ymin = Vmath::Vmin(npts, pts[1],1);
429  rng->m_ymax = Vmath::Vmax(npts, pts[1],1);
430  case 1:
431  rng->m_doXrange = true;
432  rng->m_xmin = Vmath::Vmin(npts, pts[0],1);
433  rng->m_xmax = Vmath::Vmax(npts, pts[0],1);
434  break;
435  default:
436  ASSERTL0(false,"too many values specfied in range");
437  }
438 
439  // setup rng parameters.
440  fromField->m_graph = SpatialDomains::MeshGraph::Read(fromField->m_session,rng);
441 
442  // Read in local from field partitions
443  const SpatialDomains::ExpansionMap &expansions = fromField->m_graph->GetExpansions();
444 
445 
446  fromField->m_fld = MemoryManager<LibUtilities::FieldIO>
447  ::AllocateSharedPtr(fromField->m_session->GetComm());
448 
449  Array<OneD,int> ElementGIDs(expansions.size());
450  SpatialDomains::ExpansionMap::const_iterator expIt;
451 
452  int i = 0;
453  for (expIt = expansions.begin(); expIt != expansions.end();
454  ++expIt)
455  {
456  ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
457  }
458 
459  // check to see that we do have some elmement in teh domain since
460  // possibly all points could be outside of the domain
461  ASSERTL0(i > 0, "No elements are set. Are the interpolated points "
462  "wihtin the domain given by the xml files?");
463 
464  string fromfld = m_config["fromfld"].as<string>();
465  fromField->m_fld->Import(fromfld,fromField->m_fielddef,
466  fromField->m_data,
468  ElementGIDs);
469 
470  int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
471 
472  //----------------------------------------------
473  // Set up Expansion information to use mode order from field
474  fromField->m_graph->SetExpansions(fromField->m_fielddef);
475 
476  int nfields = fromField->m_fielddef[0]->m_fields.size();
477 
478  fromField->m_exp.resize(nfields);
479  fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,true);
480 
481  m_f->m_exp.resize(nfields);
482 
483  // declare auxiliary fields.
484  for(i = 1; i < nfields; ++i)
485  {
486  fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
487  }
488 
489  // load field into expansion in fromfield.
490  for(int j = 0; j < nfields; ++j)
491  {
492  for (i = 0; i < fromField->m_fielddef.size(); i++)
493  {
494  fromField->m_exp[j]->ExtractDataToCoeffs(
495  fromField->m_fielddef[i],
496  fromField->m_data[i],
497  fromField->m_fielddef[0]->m_fields[j],
498  fromField->m_exp[j]->UpdateCoeffs());
499  }
500  fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
501  fromField->m_exp[j]->UpdatePhys());
502 
503  Array< OneD, NekDouble > newPts(m_f->m_fieldPts->GetNpoints());
504  m_f->m_fieldPts->AddField(newPts, fromField->m_fielddef[0]->m_fields[j]);
505  }
506 
507  if(rank == 0)
508  {
509  cout << "Interpolating on proc 0 [" << flush;
510  }
511 
512  NekDouble clamp_low = m_config["clamptolowervalue"].as<NekDouble>();
513  NekDouble clamp_up = m_config["clamptouppervalue"].as<NekDouble>();
514  NekDouble def_value = m_config["defaultvalue"].as<NekDouble>();
515 
516  InterpolateFieldToPts(fromField->m_exp, pts,
517  clamp_low, clamp_up, def_value, !rank);
518 
519  if(rank == 0)
520  {
521  cout << "]" << endl;
522  }
523 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:121
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.
Definition: Vmath.cpp:765
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:857
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.
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
static std::string npts
Definition: InputFld.cpp:43
boost::shared_ptr< DomainRange > DomainRangeShPtr
Definition: MeshGraph.h:157
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:263
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
Definition: ParseUtils.hpp:128
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:54
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:174

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessInterpPoints::className
static
Initial value:
=
ModuleKey(eProcessModule, "interppoints"),
"Interpolates a set of points to another, requires fromfld and "
"fromxml to be defined, a line, plane or block of points can be "
"defined")

Definition at line 56 of file ProcessInterpPoints.h.