Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessMapping.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessMapping.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: Add mapping coordinates to field
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessMapping.h"
41 
44 
45 namespace Nektar
46 {
47 namespace Utilities
48 {
49 ModuleKey ProcessMapping::className =
51  ModuleKey(eProcessModule, "mapping"),
52  ProcessMapping::create, "Add mapping coordinates to output file.");
53 
54 ProcessMapping::ProcessMapping(FieldSharedPtr f) : ProcessModule(f)
55 {
56 }
57 
59 {
60 }
61 
62 void ProcessMapping::Process(po::variables_map &vm)
63 {
64  if(m_f->m_verbose)
65  {
66  if(m_f->m_comm->GetRank() == 0)
67  {
68  cout << "ProcessMapping: Applying mapping to field..." << endl;
69  }
70  }
71 
72  // Determine dimensions of mesh, solution, etc...
73  int npoints = m_f->m_exp[0]->GetNpoints();
74  int expdim = m_f->m_graph->GetMeshDimension();
75  int spacedim = expdim;
76  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
77  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
78  {
79  spacedim = 3;
80  }
81  int nfields = m_f->m_fielddef[0]->m_fields.size();
82  int addfields = expdim;
83  m_f->m_exp.resize(nfields+addfields);
84 
85  // Load mapping
87 
88  // Convert velocity to Cartesian system
89  if (m_f->m_fieldMetaDataMap.count("MappingCartesianVel"))
90  {
91  if (m_f->m_fieldMetaDataMap["MappingCartesianVel"] == "False")
92  {
93  m_f->m_fieldMetaDataMap["MappingCartesianVel"] = "True";
94 
95  Array<OneD, Array<OneD, NekDouble> > vel (spacedim);
96  // Initialize arrays and copy velocity
97  for ( int i =0; i<spacedim; ++i )
98  {
99  vel[i] = Array<OneD, NekDouble> (npoints);
100  if (m_f->m_exp[0]->GetWaveSpace())
101  {
102  m_f->m_exp[0]->HomogeneousBwdTrans(m_f->m_exp[i]->GetPhys(),
103  vel[i]);
104  }
105  else
106  {
107  Vmath::Vcopy(npoints, m_f->m_exp[i]->GetPhys(), 1,
108  vel[i], 1);
109  }
110 
111  }
112  // Convert velocity to cartesian system
113  mapping->ContravarToCartesian(vel, vel);
114  // Copy result back
115  for ( int i =0; i<spacedim; ++i )
116  {
117  if (m_f->m_exp[0]->GetWaveSpace())
118  {
119  m_f->m_exp[0]->HomogeneousFwdTrans(vel[i],
120  m_f->m_exp[i]->UpdatePhys());
121  }
122  else
123  {
124  Vmath::Vcopy(npoints, vel[i], 1,
125  m_f->m_exp[i]->UpdatePhys(), 1);
126  }
127  m_f->m_exp[i]->FwdTrans_IterPerExp(m_f->m_exp[i]->GetPhys(),
128  m_f->m_exp[i]->UpdateCoeffs());
129  }
130  }
131  }
132 
133  // Get coordinates from mapping
135  mapping->GetCartesianCoordinates(coords[0], coords[1], coords[2]);
136 
137  // Add new information to m_f
138  string fieldNames[3] = {"xCoord", "yCoord", "zCoord"};
139  vector<string > outname;
140  for (int i = 0; i < addfields; ++i)
141  {
142  m_f->m_exp[nfields + i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
143  m_f->m_exp[nfields + i]->UpdatePhys() = coords[i];
144  m_f->m_exp[nfields + i]->FwdTrans_IterPerExp(coords[i],
145  m_f->m_exp[nfields + i]->UpdateCoeffs());
146  outname.push_back(fieldNames[i]);
147  }
148 
149  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
150  = m_f->m_exp[0]->GetFieldDefinitions();
151  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
152 
153  for (int j = 0; j < nfields + addfields; ++j)
154  {
155  for (int i = 0; i < FieldDef.size(); ++i)
156  {
157  if (j >= nfields)
158  {
159  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
160  }
161  else
162  {
163  FieldDef[i]->m_fields.push_back(m_f->m_fielddef[0]->m_fields[j]);
164  }
165  m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
166  }
167  }
168 
169  m_f->m_fielddef = FieldDef;
170  m_f->m_data = FieldData;
171 }
172 
174 {
175  // Create mapping object
177  field[0] = f->m_exp[0];
179  GlobalMapping::Mapping::Load(f->m_session,
180  field);
181 
182  // Get time from metadata
183  NekDouble time;
184  if (f->m_fieldMetaDataMap.count("Time"))
185  {
186  string s_time = f->m_fieldMetaDataMap["Time"];
187  time = atof(s_time.c_str());
188  }
189  else
190  {
191  time = 0.0;
192  }
193 
194  // Get field information
195  int npoints = f->m_exp[0]->GetNpoints();
196  int expdim = f->m_graph->GetMeshDimension();
197  int spacedim = expdim;
198  if ((f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
199  (f->m_fielddef[0]->m_numHomogeneousDir) == 2)
200  {
201  spacedim = 3;
202  }
203 
204  // Declare coordinates storage
205  Array<OneD, Array<OneD, NekDouble> > coords_new(3);
206  Array<OneD, Array<OneD, NekDouble> > coords_vel(3);
207  for (int i = 0; i < 3; i++)
208  {
209  coords_new[i] = Array<OneD, NekDouble> (npoints);
210  coords_vel[i] = Array<OneD, NekDouble> (npoints,0.0);
211  }
212 
213  string fieldNames[3] = {"x", "y", "z"};
214  string velFieldNames[3] = {"vx", "vy", "vz"};
215 
216  // Evaluate coordinates and coordinates velocity
217  if (f->m_fieldMetaDataMap.count("MappingType"))
218  {
219  if (f->m_fieldMetaDataMap["MappingType"] == "Expression")
220  {
221  // Get name of the functions
222  string funcName;
223  string velFuncName;
224  if (f->m_fieldMetaDataMap.count("MappingExpression"))
225  {
226  funcName = f->m_fieldMetaDataMap["MappingExpression"];
227  }
228  else
229  {
230  funcName = "";
231  }
232  if (f->m_fieldMetaDataMap.count("MappingVelExpression"))
233  {
234  velFuncName = f->m_fieldMetaDataMap["MappingVelExpression"];
235  }
236  else
237  {
238  velFuncName = "";
239  }
240 
241  // Get original coordinates (in case some of them are not changed)
243  for (int i = 0; i < 3; i++)
244  {
245  coords[i] = Array<OneD, NekDouble> (npoints);
246  }
247  f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
248 
249  // Load coordinates
250  std::string s_FieldStr;
251  for(int i = 0; i < 3; i++)
252  {
253  s_FieldStr = fieldNames[i];
254  if ( f->m_session->DefinesFunction(funcName, s_FieldStr))
255  {
257  f->m_session->GetFunction(funcName, s_FieldStr);
258  ffunc->Evaluate(coords[0], coords[1], coords[2],
259  time, coords_new[i]);
260  }
261  else
262  {
263  // This coordinate is not defined, so use (x^i)' = x^i
264  Vmath::Vcopy(npoints, coords[i], 1, coords_new[i], 1);
265  }
266  }
267  // Load velocities
268  if (f->m_session->DefinesFunction(velFuncName))
269  {
270  for (int i = 0; i < 3; i++)
271  {
272  s_FieldStr = velFieldNames[i];
273  if ( f->m_session->DefinesFunction(velFuncName, s_FieldStr))
274  {
276  f->m_session->GetFunction(velFuncName, s_FieldStr);
277  ffunc->Evaluate(coords[0], coords[1], coords[2],
278  time, coords_vel[i]);
279  }
280  }
281  }
282 
283  // Update mapping with coordinates
284  mapping->SetFromFunction(false);
285  mapping->UpdateMapping(time, coords_new,coords_vel);
286  }
287  else if(f->m_fieldMetaDataMap["MappingType"] == "File")
288  {
289  ASSERTL0(f->m_fieldMetaDataMap.count("FileName"),
290  "FileName parameter for Mapping missing in field file.");
291  string fileName = f->m_fieldMetaDataMap["FileName"];
292  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
293  std::vector<std::vector<NekDouble> > FieldData;
294 
295  f->m_fld->Import(fileName,
296  FieldDef,
297  FieldData);
298 
299  for (int j = 0; j < spacedim; ++j)
300  {
301  int ncoeffs = f->m_exp[0]->GetNcoeffs();
302  Array<OneD, NekDouble> fieldcoeffs(ncoeffs,0.0);
303  for (int i = 0; i < FieldData.size(); ++i)
304  {
305  f->m_exp[j]->ExtractDataToCoeffs(FieldDef[i],
306  FieldData[i],
307  fieldNames[j],
308  fieldcoeffs);
309  }
310  bool wavespace = f->m_exp[0]->GetWaveSpace();
311  f->m_exp[0]->SetWaveSpace(false);
312 
313  f->m_exp[0]->BwdTrans(fieldcoeffs,
314  coords_new[j]);
315 
316 
317 
318  // Load coordinate velocity
319  if ( std::find(FieldDef[0]->m_fields.begin(), FieldDef[0]->m_fields.end(),
320  velFieldNames[j])!=FieldDef[0]->m_fields.end())
321  {
322  for (int i = 0; i < FieldData.size(); ++i)
323  {
324  f->m_exp[j]->ExtractDataToCoeffs(FieldDef[i],
325  FieldData[i],
326  velFieldNames[j],
327  fieldcoeffs);
328  }
329  f->m_exp[0]->BwdTrans(fieldcoeffs,
330  coords_vel[j]);
331  }
332  f->m_exp[0]->SetWaveSpace(wavespace);
333  }
334  // Update mapping with coordinates
335  mapping->SetFromFunction(false);
336  mapping->UpdateMapping(time, coords_new,coords_vel);
337  }
338  }
339  else
340  {
341  // Use trivial mapping
343  Array<OneD, Array<OneD, NekDouble> > coords_vel(3);
344  for (int i = 0; i < 3; i++)
345  {
346  coords[i] = Array<OneD, NekDouble> (npoints);
347  coords_vel[i] = Array<OneD, NekDouble> (npoints, 0.0);
348  }
349  f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
350  mapping->SetFromFunction(false);
351  mapping->UpdateMapping(time, coords,coords_vel);
352  }
353 
354  return mapping;
355 
356 }
357 
358 }
359 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
pair< ModuleType, string > ModuleKey
virtual void Process()=0
STL namespace.
FieldSharedPtr m_f
Field object.
static GlobalMapping::MappingSharedPtr GetMapping(FieldSharedPtr f)
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
GLOBAL_MAPPING_EXPORT typedef boost::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:51
boost::shared_ptr< Equation > EquationSharedPtr
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:266
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