Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FilterFieldConvert.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterFieldConvert.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: Base clase for filters performing operations on samples
33 // of the field.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
38 #include <boost/program_options.hpp>
39 
40 namespace Nektar
41 {
42 namespace SolverUtils
43 {
46  "FieldConvert", FilterFieldConvert::create);
47 
50  const ParamMap &pParams)
51  : Filter(pSession)
52 {
53  ParamMap::const_iterator it;
54 
55  // OutputFile
56  it = pParams.find("OutputFile");
57  if (it == pParams.end())
58  {
59  std::stringstream outname;
60  outname << m_session->GetSessionName() << ".fld";
61  m_outputFile = outname.str();
62  }
63  else
64  {
65  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
66  if ( it->second.find_last_of('.') != string::npos)
67  {
68  m_outputFile = it->second;
69  }
70  else
71  {
72  std::stringstream outname;
73  outname << it->second << ".fld";
74  m_outputFile = outname.str();
75  }
76  }
77 
78  // Restart file
79  it = pParams.find("RestartFile");
80  if (it == pParams.end())
81  {
82  m_restartFile = "";
83  }
84  else
85  {
86  ASSERTL0(it->second.length() > 0, "Missing parameter 'RestartFile'.");
87  if ( it->second.find_last_of('.') != string::npos)
88  {
89  m_restartFile = it->second;
90  }
91  else
92  {
93  std::stringstream outname;
94  outname << it->second << ".fld";
95  m_restartFile = outname.str();
96  }
97  }
98 
99  // SampleFrequency
100  it = pParams.find("SampleFrequency");
101  if (it == pParams.end())
102  {
103  m_sampleFrequency = 1;
104  }
105  else
106  {
107  LibUtilities::Equation equ(m_session, it->second);
108  m_sampleFrequency = round(equ.Evaluate());
109  }
110 
111  // OutputFrequency
112  it = pParams.find("OutputFrequency");
113  if (it == pParams.end())
114  {
115  m_outputFrequency = m_session->GetParameter("NumSteps");
116  }
117  else
118  {
119  LibUtilities::Equation equ(m_session, it->second);
120  m_outputFrequency = round(equ.Evaluate());
121  }
122 
123  m_numSamples = 0;
124  m_index = 0;
125  m_outputIndex = 0;
126 
127  //
128  // FieldConvert modules
129  //
130  m_f = boost::shared_ptr<Field>(new Field());
131  vector<string> modcmds;
132  // Process modules
133  std::stringstream moduleStream;
134  it = pParams.find("Modules");
135  if (it != pParams.end())
136  {
137  moduleStream.str(it->second);
138  }
139  while (!moduleStream.fail())
140  {
141  std::string sMod;
142  moduleStream >> sMod;
143  if (!moduleStream.fail())
144  {
145  modcmds.push_back(sMod);
146  }
147  }
148  // Output module
149  modcmds.push_back(m_outputFile);
150  // Create modules
151  CreateModules(modcmds);
152  // Strip options from m_outputFile
153  vector<string> tmp;
154  boost::split(tmp, m_outputFile, boost::is_any_of(":"));
155  m_outputFile = tmp[0];
156 }
157 
159 {
160 }
161 
164  const NekDouble &time)
165 {
166  v_FillVariablesName(pFields);
167 
168  // m_variables need to be filled by a derived class
169  m_outFields.resize(m_variables.size());
170  int nfield;
171 
172  for (int n = 0; n < m_variables.size(); ++n)
173  {
174  // if n >= pFields.num_elements() assum we have used n=0 field
175  nfield = (n < pFields.num_elements())? n: 0;
176 
177  m_outFields[n] = Array<OneD, NekDouble>(pFields[nfield]->GetNcoeffs(), 0.0);
178  }
179 
180  m_fieldMetaData["InitialTime"] = boost::lexical_cast<std::string>(time);
181 
182  // Fill some parameters of m_f
183  m_f->m_session = m_session;
184  m_f->m_graph = pFields[0]->GetGraph();
185  m_f->m_comm = m_f->m_session->GetComm();
186 
187  // Load restart file if necessary
188  if (m_restartFile != "")
189  {
190  // Load file
191  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
192  std::vector<std::vector<NekDouble> > fieldData;
193  LibUtilities::FieldMetaDataMap fieldMetaData;
196  fld->Import(m_restartFile, fieldDef, fieldData, fieldMetaData);
197 
198  // Extract fields to output
199  int nfield,k;
200  for (int j = 0; j < m_variables.size(); ++j)
201  {
202  // see if m_variables is part of pFields definition and if
203  // so use that field for extract
204  for(k = 0; k < pFields.num_elements(); ++k)
205  {
206  if(pFields[k]->GetSession()->GetVariable(k)
207  == m_variables[j])
208  {
209  nfield = k;
210  break;
211  }
212  }
213  if(k == pFields.num_elements())
214  {
215  nfield = 0;
216  }
217 
218  for (int i = 0; i < fieldData.size(); ++i)
219  {
220  pFields[nfield]->ExtractDataToCoeffs(
221  fieldDef[i],
222  fieldData[i],
223  m_variables[j],
224  m_outFields[j]);
225  }
226  }
227 
228  // Load information for numSamples
229  if (fieldMetaData.count("NumberOfFieldDumps"))
230  {
231  m_numSamples = atoi(fieldMetaData["NumberOfFieldDumps"].c_str());
232  }
233  else
234  {
235  m_numSamples = 1;
236  }
237 
238  if(fieldMetaData.count("InitialTime"))
239  {
240  m_fieldMetaData["InitialTime"] = fieldMetaData["InitialTime"];
241  }
242 
243  // Divide by scale
244  NekDouble scale = v_GetScale();
245  for (int n = 0; n < m_outFields.size(); ++n)
246  {
247  Vmath::Smul(m_outFields[n].num_elements(),
248  1.0/scale,
249  m_outFields[n],
250  1,
251  m_outFields[n],
252  1);
253  }
254  }
255 }
256 
259 {
260  int nfield = pFields.num_elements();
261  m_variables.resize(pFields.num_elements());
262  for (int n = 0; n < nfield; ++n)
263  {
264  m_variables[n] = pFields[n]->GetSession()->GetVariable(n);
265  }
266 }
267 
270  const NekDouble &time)
271 {
272  m_index++;
273  if (m_index % m_sampleFrequency > 0)
274  {
275  return;
276  }
277 
278  m_numSamples++;
279  v_ProcessSample(pFields, time);
280 
281  if (m_index % m_outputFrequency == 0)
282  {
283  m_fieldMetaData["FinalTime"] = boost::lexical_cast<std::string>(time);
284  v_PrepareOutput(pFields, time);
285  OutputField(pFields, ++m_outputIndex);
286  }
287 }
288 
291  const NekDouble &time)
292 {
293  m_fieldMetaData["FinalTime"] = boost::lexical_cast<std::string>(time);
294  v_PrepareOutput(pFields, time);
295  OutputField(pFields);
296 }
297 
300  const NekDouble &time)
301 {
302  for(int n = 0; n < pFields.num_elements(); ++n)
303  {
304  Vmath::Vcopy(m_outFields[n].num_elements(),
305  pFields[n]->GetCoeffs(),
306  1,
307  m_outFields[n],
308  1);
309  }
310 }
311 
313  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, int dump)
314 {
315  NekDouble scale = v_GetScale();
316  for (int n = 0; n < m_outFields.size(); ++n)
317  {
318  Vmath::Smul(m_outFields[n].num_elements(),
319  scale,
320  m_outFields[n],
321  1,
322  m_outFields[n],
323  1);
324  }
325 
326  CreateFields(pFields);
327 
328  // Determine new file name
329  std::stringstream outname;
330  int dot = m_outputFile.find_last_of('.');
331  string name = m_outputFile.substr(0, dot);
332  string ext = m_outputFile.substr(dot, m_outputFile.length() - dot);
333  std::string suffix = v_GetFileSuffix();
334  if (dump == -1) // final dump
335  {
336  outname << name << suffix << ext;
337  }
338  else
339  {
340  outname << name << "_" << dump << suffix << ext;
341  }
342  m_modules[m_modules.size()-1]->RegisterConfig("outfile", outname.str());
343 
344  // Prevent checking before overwritting
345  po::options_description desc("Available options");
346  desc.add_options()
347  ("forceoutput,f",
348  "Force the output to be written without any checks");
349  po::variables_map vm;
350  vm.insert(std::make_pair("forceoutput", po::variable_value()));
351  // Run field process.
352  for (int i = 0; i < m_modules.size(); ++i)
353  {
354  m_modules[i]->Process(vm);
355  cout.flush();
356  }
357 
358  // Empty m_f to save memory
359  ClearFields();
360 
361  if (dump != -1) // not final dump so rescale
362  {
363  for (int n = 0; n < m_outFields.size(); ++n)
364  {
365  Vmath::Smul(m_outFields[n].num_elements(),
366  1.0 / scale,
367  m_outFields[n],
368  1,
369  m_outFields[n],
370  1);
371  }
372  }
373 }
374 
376 {
377  return true;
378 }
379 
380 void FilterFieldConvert::CreateModules( vector<string> &modcmds)
381 {
382  for (int i = 0; i < modcmds.size(); ++i)
383  {
384  // First split each command by the colon separator.
385  vector<string> tmp1;
386  ModuleKey module;
387  int offset = 1;
388 
389  boost::split(tmp1, modcmds[i], boost::is_any_of(":"));
390 
391  if (i == modcmds.size() - 1)
392  {
393  module.first = eOutputModule;
394 
395  // If no colon detected, automatically detect mesh type from
396  // file extension. Otherwise override and use tmp1[1] as the
397  // module to load. This also allows us to pass options to
398  // input/output modules. So, for example, to override
399  // filename.xml to be read as vtk, you use:
400  //
401  // filename.xml:vtk:opt1=arg1:opt2=arg2
402  if (tmp1.size() == 1)
403  {
404  int dot = tmp1[0].find_last_of('.') + 1;
405  string ext = tmp1[0].substr(dot, tmp1[0].length() - dot);
406 
407  module.second = ext;
408  tmp1.push_back(string("outfile=") + tmp1[0]);
409  }
410  else
411  {
412  module.second = tmp1[1];
413  tmp1.push_back(string("outfile=") + tmp1[0]);
414  offset++;
415  }
416  }
417  else
418  {
419  module.first = eProcessModule;
420  module.second = tmp1[0];
421  }
422 
423  // Create modules
424  ModuleSharedPtr mod;
425  mod = GetModuleFactory().CreateInstance(module, m_f);
426  m_modules.push_back(mod);
427 
428  // Set options for this module.
429  for (int j = offset; j < tmp1.size(); ++j)
430  {
431  vector<string> tmp2;
432  boost::split(tmp2, tmp1[j], boost::is_any_of("="));
433 
434  if (tmp2.size() == 1)
435  {
436  mod->RegisterConfig(tmp2[0], "1");
437  }
438  else if (tmp2.size() == 2)
439  {
440  mod->RegisterConfig(tmp2[0], tmp2[1]);
441  }
442  else
443  {
444  cerr << "ERROR: Invalid module configuration: format is "
445  << "either :arg or :arg=val" << endl;
446  abort();
447  }
448  }
449 
450  // Ensure configuration options have been set.
451  mod->SetDefaults();
452  }
453 
454  bool RequiresEquiSpaced = false;
455  for (int i = 0; i < m_modules.size(); ++i)
456  {
457  if(m_modules[i]->GetRequireEquiSpaced())
458  {
459  RequiresEquiSpaced = true;
460  }
461  }
462  if (RequiresEquiSpaced)
463  {
464  for (int i = 0; i < m_modules.size(); ++i)
465  {
466  m_modules[i]->SetRequireEquiSpaced(true);
467  }
468  }
469 }
470 
473 {
474  m_f->m_fieldMetaDataMap = m_fieldMetaData;
475  m_f->m_fieldPts = LibUtilities::NullPtsField;
476  // Create m_f->m_exp
477  int NumHomogeneousDir = 0;
478  if (pFields[0]->GetExpType() == MultiRegions::e3DH1D)
479  {
480  NumHomogeneousDir = 1;
481  }
482  else if (pFields[0]->GetExpType() == MultiRegions::e3DH2D)
483  {
484  NumHomogeneousDir = 2;
485  }
486 
487  m_f->m_exp.resize(m_variables.size());
488  m_f->m_exp[0] = pFields[0];
489  int nfield;
490  for (int n = 0; n < m_variables.size(); ++n)
491  {
492  // if n >= pFields.num_elements() assum we have used n=0 field
493  nfield = (n < pFields.num_elements())? n: 0;
494 
495  m_f->m_exp[n] = m_f->AppendExpList(
496  NumHomogeneousDir, m_variables[0]);
497  m_f->m_exp[n]->SetWaveSpace(false);
498 
499  ASSERTL1(pFields[nfield]->GetNcoeffs() == m_outFields[n].num_elements(),
500  "pFields[nfield] does not have the "
501  "same number of coefficients as m_outFields[n]");
502 
503  m_f->m_exp[n]->ExtractCoeffsToCoeffs(pFields[nfield], m_outFields[n],
504  m_f->m_exp[n]->UpdateCoeffs());
505 
506  m_f->m_exp[n]->BwdTrans( m_f->m_exp[n]->GetCoeffs(),
507  m_f->m_exp[n]->UpdatePhys());
508  }
509  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
510  pFields[0]->GetFieldDefinitions();
511  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
512  for (int n = 0; n < m_outFields.size(); ++n)
513  {
514  for (int i = 0; i < FieldDef.size(); ++i)
515  {
516  FieldDef[i]->m_fields.push_back(m_variables[n]);
517  m_f->m_exp[n]->AppendFieldData(FieldDef[i], FieldData[i]);
518  }
519  }
520  m_f->m_fielddef = FieldDef;
521  m_f->m_data = FieldData;
522 }
523 
525 {
526  m_f->m_fieldPts = LibUtilities::NullPtsField;
527  m_f->m_exp.clear();
528  m_f->m_fielddef = std::vector<LibUtilities::FieldDefinitionsSharedPtr>();
529  m_f->m_data = std::vector<std::vector<NekDouble> > ();
530 }
531 
532 }
533 }
void OutputField(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, int dump=-1)
LibUtilities::FieldMetaDataMap m_fieldMetaData
virtual SOLVER_UTILS_EXPORT ~FilterFieldConvert()
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static std::string className
Name of the class.
SOLVER_UTILS_EXPORT FilterFieldConvert(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
void CreateFields(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields)
pair< ModuleType, string > ModuleKey
std::vector< Array< OneD, NekDouble > > m_outFields
virtual SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:54
virtual SOLVER_UTILS_EXPORT void v_PrepareOutput(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
boost::shared_ptr< Module > ModuleSharedPtr
virtual SOLVER_UTILS_EXPORT void v_FillVariablesName(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields)
virtual SOLVER_UTILS_EXPORT std::string v_GetFileSuffix()
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
NekDouble Evaluate() const
Definition: Equation.h:102
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
virtual SOLVER_UTILS_EXPORT NekDouble v_GetScale()
void CreateModules(vector< string > &modcmds)
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:309
double NekDouble
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
virtual SOLVER_UTILS_EXPORT void v_ProcessSample(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual SOLVER_UTILS_EXPORT void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:179
static boost::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:212
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
virtual SOLVER_UTILS_EXPORT void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual SOLVER_UTILS_EXPORT bool v_IsTimeDependent()
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215