Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessInterpField.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessInterpField.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 one field to another.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 #include <string>
36 #include <iostream>
37 using namespace std;
38 
39 #include "ProcessInterpField.h"
40 
43 #include <boost/math/special_functions/fpclassify.hpp>
44 namespace Nektar
45 {
46 namespace Utilities
47 {
48 
49 ModuleKey ProcessInterpField::className =
51  ModuleKey(eProcessModule, "interpfield"),
52  ProcessInterpField::create,
53  "Interpolates one field to another, requires fromxml, "
54  "fromfld to be defined");
55 
56 
57 ProcessInterpField::ProcessInterpField(FieldSharedPtr f) : ProcessModule(f)
58 {
59 
60  m_config["fromxml"] = ConfigOption(false, "NotSet",
61  "Xml file form which to interpolate field");
62  m_config["fromfld"] = ConfigOption(false, "NotSet",
63  "Fld file form which to interpolate field");
64 
65  m_config["clamptolowervalue"] = ConfigOption(false,"-10000000",
66  "Lower bound for interpolation value");
67  m_config["clamptouppervalue"] = ConfigOption(false,"10000000",
68  "Upper bound for interpolation value");
69  m_config["defaultvalue"] = ConfigOption(false,"0",
70  "Default value if point is outside domain");
71 }
72 
74 {
75 }
76 
77 void ProcessInterpField::Process(po::variables_map &vm)
78 {
79 
80  if(m_f->m_verbose)
81  {
82  cout << "Processing interpolation" << endl;
83  }
84 
85  m_fromField = boost::shared_ptr<Field>(new Field());
86 
87  std::vector<std::string> files;
88 
89  // set up session file for from field
90  files.push_back(m_config["fromxml"].as<string>());
92  CreateInstance(0, 0, files);
93 
94  // Set up range based on min and max of local parallel partition
97 
98  int coordim = m_f->m_exp[0]->GetCoordim(0);
99  int npts = m_f->m_exp[0]->GetTotPoints();
101 
102 
103  for(int i = 0; i < coordim; ++i)
104  {
105  coords[i] = Array<OneD, NekDouble>(npts);
106  }
107 
108  for(int i = coordim; i < 3; ++i)
109  {
110  coords[i] = NullNekDouble1DArray;
111  }
112 
113  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
114 
115  rng->m_checkShape = false;
116  switch(coordim)
117  {
118  case 3:
119  rng->m_doZrange = true;
120  rng->m_zmin = Vmath::Vmin(npts,coords[2],1);
121  rng->m_zmax = Vmath::Vmax(npts,coords[2],1);
122  case 2:
123  rng->m_doYrange = true;
124  rng->m_ymin = Vmath::Vmin(npts,coords[1],1);
125  rng->m_ymax = Vmath::Vmax(npts,coords[1],1);
126  case 1:
127  rng->m_doXrange = true;
128  rng->m_xmin = Vmath::Vmin(npts,coords[0],1);
129  rng->m_xmax = Vmath::Vmax(npts,coords[0],1);
130  break;
131  default:
132  ASSERTL0(false,"too many values specfied in range");
133  }
134 
135  // setup rng parameters.
136  m_fromField->m_graph =
138 
139  // Read in local from field partitions
140  const SpatialDomains::ExpansionMap &expansions =
141  m_fromField->m_graph->GetExpansions();
142 
143  // check for case where no elements are specified on this
144  // parallel partition
145  if(!expansions.size())
146  {
147  return;
148  }
149 
150  Array<OneD,int> ElementGIDs(expansions.size());
151  SpatialDomains::ExpansionMap::const_iterator expIt;
152 
153  int i = 0;
154  for (expIt = expansions.begin(); expIt != expansions.end();
155  ++expIt)
156  {
157  ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
158  }
159 
160  string fromfld = m_config["fromfld"].as<string>();
161  m_f->m_fld->Import(fromfld,m_fromField->m_fielddef,
162  m_fromField->m_data,
164  ElementGIDs);
165 
166  int NumHomogeneousDir = m_fromField->m_fielddef[0]->m_numHomogeneousDir;
167 
168  //----------------------------------------------
169  // Set up Expansion information to use mode order from field
170  m_fromField->m_graph->SetExpansions(m_fromField->m_fielddef);
171 
172  int nfields = m_fromField->m_fielddef[0]->m_fields.size();
173 
174  m_fromField->m_exp.resize(nfields);
175  m_fromField->m_exp[0] = m_fromField->SetUpFirstExpList(NumHomogeneousDir,true);
176 
177  m_f->m_exp.resize(nfields);
178 
179  // declare auxiliary fields.
180  for(i = 1; i < nfields; ++i)
181  {
182  m_f->m_exp[i] = m_f->AppendExpList(NumHomogeneousDir);
183  m_fromField->m_exp[i] = m_fromField->AppendExpList(NumHomogeneousDir);
184  }
185 
186 
187  // load field into expansion in fromfield.
188  for(int j = 0; j < nfields; ++j)
189  {
190  for (i = 0; i < m_fromField->m_fielddef.size(); i++)
191  {
192  m_fromField->m_exp[j]->ExtractDataToCoeffs(
193  m_fromField->m_fielddef[i],
194  m_fromField->m_data[i],
195  m_fromField->m_fielddef[0]->m_fields[j],
196  m_fromField->m_exp[j]->UpdateCoeffs());
197  }
198  m_fromField->m_exp[j]->BwdTrans(m_fromField->m_exp[j]->GetCoeffs(),
199  m_fromField->m_exp[j]->UpdatePhys());
200 
201  }
202 
203  int nq1 = m_f->m_exp[0]->GetTotPoints();
204 
205  Array<OneD, NekDouble> x1(nq1);
206  Array<OneD, NekDouble> y1(nq1);
207  Array<OneD, NekDouble> z1(nq1);
208 
209  if (coordim == 2)
210  {
211  m_f->m_exp[0]->GetCoords(x1, y1);
212  }
213  else if (coordim == 3)
214  {
215  m_f->m_exp[0]->GetCoords(x1, y1, z1);
216  }
217 
218  if(m_f->m_session->GetComm()->TreatAsRankZero())
219  {
220  cout << "Interpolating [" << flush;
221  }
222 
223  NekDouble clamp_low = m_config["clamptolowervalue"].as<NekDouble>();
224  NekDouble clamp_up = m_config["clamptouppervalue"].as<NekDouble>();
225  NekDouble def_value = m_config["defaultvalue"].as<NekDouble>();
226 
227  InterpolateField(m_fromField->m_exp, m_f->m_exp,
228  x1, y1, z1, clamp_low, clamp_up,def_value);
229 
230  if(m_f->m_session->GetComm()->TreatAsRankZero())
231  {
232  cout << "]" << endl;
233  }
234 
235 
236  // put field into field data for output
237  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
238  = m_f->m_exp[0]->GetFieldDefinitions();
239  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
240 
241  for (int j = 0; j < nfields; ++j)
242  {
243  m_f->m_exp[j]->FwdTrans(m_f->m_exp[j]->GetPhys(),
244  m_f->m_exp[j]->UpdateCoeffs());
245  for (i = 0; i < FieldDef.size(); ++i)
246  {
247  FieldDef[i]->m_fields.push_back(m_fromField->m_fielddef[0]->m_fields[j]);
248  m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
249  }
250  }
251 
252  m_f->m_fielddef = FieldDef;
253  m_f->m_data = FieldData;
254 }
255 
257  vector<MultiRegions::ExpListSharedPtr> &field0,
258  vector<MultiRegions::ExpListSharedPtr> &field1,
262  NekDouble clamp_low,
263  NekDouble clamp_up,
264  NekDouble def_value)
265 {
266  int expdim = field0[0]->GetCoordim(0);
267 
268  Array<OneD, NekDouble> coords(expdim), Lcoords(expdim);
269  int nq1 = field1[0]->GetTotPoints();
270  int elmtid, offset;
271  int r, f;
272  static int intpts = 0;
273 
274  ASSERTL0(field0.size() == field1.size(),
275  "Input field dimension must be same as output dimension");
276 
277  for (r = 0; r < nq1; r++)
278  {
279  coords[0] = x[r];
280  coords[1] = y[r];
281  if (expdim == 3)
282  {
283  coords[2] = z[r];
284  }
285 
286  // Obtain nearest Element and LocalCoordinate to interpolate
287  elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3,true);
288 
289  if(elmtid >= 0)
290  {
291  offset = field0[0]->GetPhys_Offset(field0[0]->
292  GetOffset_Elmt_Id(elmtid));
293 
294  for (f = 0; f < field1.size(); ++f)
295  {
296  NekDouble value;
297  value = field0[f]->GetExp(elmtid)->
298  StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
299 
300  if ((boost::math::isnan)(value))
301  {
302  ASSERTL0(false, "new value is not a number");
303  }
304  else
305  {
306  if(value > clamp_up)
307  {
308  value = clamp_up;
309  }
310  else if( value < clamp_low)
311  {
312  value = clamp_low;
313  }
314 
315  field1[f]->UpdatePhys()[r] = value;
316  }
317  }
318  }
319  else
320  {
321  for (f = 0; f < field1.size(); ++f)
322  {
323  field1[f]->UpdatePhys()[r] = def_value;
324  }
325  }
326 
327  if (intpts%1000 == 0)
328  {
329  cout <<"." << flush;
330  }
331  intpts ++;
332  }
333 }
334 
335 }
336 }
337 
338 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
static Array< OneD, NekDouble > NullNekDouble1DArray
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
void InterpolateField(vector< MultiRegions::ExpListSharedPtr > &field0, vector< MultiRegions::ExpListSharedPtr > &field1, Array< OneD, NekDouble > x, Array< OneD, NekDouble > y, Array< OneD, NekDouble > z, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value)
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
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:157
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
Represents a command-line configuration option.
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:54
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:174
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215