Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
OutputFld.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: OutputFld.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: FLD file format output.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <set>
37 #include <string>
38 using namespace std;
39 
40 #include "OutputFld.h"
42 
43 namespace Nektar
44 {
45 namespace Utilities
46 {
47 
48 ModuleKey OutputFld::m_className[2] = {
50  ModuleKey(eOutputModule, "fld"), OutputFld::create,
51  "Writes a Fld file."),
53  ModuleKey(eOutputModule, "chk"), OutputFld::create,
54  "Writes a Fld file."),
55 };
56 
57 OutputFld::OutputFld(FieldSharedPtr f) : OutputModule(f)
58 {
59 }
60 
62 {
63 }
64 
65 void OutputFld::Process(po::variables_map &vm)
66 {
67  // Extract the output filename and extension
68  string filename = m_config["outfile"].as<string>();
69 
70  if (m_f->m_writeBndFld)
71  {
72  ModuleKey module;
73 
74  if (m_f->m_verbose)
75  {
76  if(m_f->m_comm->GetRank() == 0)
77  {
78  cout << "OutputFld: Writing boundary file(s): ";
79  for(int i = 0; i < m_f->m_bndRegionsToWrite.size(); ++i)
80  {
81  cout << m_f->m_bndRegionsToWrite[i];
82  if(i < m_f->m_bndRegionsToWrite.size()-1)
83  {
84  cout << ",";
85  }
86  }
87  cout << endl;
88  }
89  }
90 
91  // Extract data to boundaryconditions
92  if (m_f->m_fldToBnd) {
93  for (int i = 0; i < m_f->m_exp.size(); ++i)
94  {
95  m_f->m_exp[i]->FillBndCondFromField();
96  }
97  }
98 
99 
100  int nfields = m_f->m_exp.size();
102  BndExp(nfields);
103  for (int i = 0; i < nfields; ++i)
104  {
105  BndExp[i] = m_f->m_exp[i]->GetBndCondExpansions();
106  }
107 
108  // get hold of partition boundary regions so we can match it to desired
109  // region extraction
111  m_f->m_exp[0]->GetGraph());
113  bcs.GetBoundaryRegions();
114  SpatialDomains::BoundaryRegionCollection::const_iterator breg_it;
115  map<int,int> BndRegionMap;
116  int cnt =0;
117  for(breg_it = bregions.begin(); breg_it != bregions.end();
118  ++breg_it, ++cnt)
119  {
120  BndRegionMap[breg_it->first] = cnt;
121  }
122 
123  // find ending of output file and insert _b1, _b2
124  int dot = filename.find_last_of('.') + 1;
125  string ext = filename.substr(dot, filename.length() - dot);
126  string name = filename.substr(0, dot-1);
127 
128  for(int i = 0; i < m_f->m_bndRegionsToWrite.size(); ++i)
129  {
130  string outname = name + "_b"
131  + boost::lexical_cast<string>(m_f->m_bndRegionsToWrite[i])
132  + "." + ext;
133 
134  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
135  std::vector<std::vector<NekDouble> > FieldData;
136 
137  if(BndRegionMap.count(m_f->m_bndRegionsToWrite[i]) == 1)
138  {
139  int Border = BndRegionMap[m_f->m_bndRegionsToWrite[i]];
140 
141  FieldDef = BndExp[0][Border]->GetFieldDefinitions();
142  FieldData.resize(FieldDef.size());
143 
144  for (int j = 0; j < nfields; ++j)
145  {
146  for (int k = 0; k < FieldDef.size(); ++k)
147  {
148  BndExp[j][Border]->AppendFieldData(FieldDef[k],
149  FieldData[k]);
150 
151  if (m_f->m_fielddef.size() > 0)
152  {
153  FieldDef[k]->m_fields.push_back(
154  m_f->m_fielddef[0]->m_fields[j]);
155  }
156  else
157  {
158  FieldDef[k]->m_fields.push_back(
159  m_f->m_session->GetVariable(j));
160  }
161  }
162  }
163 
164  if(m_f->m_addNormals)
165  {
166  int normdim = m_f->m_graph->GetMeshDimension();
167  string normstr[3] = {"Norm_x","Norm_y","Norm_z"};
168 
169  // Add normal information
171  Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
172 
173  m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
174  BoundarytoTraceID);
175 
176  // determine offset of this Bnd Expansion Border
177  int cnt = 0;
178  for(int n = 0; n < Border; ++n)
179  {
180  cnt += BndExp[0][n]->GetExpSize();
181  }
182 
183  Array<OneD, NekDouble> tmp_array;
184  Array<OneD, Array<OneD, NekDouble> > NormPhys(normdim);
185 
186  for(int j = 0; j < normdim; ++j)
187  {
188  NormPhys[j] = Array<OneD, NekDouble>(BndExp[0][Border]->GetTotPoints(),0.0);
189  }
190 
191  // setup phys arrays of normals
192  for(int j=0; j < BndExp[0][Border]->GetExpSize(); ++j,++cnt)
193  {
194  int elmtid = BoundarytoElmtID[cnt];
195 
196  elmt = m_f->m_exp[0]->GetExp(elmtid);
197 
198  //identify boundary of element looking at.
199  int boundary = BoundarytoTraceID[cnt];
200 
201  // Dimension specific part
202  switch(normdim)
203  {
204  case 2:
205  {
206  // Get edge 1D expansion from element expansion
208  bc = boost::dynamic_pointer_cast
210  (BndExp[0][Border]->GetExp(j));
211  // Get edge normals
213  normals = elmt->GetEdgeNormal(boundary);
214 
215  for(int k = 0; k < normdim; ++k)
216  {
217  Vmath::Vcopy(bc->GetTotPoints(),
218  normals[k], 1,
219  tmp_array = NormPhys[k]+
220  BndExp[0][Border]->
221  GetPhys_Offset(j), 1);
222  }
223  }
224  break;
225 
226  case 3:
227  {
228  // Get face 2D expansion from element expansion
230  bc = boost::dynamic_pointer_cast
232  (BndExp[0][Border]->GetExp(j));
233  //Get face normals
235  normals = elmt->GetFaceNormal(boundary);
236 
237  for(int k = 0; k < normdim; ++k)
238  {
239  Vmath::Vcopy(bc->GetTotPoints(),
240  normals[k], 1,
241  tmp_array = NormPhys[k]+
242  BndExp[0][Border]->
243  GetPhys_Offset(j), 1);
244  }
245  }
246  break;
247 
248  default:
249  ASSERTL0(false, "Addnormals requires expdim >=2.");
250  break;
251  }
252  }
253 
254  // add normal coefficients to list to be dumped
255  for (int j = 0; j < normdim; ++j)
256  {
257  BndExp[0][Border]->FwdTrans( NormPhys[j],
258  BndExp[0][Border]->UpdateCoeffs());
259 
260  for (int k = 0; k < FieldDef.size(); ++k)
261  {
262  BndExp[0][Border]->AppendFieldData(FieldDef[k],
263  FieldData[k]);
264  FieldDef[k]->m_fields.push_back(normstr[j]);
265  }
266  }
267  }
268 
269  // output error for regression checking.
270  if (vm.count("error"))
271  {
272  int rank = m_f->m_session->GetComm()->GetRank();
273 
274  for (int j = 0; j < nfields; ++j)
275  {
276  BndExp[j][Border]->BwdTrans(BndExp[j][Border]->GetCoeffs(),
277  BndExp[j][Border]->UpdatePhys());
278 
279  //Note currently these calls will
280  //hange since not all partitions will
281  //call error.
282  NekDouble l2err = BndExp[j][Border]
283  ->L2(BndExp[j][Border]->GetPhys());
284 
285  NekDouble linferr = BndExp[j][Border]
286  ->Linf(BndExp[j][Border]->GetPhys());
287 
288  if (rank == 0)
289  {
290  cout << "L 2 error (variable "
291  << FieldDef[0]->m_fields[j]
292  << ") : " << l2err << endl;
293 
294  cout << "L inf error (variable "
295  << FieldDef[0]->m_fields[j]
296  << ") : " << linferr << endl;
297  }
298  }
299  }
300  }
301 
302  m_f->m_fld->Write(outname, FieldDef, FieldData,
303  m_f->m_fieldMetaDataMap);
304 
305  }
306  }
307  else
308  {
309  if (m_f->m_verbose)
310  {
311  if(m_f->m_comm->GetRank() == 0)
312  {
313  cout << "OutputFld: Writing file..." << endl;
314  }
315  }
316 
317  fs::path writefile(filename);
318  int writefld = 1;
319  if(fs::exists(writefile)&&(vm.count("forceoutput") == 0))
320  {
321  LibUtilities::CommSharedPtr comm = m_f->m_session->GetComm();
322  int rank = comm->GetRank();
323  writefld = 0; // set to zero for reduce all to be correct.
324 
325  if(rank == 0)
326  {
327  string answer;
328  cout << "Did you wish to overwrite " << filename << " (y/n)? ";
329  getline(cin,answer);
330  if(answer.compare("y") == 0)
331  {
332  writefld = 1;
333  }
334  else
335  {
336  cout << "Not writing file " << filename << " because it already exists" << endl;
337  }
338  }
339 
340  comm->AllReduce(writefld,LibUtilities::ReduceSum);
341 
342  }
343 
344  if(writefld)
345  {
346  m_f->m_fld->Write(filename, m_f->m_fielddef, m_f->m_data,
347  m_f->m_fieldMetaDataMap);
348  }
349 
350  // output error for regression checking.
351  if (vm.count("error"))
352  {
353  int rank = m_f->m_session->GetComm()->GetRank();
354 
355  for (int j = 0; j < m_f->m_exp.size(); ++j)
356  {
357  if (m_f->m_exp[j]->GetPhysState() == false)
358  {
359  m_f->m_exp[j]->BwdTrans(
360  m_f->m_exp[j]->GetCoeffs(),
361  m_f->m_exp[j]->UpdatePhys());
362  }
363 
364  NekDouble l2err = m_f->m_exp[j]->L2(
365  m_f->m_exp[j]->GetPhys());
366 
367  NekDouble linferr = m_f->m_exp[j]->Linf(
368  m_f->m_exp[j]->GetPhys());
369  if (rank == 0)
370  {
371  cout << "L 2 error (variable "
372  << m_f->m_fielddef[0]->m_fields[j]
373  << ") : " << l2err << endl;
374 
375  cout << "L inf error (variable "
376  << m_f->m_fielddef[0]->m_fields[j]
377  << ") : " << linferr << endl;
378  }
379  }
380  }
381  }
382 }
383 
384 }
385 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
pair< ModuleType, string > ModuleKey
Abstract base class for output modules.
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
FieldSharedPtr m_f
Field object.
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:206
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:227
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215