Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 <iostream>
37 #include <string>
38 using namespace std;
39 
40 #include "ProcessMapping.h"
41 
44 
45 namespace Nektar
46 {
47 namespace FieldUtils
48 {
49 ModuleKey ProcessMapping::className =
51  ModuleKey(eProcessModule, "mapping"),
52  ProcessMapping::create,
53  "Add mapping coordinates to output file.");
54 
55 ProcessMapping::ProcessMapping(FieldSharedPtr f) : ProcessModule(f)
56 {
57 }
58 
60 {
61 }
62 
63 void ProcessMapping::Process(po::variables_map &vm)
64 {
65  if (m_f->m_verbose)
66  {
67  if (m_f->m_comm->TreatAsRankZero())
68  {
69  cout << "ProcessMapping: Applying mapping to field..." << endl;
70  }
71  }
72 
73  // Determine dimensions of mesh, solution, etc...
74  int npoints = m_f->m_exp[0]->GetNpoints();
75  int expdim = m_f->m_graph->GetMeshDimension();
76  int spacedim = expdim;
77  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
78  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
79  {
80  spacedim = 3;
81  }
82  int nfields = m_f->m_fielddef[0]->m_fields.size();
83  int addfields = expdim;
84  m_f->m_exp.resize(nfields + addfields);
85 
86  // Load mapping
88 
89  // Convert velocity to Cartesian system
90  if (m_f->m_fieldMetaDataMap.count("MappingCartesianVel"))
91  {
92  if (m_f->m_fieldMetaDataMap["MappingCartesianVel"] == "False")
93  {
94  m_f->m_fieldMetaDataMap["MappingCartesianVel"] = "True";
95 
96  Array<OneD, Array<OneD, NekDouble> > vel(spacedim);
97  // Initialize arrays and copy velocity
98  for (int i = 0; i < spacedim; ++i)
99  {
100  vel[i] = Array<OneD, NekDouble>(npoints);
101  if (m_f->m_exp[0]->GetWaveSpace())
102  {
103  m_f->m_exp[0]->HomogeneousBwdTrans(m_f->m_exp[i]->GetPhys(),
104  vel[i]);
105  }
106  else
107  {
108  Vmath::Vcopy(npoints, m_f->m_exp[i]->GetPhys(), 1, vel[i],
109  1);
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(
120  vel[i], 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(
128  m_f->m_exp[i]->GetPhys(), 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] =
143  m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
144  m_f->m_exp[nfields + i]->UpdatePhys() = coords[i];
145  m_f->m_exp[nfields + i]->FwdTrans_IterPerExp(
146  coords[i], m_f->m_exp[nfields + i]->UpdateCoeffs());
147  outname.push_back(fieldNames[i]);
148  }
149 
150  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
151  m_f->m_exp[0]->GetFieldDefinitions();
152  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
153 
154  for (int j = 0; j < nfields + addfields; ++j)
155  {
156  for (int i = 0; i < FieldDef.size(); ++i)
157  {
158  if (j >= nfields)
159  {
160  FieldDef[i]->m_fields.push_back(outname[j - nfields]);
161  }
162  else
163  {
164  FieldDef[i]->m_fields.push_back(
165  m_f->m_fielddef[0]->m_fields[j]);
166  }
167  m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
168  }
169  }
170 
171  m_f->m_fielddef = FieldDef;
172  m_f->m_data = FieldData;
173 }
174 
176 {
177  // Create mapping object
179  field[0] = f->m_exp[0];
181  GlobalMapping::Mapping::Load(f->m_session, field);
182 
183  // Get time from metadata
184  NekDouble time;
185  if (f->m_fieldMetaDataMap.count("Time"))
186  {
187  string s_time = f->m_fieldMetaDataMap["Time"];
188  time = atof(s_time.c_str());
189  }
190  else
191  {
192  time = 0.0;
193  }
194 
195  // Get field information
196  int npoints = f->m_exp[0]->GetNpoints();
197  int expdim = f->m_graph->GetMeshDimension();
198  int spacedim = expdim;
199  if ((f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
200  (f->m_fielddef[0]->m_numHomogeneousDir) == 2)
201  {
202  spacedim = 3;
203  }
204 
205  // Declare coordinates storage
206  Array<OneD, Array<OneD, NekDouble> > coords_new(3);
207  Array<OneD, Array<OneD, NekDouble> > coords_vel(3);
208  for (int i = 0; i < 3; i++)
209  {
210  coords_new[i] = Array<OneD, NekDouble>(npoints);
211  coords_vel[i] = Array<OneD, NekDouble>(npoints, 0.0);
212  }
213 
214  string fieldNames[3] = {"x", "y", "z"};
215  string velFieldNames[3] = {"vx", "vy", "vz"};
216 
217  // Evaluate coordinates and coordinates velocity
218  if (f->m_fieldMetaDataMap.count("MappingType"))
219  {
220  if (f->m_fieldMetaDataMap["MappingType"] == "Expression")
221  {
222  // Get name of the functions
223  string funcName;
224  string velFuncName;
225  if (f->m_fieldMetaDataMap.count("MappingExpression"))
226  {
227  funcName = f->m_fieldMetaDataMap["MappingExpression"];
228  }
229  else
230  {
231  funcName = "";
232  }
233  if (f->m_fieldMetaDataMap.count("MappingVelExpression"))
234  {
235  velFuncName = f->m_fieldMetaDataMap["MappingVelExpression"];
236  }
237  else
238  {
239  velFuncName = "";
240  }
241 
242  // Get original coordinates (in case some of them are not changed)
244  for (int i = 0; i < 3; i++)
245  {
246  coords[i] = Array<OneD, NekDouble>(npoints);
247  }
248  f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
249 
250  // Load coordinates
251  std::string s_FieldStr;
252  for (int i = 0; i < 3; i++)
253  {
254  s_FieldStr = fieldNames[i];
255  if (f->m_session->DefinesFunction(funcName, s_FieldStr))
256  {
258  f->m_session->GetFunction(funcName, s_FieldStr);
259  ffunc->Evaluate(coords[0], coords[1], coords[2], time,
260  coords_new[i]);
261  }
262  else
263  {
264  // This coordinate is not defined, so use (x^i)' = x^i
265  Vmath::Vcopy(npoints, coords[i], 1, coords_new[i], 1);
266  }
267  }
268  // Load velocities
269  if (f->m_session->DefinesFunction(velFuncName))
270  {
271  for (int i = 0; i < 3; i++)
272  {
273  s_FieldStr = velFieldNames[i];
274  if (f->m_session->DefinesFunction(velFuncName, s_FieldStr))
275  {
277  f->m_session->GetFunction(velFuncName, s_FieldStr);
278  ffunc->Evaluate(coords[0], coords[1], coords[2], time,
279  coords_vel[i]);
280  }
281  }
282  }
283 
284  // Update mapping with coordinates
285  mapping->SetFromFunction(false);
286  mapping->UpdateMapping(time, coords_new, coords_vel);
287  }
288  else if (f->m_fieldMetaDataMap["MappingType"] == "File")
289  {
290  ASSERTL0(f->m_fieldMetaDataMap.count("FileName"),
291  "FileName parameter for Mapping missing in field file.");
292  string fileName = f->m_fieldMetaDataMap["FileName"];
293  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
294  std::vector<std::vector<NekDouble> > FieldData;
295 
296  f->FieldIOForFile(fileName)->Import(fileName, FieldDef, FieldData);
297 
298  for (int j = 0; j < spacedim; ++j)
299  {
300  int ncoeffs = f->m_exp[0]->GetNcoeffs();
301  Array<OneD, NekDouble> fieldcoeffs(ncoeffs, 0.0);
302  for (int i = 0; i < FieldData.size(); ++i)
303  {
304  f->m_exp[j]->ExtractDataToCoeffs(
305  FieldDef[i], FieldData[i], fieldNames[j], fieldcoeffs);
306  }
307  bool wavespace = f->m_exp[0]->GetWaveSpace();
308  f->m_exp[0]->SetWaveSpace(false);
309 
310  f->m_exp[0]->BwdTrans(fieldcoeffs, coords_new[j]);
311 
312  // Load coordinate velocity
313  if (std::find(FieldDef[0]->m_fields.begin(),
314  FieldDef[0]->m_fields.end(),
315  velFieldNames[j]) != FieldDef[0]->m_fields.end())
316  {
317  for (int i = 0; i < FieldData.size(); ++i)
318  {
319  f->m_exp[j]->ExtractDataToCoeffs(
320  FieldDef[i], FieldData[i], velFieldNames[j],
321  fieldcoeffs);
322  }
323  f->m_exp[0]->BwdTrans(fieldcoeffs, 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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
STL namespace.
pair< ModuleType, string > ModuleKey
static GlobalMapping::MappingSharedPtr GetMapping(FieldSharedPtr f)
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:767
double NekDouble
GLOBAL_MAPPING_EXPORT typedef boost::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:51
virtual void Process(po::variables_map &vm)
Write mesh to output file.
boost::shared_ptr< Equation > EquationSharedPtr
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
Abstract base class for processing modules.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:265
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215