Nektar++
ProcessInterpPoints.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessInterpPoints.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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Interpolate field to a series of specified points.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 #include <string>
36 #include <iostream>
37 using namespace std;
38 
39 #include "ProcessInterpPoints.h"
40 
43 #include <boost/math/special_functions/fpclassify.hpp>
44 namespace Nektar
45 {
46 namespace Utilities
47 {
48 
49 ModuleKey ProcessInterpPoints::className =
51  ModuleKey(eProcessModule, "interppoints"),
52  ProcessInterpPoints::create,
53  "Interpolates a set of points to another, requires fromfld and "
54  "fromxml to be defined, a line or plane of points can be "
55  "defined");
56 
57 
58 ProcessInterpPoints::ProcessInterpPoints(FieldSharedPtr f) : 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 
89  m_config["plane"] =
90  ConfigOption(false, "NotSet",
91  "Specify a plane of N1 x N1 points using "
92  "plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,"
93  "y3,z3");
94 }
95 
97 {
98 }
99 
100 void ProcessInterpPoints::Process(po::variables_map &vm)
101 {
102 
103  if(m_f->m_verbose)
104  {
105  cout << "Processing point interpolation" << endl;
106  }
107 
108 
109  // Check for command line point specification if no .pts file specified
110  if(m_f->m_fieldPts == LibUtilities::NullPtsField)
111  {
112  if(m_config["line"].as<string>().compare("NotSet") != 0)
113  {
114  string help = m_config["line"].as<string>();
115  vector<NekDouble> values;
117  m_config["line"].as<string>().c_str(),values),
118  "Failed to interpret line string");
119 
120  ASSERTL0(values.size() > 3,
121  "line string should contain 2Dim+1 values "
122  "N,x0,y0,z0,x1,y1,z1");
123 
124  int dim = (values.size()-1)/2;
125  int npts = values[0];
127 
128  for(int i = 0; i < dim; ++i)
129  {
130  pts[i] = Array<OneD,NekDouble>(npts);
131  }
132 
133  for(int i = 0; i < npts; ++i)
134  {
135  pts[0][i] = values[1]
136  + i/((NekDouble)(npts-1))*(values[dim+1] - values[1]);
137  if(dim > 1)
138  {
139  pts[1][i] = values[2]
140  + i/((NekDouble)(npts-1))*(values[dim+2] - values[2]);
141 
142  if(dim > 2)
143  {
144  pts[2][i] = values[3]
145  + i/((NekDouble)(npts-1))*(values[dim+3] - values[3]);
146  }
147  }
148  }
149 
150  vector<int> ppe;
151  ppe.push_back(npts);
153  m_f->m_fieldPts->SetPointsPerEdge(ppe);
154  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsLine);
155 
156  }
157  else if(m_config["plane"].as<string>().compare("NotSet") != 0)
158  {
159  string help = m_config["plane"].as<string>();
160  vector<NekDouble> values;
162  m_config["plane"].as<string>().c_str(),values),
163  "Failed to interpret plane string");
164 
165  ASSERTL0(values.size() > 3,
166  "line string should contain 2Dim+1 values "
167  "N,x0,y0,z0,x1,y1,z1");
168 
169 
170  int dim = (values.size()-1)/4;
171 
172  int npts1 = values[0];
173  int npts2 = values[1];
174 
176 
177  for(int i = 0; i < dim; ++i)
178  {
179  pts[i] = Array<OneD,NekDouble>(npts1*npts2);
180  }
181 
182  for(int j = 0; j < npts2; ++j)
183  {
184  for(int i = 0; i < npts1; ++i)
185  {
186  pts[0][i+j*npts1] =
187  (values[2] + i/((NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((NekDouble)(npts2-1))) +
188  (values[3*dim+2] + i/((NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((NekDouble)(npts2-1)));
189  if(dim > 1)
190  {
191  pts[1][i+j*npts1] =
192  (values[3] + i/((NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((NekDouble)(npts2-1))) +
193  (values[3*dim+3] + i/((NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((NekDouble)(npts2-1)));
194 
195  if(dim > 2)
196  {
197  pts[2][i+j*npts1] =
198  (values[4] + i/((NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((NekDouble)(npts2-1))) +
199  (values[3*dim+4] + i/((NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((NekDouble)(npts2-1)));
200  }
201  }
202  }
203  }
204 
205  vector<int> ppe;
206  ppe.push_back(npts1);
207  ppe.push_back(npts2);
209  m_f->m_fieldPts->SetPointsPerEdge(ppe);
210  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsPlane);
211 
212  }
213  }
214 
215 
216  FieldSharedPtr fromField = boost::shared_ptr<Field>(new Field());
217 
218  std::vector<std::string> files;
219  // set up session file for from field
220  files.push_back(m_config["fromxml"].as<string>());
221  fromField->m_session = LibUtilities::SessionReader::
222  CreateInstance(0, 0, files);
223 
224  // Set up range based on min and max of local parallel partition
226 
227  int coordim = m_f->m_fieldPts->GetDim();
228  int npts = m_f->m_fieldPts->GetNpoints();
230  m_f->m_fieldPts->GetPts(pts);
231 
232  rng->m_checkShape = false;
233  switch(coordim)
234  {
235  case 3:
236  rng->m_doZrange = true;
237  rng->m_zmin = Vmath::Vmin(npts, pts[2],1);
238  rng->m_zmax = Vmath::Vmax(npts, pts[2],1);
239  case 2:
240  rng->m_doYrange = true;
241  rng->m_ymin = Vmath::Vmin(npts, pts[1],1);
242  rng->m_ymax = Vmath::Vmax(npts, pts[1],1);
243  case 1:
244  rng->m_doXrange = true;
245  rng->m_xmin = Vmath::Vmin(npts, pts[0],1);
246  rng->m_xmax = Vmath::Vmax(npts, pts[0],1);
247  break;
248  default:
249  ASSERTL0(false,"too many values specfied in range");
250  }
251 
252  // setup rng parameters.
253  fromField->m_graph = SpatialDomains::MeshGraph::Read(fromField->m_session,rng);
254 
255  // Read in local from field partitions
256  const SpatialDomains::ExpansionMap &expansions = fromField->m_graph->GetExpansions();
257 
258  fromField->m_fld = MemoryManager<LibUtilities::FieldIO>
259  ::AllocateSharedPtr(fromField->m_session->GetComm());
260 
261  Array<OneD,int> ElementGIDs(expansions.size());
262  SpatialDomains::ExpansionMap::const_iterator expIt;
263 
264  int i = 0;
265  for (expIt = expansions.begin(); expIt != expansions.end();
266  ++expIt)
267  {
268  ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
269  }
270 
271  string fromfld = m_config["fromfld"].as<string>();
272  fromField->m_fld->Import(fromfld,fromField->m_fielddef,
273  fromField->m_data,
275  ElementGIDs);
276 
277  int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
278 
279  //----------------------------------------------
280  // Set up Expansion information to use mode order from field
281  fromField->m_graph->SetExpansions(fromField->m_fielddef);
282 
283  int nfields = fromField->m_fielddef[0]->m_fields.size();
284 
285  fromField->m_exp.resize(nfields);
286  fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,true);
287 
288  m_f->m_exp.resize(nfields);
289 
290  // declare auxiliary fields.
291  for(i = 1; i < nfields; ++i)
292  {
293  fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
294  }
295 
296  // load field into expansion in fromfield.
297  for(int j = 0; j < nfields; ++j)
298  {
299  for (i = 0; i < fromField->m_fielddef.size(); i++)
300  {
301  fromField->m_exp[j]->ExtractDataToCoeffs(
302  fromField->m_fielddef[i],
303  fromField->m_data[i],
304  fromField->m_fielddef[0]->m_fields[j],
305  fromField->m_exp[j]->UpdateCoeffs());
306  }
307  fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
308  fromField->m_exp[j]->UpdatePhys());
309 
310  Array< OneD, NekDouble > newPts(m_f->m_fieldPts->GetNpoints());
311  m_f->m_fieldPts->AddField(newPts, fromField->m_fielddef[0]->m_fields[j]);
312  }
313 
314  if(fromField->m_session->GetComm()->GetRank() == 0)
315  {
316  cout << "Interpolating [" << flush;
317  }
318 
319  NekDouble clamp_low = m_config["clamptolowervalue"].as<NekDouble>();
320  NekDouble clamp_up = m_config["clamptouppervalue"].as<NekDouble>();
321  NekDouble def_value = m_config["defaultvalue"].as<NekDouble>();
322 
323  InterpolateFieldToPts(fromField->m_exp, pts,
324  clamp_low, clamp_up, def_value);
325 
326  if(fromField->m_session->GetComm()->GetRank() == 0)
327  {
328  cout << "]" << endl;
329  }
330 
331 }
332 
334  vector<MultiRegions::ExpListSharedPtr> &field0,
336  NekDouble clamp_low,
337  NekDouble clamp_up,
338  NekDouble def_value)
339 {
340  int expdim = field0[0]->GetCoordim(0);
341 
342  Array<OneD, NekDouble> coords(expdim), Lcoords(expdim);
343  int nq1 = pts[0].num_elements();
344  int elmtid, offset;
345  int r, f;
346  int intpts = 0;
347 
348  // resize data field
349  m_f->m_data.resize(field0.size());
350 
351  for (f = 0; f < field0.size(); ++f)
352  {
353  m_f->m_data[f].resize(nq1);
354  }
355 
356  for (r = 0; r < nq1; r++)
357  {
358  coords[0] = pts[0][r];
359  coords[1] = pts[1][r];
360  if (expdim == 3)
361  {
362  coords[2] = pts[2][r];
363  }
364 
365  // Obtain Element and LocalCoordinate to interpolate
366  elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3);
367 
368  if(elmtid >= 0)
369  {
370  offset = field0[0]->GetPhys_Offset(field0[0]->
371  GetOffset_Elmt_Id(elmtid));
372 
373  for (f = 0; f < field0.size(); ++f)
374  {
375  NekDouble value;
376  value = field0[f]->GetExp(elmtid)->
377  StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
378 
379  if ((boost::math::isnan)(value))
380  {
381  ASSERTL0(false, "new value is not a number");
382  }
383  else
384  {
385  value = (value > clamp_up)? clamp_up :
386  ((value < clamp_low)? clamp_low :
387  value);
388 
389  m_f->m_data[f][r] = value;
390  }
391  }
392  }
393  else
394  {
395  for (f = 0; f < field0.size(); ++f)
396  {
397  m_f->m_data[f][r] = def_value;
398  }
399  }
400 
401  if (intpts%1000 == 0)
402  {
403  cout <<"." << flush;
404  }
405  intpts ++;
406  }
407 }
408 
409 }
410 }
411 
412 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:119
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:756
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:848
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
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:154
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:677
Represents a command-line configuration option.
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:249
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
Definition: ParseUtils.hpp:127
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:66
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:171
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
void InterpolateFieldToPts(vector< MultiRegions::ExpListSharedPtr > &field0, Array< OneD, Array< OneD, NekDouble > > &pts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215