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  // Determine dimensions of mesh, solution, etc...
65  int npoints = m_f->m_exp[0]->GetNpoints();
66  int expdim = m_f->m_graph->GetMeshDimension();
67  int spacedim = expdim;
68  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
69  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
70  {
71  spacedim = 3;
72  }
73  int nfields = m_f->m_fielddef[0]->m_fields.size();
74  int addfields = expdim;
75  m_f->m_exp.resize(nfields+addfields);
76 
77  // Load mapping
79 
80  // Convert velocity to Cartesian system
81  if (m_f->m_fieldMetaDataMap.count("MappingCartesianVel"))
82  {
83  if (m_f->m_fieldMetaDataMap["MappingCartesianVel"] == "False")
84  {
85  m_f->m_fieldMetaDataMap["MappingCartesianVel"] = "True";
86 
87  Array<OneD, Array<OneD, NekDouble> > vel (spacedim);
88  // Initialize arrays and copy velocity
89  for ( int i =0; i<spacedim; ++i )
90  {
91  vel[i] = Array<OneD, NekDouble> (npoints);
92  if (m_f->m_exp[0]->GetWaveSpace())
93  {
94  m_f->m_exp[0]->HomogeneousBwdTrans(m_f->m_exp[i]->GetPhys(),
95  vel[i]);
96  }
97  else
98  {
99  Vmath::Vcopy(npoints, m_f->m_exp[i]->GetPhys(), 1,
100  vel[i], 1);
101  }
102 
103  }
104  // Convert velocity to cartesian system
105  mapping->ContravarToCartesian(vel, vel);
106  // Copy result back
107  for ( int i =0; i<spacedim; ++i )
108  {
109  if (m_f->m_exp[0]->GetWaveSpace())
110  {
111  m_f->m_exp[0]->HomogeneousFwdTrans(vel[i],
112  m_f->m_exp[i]->UpdatePhys());
113  }
114  else
115  {
116  Vmath::Vcopy(npoints, vel[i], 1,
117  m_f->m_exp[i]->UpdatePhys(), 1);
118  }
119  m_f->m_exp[i]->FwdTrans_IterPerExp(m_f->m_exp[i]->GetPhys(),
120  m_f->m_exp[i]->UpdateCoeffs());
121  }
122  }
123  }
124 
125  // Get coordinates from mapping
127  mapping->GetCartesianCoordinates(coords[0], coords[1], coords[2]);
128 
129  // Add new information to m_f
130  string fieldNames[3] = {"xCoord", "yCoord", "zCoord"};
131  vector<string > outname;
132  for (int i = 0; i < addfields; ++i)
133  {
134  m_f->m_exp[nfields + i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
135  m_f->m_exp[nfields + i]->UpdatePhys() = coords[i];
136  m_f->m_exp[nfields + i]->FwdTrans_IterPerExp(coords[i],
137  m_f->m_exp[nfields + i]->UpdateCoeffs());
138  outname.push_back(fieldNames[i]);
139  }
140 
141  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
142  = m_f->m_exp[0]->GetFieldDefinitions();
143  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
144 
145  for (int j = 0; j < nfields + addfields; ++j)
146  {
147  for (int i = 0; i < FieldDef.size(); ++i)
148  {
149  if (j >= nfields)
150  {
151  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
152  }
153  else
154  {
155  FieldDef[i]->m_fields.push_back(m_f->m_fielddef[0]->m_fields[j]);
156  }
157  m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
158  }
159  }
160 
161  m_f->m_fielddef = FieldDef;
162  m_f->m_data = FieldData;
163 
164 }
165 
167 {
168  // Create mapping object
170  field[0] = f->m_exp[0];
172  GlobalMapping::Mapping::Load(f->m_session,
173  field);
174 
175  // Get time from metadata
176  NekDouble time;
177  if (f->m_fieldMetaDataMap.count("Time"))
178  {
179  string s_time = f->m_fieldMetaDataMap["Time"];
180  time = atof(s_time.c_str());
181  }
182  else
183  {
184  time = 0.0;
185  }
186 
187  // Get field information
188  int npoints = f->m_exp[0]->GetNpoints();
189  int expdim = f->m_graph->GetMeshDimension();
190  int spacedim = expdim;
191  if ((f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
192  (f->m_fielddef[0]->m_numHomogeneousDir) == 2)
193  {
194  spacedim = 3;
195  }
196 
197  // Declare coordinates storage
198  Array<OneD, Array<OneD, NekDouble> > coords_new(3);
199  Array<OneD, Array<OneD, NekDouble> > coords_vel(3);
200  for (int i = 0; i < 3; i++)
201  {
202  coords_new[i] = Array<OneD, NekDouble> (npoints);
203  coords_vel[i] = Array<OneD, NekDouble> (npoints,0.0);
204  }
205 
206  string fieldNames[3] = {"x", "y", "z"};
207  string velFieldNames[3] = {"vx", "vy", "vz"};
208 
209  // Evaluate coordinates and coordinates velocity
210  if (f->m_fieldMetaDataMap.count("MappingType"))
211  {
212  if (f->m_fieldMetaDataMap["MappingType"] == "Expression")
213  {
214  // Get name of the functions
215  string funcName;
216  string velFuncName;
217  if (f->m_fieldMetaDataMap.count("MappingExpression"))
218  {
219  funcName = f->m_fieldMetaDataMap["MappingExpression"];
220  }
221  else
222  {
223  funcName = "";
224  }
225  if (f->m_fieldMetaDataMap.count("MappingVelExpression"))
226  {
227  velFuncName = f->m_fieldMetaDataMap["MappingVelExpression"];
228  }
229  else
230  {
231  velFuncName = "";
232  }
233 
234  // Get original coordinates (in case some of them are not changed)
236  for (int i = 0; i < 3; i++)
237  {
238  coords[i] = Array<OneD, NekDouble> (npoints);
239  }
240  f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
241 
242  // Load coordinates
243  std::string s_FieldStr;
244  for(int i = 0; i < 3; i++)
245  {
246  s_FieldStr = fieldNames[i];
247  if ( f->m_session->DefinesFunction(funcName, s_FieldStr))
248  {
250  f->m_session->GetFunction(funcName, s_FieldStr);
251  ffunc->Evaluate(coords[0], coords[1], coords[2],
252  time, coords_new[i]);
253  }
254  else
255  {
256  // This coordinate is not defined, so use (x^i)' = x^i
257  Vmath::Vcopy(npoints, coords[i], 1, coords_new[i], 1);
258  }
259  }
260  // Load velocities
261  if (f->m_session->DefinesFunction(velFuncName))
262  {
263  for (int i = 0; i < 3; i++)
264  {
265  s_FieldStr = velFieldNames[i];
266  if ( f->m_session->DefinesFunction(velFuncName, s_FieldStr))
267  {
269  f->m_session->GetFunction(velFuncName, s_FieldStr);
270  ffunc->Evaluate(coords[0], coords[1], coords[2],
271  time, coords_vel[i]);
272  }
273  }
274  }
275 
276  // Update mapping with coordinates
277  mapping->SetFromFunction(false);
278  mapping->UpdateMapping(time, coords_new,coords_vel);
279  }
280  else if(f->m_fieldMetaDataMap["MappingType"] == "File")
281  {
282  ASSERTL0(f->m_fieldMetaDataMap.count("FileName"),
283  "FileName parameter for Mapping missing in field file.");
284  string fileName = f->m_fieldMetaDataMap["FileName"];
285  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
286  std::vector<std::vector<NekDouble> > FieldData;
287 
288  f->m_fld->Import(fileName,
289  FieldDef,
290  FieldData);
291 
292  for (int j = 0; j < spacedim; ++j)
293  {
294  int ncoeffs = f->m_exp[0]->GetNcoeffs();
295  Array<OneD, NekDouble> fieldcoeffs(ncoeffs,0.0);
296  for (int i = 0; i < FieldData.size(); ++i)
297  {
298  f->m_exp[j]->ExtractDataToCoeffs(FieldDef[i],
299  FieldData[i],
300  fieldNames[j],
301  fieldcoeffs);
302  }
303  bool wavespace = f->m_exp[0]->GetWaveSpace();
304  f->m_exp[0]->SetWaveSpace(false);
305 
306  f->m_exp[0]->BwdTrans(fieldcoeffs,
307  coords_new[j]);
308 
309 
310 
311  // Load coordinate velocity
312  if ( std::find(FieldDef[0]->m_fields.begin(), FieldDef[0]->m_fields.end(),
313  velFieldNames[j])!=FieldDef[0]->m_fields.end())
314  {
315  for (int i = 0; i < FieldData.size(); ++i)
316  {
317  f->m_exp[j]->ExtractDataToCoeffs(FieldDef[i],
318  FieldData[i],
319  velFieldNames[j],
320  fieldcoeffs);
321  }
322  f->m_exp[0]->BwdTrans(fieldcoeffs,
323  coords_vel[j]);
324  }
325  f->m_exp[0]->SetWaveSpace(wavespace);
326  }
327  // Update mapping with coordinates
328  mapping->SetFromFunction(false);
329  mapping->UpdateMapping(time, coords_new,coords_vel);
330  }
331  }
332  else
333  {
334  // Use trivial mapping
336  Array<OneD, Array<OneD, NekDouble> > coords_vel(3);
337  for (int i = 0; i < 3; i++)
338  {
339  coords[i] = Array<OneD, NekDouble> (npoints);
340  coords_vel[i] = Array<OneD, NekDouble> (npoints, 0.0);
341  }
342  f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
343  mapping->SetFromFunction(false);
344  mapping->UpdateMapping(time, coords,coords_vel);
345  }
346 
347  return mapping;
348 
349 }
350 
351 }
352 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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