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  if(m_f->m_verbose)
80  {
81  if(m_f->m_comm->GetRank() == 0)
82  {
83  cout << "ProcessInterpField: Interpolating field..." << endl;
84  }
85  }
86 
87  m_fromField = boost::shared_ptr<Field>(new Field());
88 
89  std::vector<std::string> files;
90 
91  // set up session file for from field
92  files.push_back(m_config["fromxml"].as<string>());
94  CreateInstance(0, 0, files);
95 
96  // Set up range based on min and max of local parallel partition
99 
100  int coordim = m_f->m_exp[0]->GetCoordim(0);
101  int npts = m_f->m_exp[0]->GetTotPoints();
103 
104 
105  for(int i = 0; i < coordim; ++i)
106  {
107  coords[i] = Array<OneD, NekDouble>(npts);
108  }
109 
110  for(int i = coordim; i < 3; ++i)
111  {
112  coords[i] = NullNekDouble1DArray;
113  }
114 
115  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
116 
117  rng->m_checkShape = false;
118  switch(coordim)
119  {
120  case 3:
121  rng->m_doZrange = true;
122  rng->m_zmin = Vmath::Vmin(npts,coords[2],1);
123  rng->m_zmax = Vmath::Vmax(npts,coords[2],1);
124  case 2:
125  rng->m_doYrange = true;
126  rng->m_ymin = Vmath::Vmin(npts,coords[1],1);
127  rng->m_ymax = Vmath::Vmax(npts,coords[1],1);
128  case 1:
129  rng->m_doXrange = true;
130  rng->m_xmin = Vmath::Vmin(npts,coords[0],1);
131  rng->m_xmax = Vmath::Vmax(npts,coords[0],1);
132  break;
133  default:
134  ASSERTL0(false,"too many values specfied in range");
135  }
136 
137  // setup rng parameters.
138  m_fromField->m_graph =
140 
141  // Read in local from field partitions
142  const SpatialDomains::ExpansionMap &expansions =
143  m_fromField->m_graph->GetExpansions();
144 
145  // check for case where no elements are specified on this
146  // parallel partition
147  if(!expansions.size())
148  {
149  return;
150  }
151 
152  Array<OneD,int> ElementGIDs(expansions.size());
153  SpatialDomains::ExpansionMap::const_iterator expIt;
154 
155  int i = 0;
156  for (expIt = expansions.begin(); expIt != expansions.end();
157  ++expIt)
158  {
159  ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
160  }
161 
162  string fromfld = m_config["fromfld"].as<string>();
163  m_f->m_fld->Import(fromfld,m_fromField->m_fielddef,
164  m_fromField->m_data,
166  ElementGIDs);
167 
168  int NumHomogeneousDir = m_fromField->m_fielddef[0]->m_numHomogeneousDir;
169 
170  //----------------------------------------------
171  // Set up Expansion information to use mode order from field
172  m_fromField->m_graph->SetExpansions(m_fromField->m_fielddef);
173 
174  int nfields = m_fromField->m_fielddef[0]->m_fields.size();
175 
176  m_fromField->m_exp.resize(nfields);
177  m_fromField->m_exp[0] = m_fromField->SetUpFirstExpList(NumHomogeneousDir,true);
178 
179  m_f->m_exp.resize(nfields);
180 
181  // declare auxiliary fields.
182  for(i = 1; i < nfields; ++i)
183  {
184  m_f->m_exp[i] = m_f->AppendExpList(NumHomogeneousDir);
185  m_fromField->m_exp[i] = m_fromField->AppendExpList(NumHomogeneousDir);
186  }
187 
188 
189  // load field into expansion in fromfield.
190  for(int j = 0; j < nfields; ++j)
191  {
192  for (i = 0; i < m_fromField->m_fielddef.size(); i++)
193  {
194  m_fromField->m_exp[j]->ExtractDataToCoeffs(
195  m_fromField->m_fielddef[i],
196  m_fromField->m_data[i],
197  m_fromField->m_fielddef[0]->m_fields[j],
198  m_fromField->m_exp[j]->UpdateCoeffs());
199  }
200  m_fromField->m_exp[j]->BwdTrans(m_fromField->m_exp[j]->GetCoeffs(),
201  m_fromField->m_exp[j]->UpdatePhys());
202 
203  }
204 
205  int nq1 = m_f->m_exp[0]->GetTotPoints();
206 
207  Array<OneD, NekDouble> x1(nq1);
208  Array<OneD, NekDouble> y1(nq1);
209  Array<OneD, NekDouble> z1(nq1);
210 
211  if (coordim == 2)
212  {
213  m_f->m_exp[0]->GetCoords(x1, y1);
214  }
215  else if (coordim == 3)
216  {
217  m_f->m_exp[0]->GetCoords(x1, y1, z1);
218  }
219 
220  if(m_f->m_session->GetComm()->TreatAsRankZero())
221  {
222  cout << "Interpolating [" << flush;
223  }
224 
225  NekDouble clamp_low = m_config["clamptolowervalue"].as<NekDouble>();
226  NekDouble clamp_up = m_config["clamptouppervalue"].as<NekDouble>();
227  NekDouble def_value = m_config["defaultvalue"].as<NekDouble>();
228 
229  InterpolateField(m_fromField->m_exp, m_f->m_exp,
230  x1, y1, z1, clamp_low, clamp_up,def_value);
231 
232  if(m_f->m_session->GetComm()->TreatAsRankZero())
233  {
234  cout << "]" << endl;
235  }
236 
237 
238  // put field into field data for output
239  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
240  = m_f->m_exp[0]->GetFieldDefinitions();
241  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
242 
243  for (int j = 0; j < nfields; ++j)
244  {
245  m_f->m_exp[j]->FwdTrans(m_f->m_exp[j]->GetPhys(),
246  m_f->m_exp[j]->UpdateCoeffs());
247  for (i = 0; i < FieldDef.size(); ++i)
248  {
249  FieldDef[i]->m_fields.push_back(m_fromField->m_fielddef[0]->m_fields[j]);
250  m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
251  }
252  }
253 
254  m_f->m_fielddef = FieldDef;
255  m_f->m_data = FieldData;
256 }
257 
259  vector<MultiRegions::ExpListSharedPtr> &field0,
260  vector<MultiRegions::ExpListSharedPtr> &field1,
264  NekDouble clamp_low,
265  NekDouble clamp_up,
266  NekDouble def_value)
267 {
268  int expdim = field0[0]->GetCoordim(0);
269 
270  Array<OneD, NekDouble> coords(expdim), Lcoords(expdim);
271  int nq1 = field1[0]->GetTotPoints();
272  int elmtid, offset;
273  int r, f;
274  static int intpts = 0;
275 
276  ASSERTL0(field0.size() == field1.size(),
277  "Input field dimension must be same as output dimension");
278 
279  for (r = 0; r < nq1; r++)
280  {
281  coords[0] = x[r];
282  coords[1] = y[r];
283  if (expdim == 3)
284  {
285  coords[2] = z[r];
286  }
287 
288  // Obtain nearest Element and LocalCoordinate to interpolate
289  elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3,true);
290 
291  if(elmtid >= 0)
292  {
293  offset = field0[0]->GetPhys_Offset(field0[0]->
294  GetOffset_Elmt_Id(elmtid));
295 
296  for (f = 0; f < field1.size(); ++f)
297  {
298  NekDouble value;
299  value = field0[f]->GetExp(elmtid)->
300  StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
301 
302  if ((boost::math::isnan)(value))
303  {
304  ASSERTL0(false, "new value is not a number");
305  }
306  else
307  {
308  if(value > clamp_up)
309  {
310  value = clamp_up;
311  }
312  else if( value < clamp_low)
313  {
314  value = clamp_low;
315  }
316 
317  field1[f]->UpdatePhys()[r] = value;
318  }
319  }
320  }
321  else
322  {
323  for (f = 0; f < field1.size(); ++f)
324  {
325  field1[f]->UpdatePhys()[r] = def_value;
326  }
327  }
328 
329  if (intpts%1000 == 0)
330  {
331  cout <<"." << flush;
332  }
333  intpts ++;
334  }
335 }
336 
337 }
338 }
339 
340 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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