Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FieldIO.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: FieldIO.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: I/O routines relating to Fields
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/format.hpp>
37 #include <boost/date_time/posix_time/posix_time.hpp>
38 #include <boost/date_time/posix_time/posix_time_io.hpp>
39 #include <boost/asio/ip/host_name.hpp>
40 
41 #include <LibUtilities/BasicUtils/FieldIO.h>
42 #include <LibUtilities/BasicUtils/FileSystem.h>
43 #include <LibUtilities/BasicConst/GitRevision.h>
44 
45 #include "zlib.h"
46 #include <set>
47 
48 #ifdef NEKTAR_USE_MPI
49 #include <mpi.h>
50 #endif
51 
52 // Buffer size for zlib compression/decompression
53 #define CHUNK 16384
54 
55 #ifndef NEKTAR_VERSION
56 #define NEKTAR_VERSION "Unknown"
57 #endif
58 
59 namespace ptime = boost::posix_time;
60 namespace ip = boost::asio::ip;
61 namespace berrc = boost::system::errc;
62 
63 namespace Nektar
64 {
65  namespace LibUtilities
66  {
67  /**
68  * This function allows for data to be written to an FLD file when a
69  * session and/or communicator is not instantiated. Typically used in
70  * utilities which do not take XML input and operate in serial only.
71  */
72  void Write( const std::string &outFile,
73  std::vector<FieldDefinitionsSharedPtr> &fielddefs,
74  std::vector<std::vector<NekDouble> > &fielddata,
75  const FieldMetaDataMap &fieldinfomap)
76  {
77 #ifdef NEKTAR_USE_MPI
78  int size;
79  int init;
80  MPI_Initialized(&init);
81 
82  // If MPI has been initialised we can check the number of processes
83  // and, if > 1, tell the user he should not be running this
84  // function in parallel. If it is not initialised, we do not
85  // initialise it here, and assume the user knows what they are
86  // doing.
87  if (init)
88  {
89  MPI_Comm_size( MPI_COMM_WORLD, &size );
90  ASSERTL0(size == 1,
91  "This static function is not available in parallel. Please"
92  "instantiate a FieldIO object for parallel use.");
93  }
94 #endif
95  CommSharedPtr c = GetCommFactory().CreateInstance("Serial", 0, 0);
96  FieldIO f(c);
97  f.Write(outFile, fielddefs, fielddata, fieldinfomap);
98  }
99 
100 
101  /**
102  * This function allows for data to be imported from an FLD file when
103  * a session and/or communicator is not instantiated. Typically used in
104  * utilities which only operate in serial.
105  */
106  LIB_UTILITIES_EXPORT void Import(
107  const std::string& infilename,
108  std::vector<FieldDefinitionsSharedPtr> &fielddefs,
109  std::vector<std::vector<NekDouble> > &fielddata,
110  FieldMetaDataMap &fieldinfomap,
111  const Array<OneD, int> ElementiDs)
112  {
113 #ifdef NEKTAR_USE_MPI
114  int size;
115  int init;
116  MPI_Initialized(&init);
117 
118  // If MPI has been initialised we can check the number of processes
119  // and, if > 1, tell the user he should not be running this
120  // function in parallel. If it is not initialised, we do not
121  // initialise it here, and assume the user knows what they are
122  // doing.
123  if (init)
124  {
125  MPI_Comm_size( MPI_COMM_WORLD, &size );
126  ASSERTL0(size == 1,
127  "This static function is not available in parallel. Please"
128  "instantiate a FieldIO object for parallel use.");
129  }
130 #endif
131  CommSharedPtr c = GetCommFactory().CreateInstance("Serial", 0, 0);
132  FieldIO f(c);
133  f.Import(infilename, fielddefs, fielddata, fieldinfomap, ElementiDs);
134  }
135 
136 
137  /**
138  *
139  */
140  FieldIO::FieldIO(
141  LibUtilities::CommSharedPtr pComm)
142  : m_comm(pComm)
143  {
144  }
145 
146 
147  /**
148  *
149  */
150  void FieldIO::Write(const std::string &outFile,
151  std::vector<FieldDefinitionsSharedPtr> &fielddefs,
152  std::vector<std::vector<NekDouble> > &fielddata,
153  const FieldMetaDataMap &fieldmetadatamap)
154  {
155  // Check everything seems sensible
156  ASSERTL1(fielddefs.size() == fielddata.size(),
157  "Length of fielddefs and fielddata incompatible");
158  for (int f = 0; f < fielddefs.size(); ++f)
159  {
160  ASSERTL1(fielddata[f].size() > 0,
161  "Fielddata vector must contain at least one value.");
162 
163  ASSERTL1(fielddata[f].size() ==
164  fielddefs[f]->m_fields.size() *
165  CheckFieldDefinition(fielddefs[f]),
166  "Invalid size of fielddata vector.");
167  }
168 
169  // Prepare to write out data. In parallel, we must create directory
170  // and determine the full pathname to the file to write out.
171  // Any existing file/directory which is in the way is removed.
172  std::string filename = SetUpOutput(outFile, fielddefs, fieldmetadatamap);
173 
174  // Create the file (partition)
175  TiXmlDocument doc;
176  TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
177  doc.LinkEndChild(decl);
178 
179  TiXmlElement * root = new TiXmlElement("NEKTAR");
180  doc.LinkEndChild(root);
181 
182  AddInfoTag(root,fieldmetadatamap);
183 
184  for (int f = 0; f < fielddefs.size(); ++f)
185  {
186  //---------------------------------------------
187  // Write ELEMENTS
188  TiXmlElement * elemTag = new TiXmlElement("ELEMENTS");
189  root->LinkEndChild(elemTag);
190 
191  // Write FIELDS
192  std::string fieldsString;
193  {
194  std::stringstream fieldsStringStream;
195  bool first = true;
196  for (std::vector<int>::size_type i = 0; i
197  < fielddefs[f]->m_fields.size(); i++)
198  {
199  if (!first)
200  fieldsStringStream << ",";
201  fieldsStringStream << fielddefs[f]->m_fields[i];
202  first = false;
203  }
204  fieldsString = fieldsStringStream.str();
205  }
206  elemTag->SetAttribute("FIELDS", fieldsString);
207 
208  // Write SHAPE
209  std::string shapeString;
210  {
211  std::stringstream shapeStringStream;
212  shapeStringStream << ShapeTypeMap[fielddefs[f]->m_shapeType];
213  if(fielddefs[f]->m_numHomogeneousDir == 1)
214  {
215  shapeStringStream << "-HomogenousExp1D";
216  }
217  else if (fielddefs[f]->m_numHomogeneousDir == 2)
218  {
219  shapeStringStream << "-HomogenousExp2D";
220  }
221 
222  shapeString = shapeStringStream.str();
223  }
224  elemTag->SetAttribute("SHAPE", shapeString);
225 
226  // Write BASIS
227  std::string basisString;
228  {
229  std::stringstream basisStringStream;
230  bool first = true;
231  for (std::vector<BasisType>::size_type i = 0; i < fielddefs[f]->m_basis.size(); i++)
232  {
233  if (!first)
234  basisStringStream << ",";
235  basisStringStream
236  << BasisTypeMap[fielddefs[f]->m_basis[i]];
237  first = false;
238  }
239  basisString = basisStringStream.str();
240  }
241  elemTag->SetAttribute("BASIS", basisString);
242 
243  // Write homogeneuous length details
244  if(fielddefs[f]->m_numHomogeneousDir)
245  {
246  std::string homoLenString;
247  {
248  std::stringstream homoLenStringStream;
249  bool first = true;
250  for (int i = 0; i < fielddefs[f]->m_numHomogeneousDir; ++i)
251  {
252  if (!first)
253  homoLenStringStream << ",";
254  homoLenStringStream
255  << fielddefs[f]->m_homogeneousLengths[i];
256  first = false;
257  }
258  homoLenString = homoLenStringStream.str();
259  }
260  elemTag->SetAttribute("HOMOGENEOUSLENGTHS", homoLenString);
261  }
262 
263  // Write homogeneuous planes/lines details
264  if(fielddefs[f]->m_numHomogeneousDir)
265  {
266  if(fielddefs[f]->m_homogeneousYIDs.size() > 0)
267  {
268  std::string homoYIDsString;
269  {
270  std::stringstream homoYIDsStringStream;
271  bool first = true;
272  for(int i = 0; i < fielddefs[f]->m_homogeneousYIDs.size(); i++)
273  {
274  if (!first)
275  homoYIDsStringStream << ",";
276  homoYIDsStringStream << fielddefs[f]->m_homogeneousYIDs[i];
277  first = false;
278  }
279  homoYIDsString = homoYIDsStringStream.str();
280  }
281  elemTag->SetAttribute("HOMOGENEOUSYIDS", homoYIDsString);
282  }
283 
284  if(fielddefs[f]->m_homogeneousZIDs.size() > 0)
285  {
286  std::string homoZIDsString;
287  {
288  std::stringstream homoZIDsStringStream;
289  bool first = true;
290  for(int i = 0; i < fielddefs[f]->m_homogeneousZIDs.size(); i++)
291  {
292  if (!first)
293  homoZIDsStringStream << ",";
294  homoZIDsStringStream << fielddefs[f]->m_homogeneousZIDs[i];
295  first = false;
296  }
297  homoZIDsString = homoZIDsStringStream.str();
298  }
299  elemTag->SetAttribute("HOMOGENEOUSZIDS", homoZIDsString);
300  }
301  }
302 
303  // Write NUMMODESPERDIR
304  std::string numModesString;
305  {
306  std::stringstream numModesStringStream;
307 
308  if (fielddefs[f]->m_uniOrder)
309  {
310  numModesStringStream << "UNIORDER:";
311  // Just dump single definition
312  bool first = true;
313  for (std::vector<int>::size_type i = 0; i
314  < fielddefs[f]->m_basis.size(); i++)
315  {
316  if (!first)
317  numModesStringStream << ",";
318  numModesStringStream << fielddefs[f]->m_numModes[i];
319  first = false;
320  }
321  }
322  else
323  {
324  numModesStringStream << "MIXORDER:";
325  bool first = true;
326  for (std::vector<int>::size_type i = 0; i
327  < fielddefs[f]->m_numModes.size(); i++)
328  {
329  if (!first)
330  numModesStringStream << ",";
331  numModesStringStream << fielddefs[f]->m_numModes[i];
332  first = false;
333  }
334  }
335 
336  numModesString = numModesStringStream.str();
337  }
338  elemTag->SetAttribute("NUMMODESPERDIR", numModesString);
339 
340  // Write ID
341  // Should ideally look at ways of compressing this stream
342  // if just sequential;
343  std::string idString;
344  {
345  std::stringstream idStringStream;
346  GenerateSeqString(fielddefs[f]->m_elementIDs,idString);
347  }
348  elemTag->SetAttribute("ID", idString);
349 
350  std::string compressedDataString;
351  ASSERTL0(Z_OK == Deflate(fielddata[f], compressedDataString),
352  "Failed to compress field data.");
353 
354  // If the string length is not divisible by 3,
355  // pad it. There is a bug in transform_width
356  // that will make it reference past the end
357  // and crash.
358  switch (compressedDataString.length() % 3)
359  {
360  case 1:
361  compressedDataString += '\0';
362  case 2:
363  compressedDataString += '\0';
364  break;
365  }
366 
367  // Convert from binary to base64.
368  typedef boost::archive::iterators::base64_from_binary<
369  boost::archive::iterators::transform_width<
370  std::string::const_iterator, 6, 8> > base64_t;
371  std::string base64string(base64_t(compressedDataString.begin()),
372  base64_t(compressedDataString.end()));
373  elemTag->LinkEndChild(new TiXmlText(base64string));
374 
375  }
376  doc.SaveFile(filename);
377  }
378 
379 
380  /**
381  *
382  */
383  void FieldIO::Import(const std::string& infilename,
384  std::vector<FieldDefinitionsSharedPtr> &fielddefs,
385  std::vector<std::vector<NekDouble> > &fielddata,
386  FieldMetaDataMap &fieldmetadatamap,
387  const Array<OneD, int> ElementIDs)
388  {
389 
390  std::string infile = infilename;
391 
392  fs::path pinfilename(infilename);
393 
394  if(fs::is_directory(pinfilename)) // check to see that infile is a directory
395  {
396  fs::path infofile("Info.xml");
397  fs::path fullpath = pinfilename / infofile;
398  infile = PortablePath(fullpath);
399 
400  std::vector<std::string> filenames;
401  std::vector<std::vector<unsigned int> > elementIDs_OnPartitions;
402 
403 
404  ImportMultiFldFileIDs(infile,filenames, elementIDs_OnPartitions,
405  fieldmetadatamap);
406 
407  if(ElementIDs == NullInt1DArray) //load all fields
408  {
409  for(int i = 0; i < filenames.size(); ++i)
410  {
411  fs::path pfilename(filenames[i]);
412  fullpath = pinfilename / pfilename;
413  string fname = PortablePath(fullpath);
414 
415  TiXmlDocument doc1(fname);
416  bool loadOkay1 = doc1.LoadFile();
417 
418  std::stringstream errstr;
419  errstr << "Unable to load file: " << fname << std::endl;
420  errstr << "Reason: " << doc1.ErrorDesc() << std::endl;
421  errstr << "Position: Line " << doc1.ErrorRow() << ", Column " << doc1.ErrorCol() << std::endl;
422  ASSERTL0(loadOkay1, errstr.str());
423 
424  ImportFieldDefs(doc1, fielddefs, false);
425  if(fielddata != NullVectorNekDoubleVector)
426  {
427  ImportFieldData(doc1, fielddefs, fielddata);
428  }
429  }
430 
431  }
432  else // only load relevant partitions
433  {
434  int i,j;
435  map<int,vector<int> > FileIDs;
436  map<int,vector<int> >::iterator it;
437  set<int> LoadFile;
438 
439  for(i = 0; i < elementIDs_OnPartitions.size(); ++i)
440  {
441  for(j = 0; j < elementIDs_OnPartitions[i].size(); ++j)
442  {
443  FileIDs[elementIDs_OnPartitions[i][j]].push_back(i);
444  }
445  }
446 
447  for(i = 0; i < ElementIDs.num_elements(); ++i)
448  {
449  it = FileIDs.find(ElementIDs[i]);
450  if (it != FileIDs.end())
451  {
452  for (j = 0; j < it->second.size(); ++j)
453  {
454  LoadFile.insert(it->second[j]);
455  }
456  }
457  }
458 
459  set<int>::iterator iter;
460  for(iter = LoadFile.begin(); iter != LoadFile.end(); ++iter)
461  {
462  fs::path pfilename(filenames[*iter]);
463  fullpath = pinfilename / pfilename;
464  string fname = PortablePath(fullpath);
465  TiXmlDocument doc1(fname);
466  bool loadOkay1 = doc1.LoadFile();
467 
468  std::stringstream errstr;
469  errstr << "Unable to load file: " << fname << std::endl;
470  errstr << "Reason: " << doc1.ErrorDesc() << std::endl;
471  errstr << "Position: Line " << doc1.ErrorRow() << ", Column " << doc1.ErrorCol() << std::endl;
472  ASSERTL0(loadOkay1, errstr.str());
473 
474  ImportFieldDefs(doc1, fielddefs, false);
475  if(fielddata != NullVectorNekDoubleVector)
476  {
477  ImportFieldData(doc1, fielddefs, fielddata);
478  }
479  }
480  }
481  }
482  else // serial format case
483  {
484 
485  TiXmlDocument doc(infile);
486  bool loadOkay = doc.LoadFile();
487 
488  std::stringstream errstr;
489  errstr << "Unable to load file: " << infile << std::endl;
490  errstr << "Reason: " << doc.ErrorDesc() << std::endl;
491  errstr << "Position: Line " << doc.ErrorRow() << ", Column " <<
492  doc.ErrorCol() << std::endl;
493  ASSERTL0(loadOkay, errstr.str());
494 
495  ImportFieldMetaData(doc,fieldmetadatamap);
496  ImportFieldDefs(doc, fielddefs, false);
497  if(fielddata != NullVectorNekDoubleVector)
498  {
499  ImportFieldData(doc, fielddefs, fielddata);
500  }
501  }
502  }
503 
504 
505  /**
506  *
507  */
508  void FieldIO::WriteMultiFldFileIDs(const std::string &outFile,
509  const std::vector<std::string> fileNames,
510  std::vector<std::vector<unsigned int> > &elementList,
511  const FieldMetaDataMap &fieldmetadatamap)
512  {
513  TiXmlDocument doc;
514  TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
515  doc.LinkEndChild(decl);
516 
517  ASSERTL0(fileNames.size() == elementList.size(),"Outfile names and list of elements ids does not match");
518 
519  TiXmlElement * root = new TiXmlElement("NEKTAR");
520  doc.LinkEndChild(root);
521 
522  AddInfoTag(root,fieldmetadatamap);
523 
524  for (int t = 0; t < fileNames.size(); ++t)
525  {
526  if(elementList[t].size())
527  {
528  TiXmlElement * elemIDs = new TiXmlElement("Partition");
529  root->LinkEndChild(elemIDs);
530 
531  elemIDs->SetAttribute("FileName",fileNames[t]);
532 
533  string IDstring;
534 
535  GenerateSeqString(elementList[t],IDstring);
536 
537  elemIDs->LinkEndChild(new TiXmlText(IDstring));
538  }
539  }
540 
541  doc.SaveFile(outFile);
542  }
543 
544 
545  /**
546  *
547  */
548  void FieldIO::ImportMultiFldFileIDs(const std::string &inFile,
549  std::vector<std::string> &fileNames,
550  std::vector<std::vector<unsigned int> > &elementList,
551  FieldMetaDataMap &fieldmetadatamap)
552  {
553  TiXmlDocument doc(inFile);
554  bool loadOkay = doc.LoadFile();
555 
556 
557  std::stringstream errstr;
558  errstr << "Unable to load file: " << inFile<< std::endl;
559  errstr << "Reason: " << doc.ErrorDesc() << std::endl;
560  errstr << "Position: Line " << doc.ErrorRow() << ", Column " << doc.ErrorCol() << std::endl;
561  ASSERTL0(loadOkay, errstr.str());
562 
563  // Handle on XML document
564  TiXmlHandle docHandle(&doc);
565 
566  // Retrieve main NEKTAR tag - XML specification states one
567  // top-level element tag per file.
568  TiXmlElement* master = doc.FirstChildElement("NEKTAR");
569  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
570 
571  // Partition element tag name
572  std::string strPartition = "Partition";
573 
574  // First attempt to get the first Partition element
575  TiXmlElement* fldfileIDs = master->FirstChildElement(strPartition.c_str());
576  if (!fldfileIDs)
577  {
578  // If this files try previous name
579  strPartition = "MultipleFldFiles";
580  fldfileIDs = master->FirstChildElement("MultipleFldFiles");
581  }
582  ASSERTL0(fldfileIDs,
583  "Unable to find 'Partition' or 'MultipleFldFiles' tag "
584  "within nektar tag.");
585 
586  while (fldfileIDs)
587  {
588  // Read file name of partition file
589  const char *attr = fldfileIDs->Attribute("FileName");
590  ASSERTL0(attr, "'FileName' not provided as an attribute of '"
591  + strPartition + "' tag.");
592  fileNames.push_back(std::string(attr));
593 
594  const char* elementIDs = fldfileIDs->GetText();
595  ASSERTL0(elementIDs, "Element IDs not specified.");
596 
597  std::string elementIDsStr(elementIDs);
598 
599  std::vector<unsigned int> idvec;
600  ParseUtils::GenerateSeqVector(elementIDsStr.c_str(),idvec);
601 
602  elementList.push_back(idvec);
603 
604  fldfileIDs = fldfileIDs->NextSiblingElement(strPartition.c_str());
605  }
606  }
607 
608  void FieldIO::ImportFieldMetaData(std::string filename,
609  FieldMetaDataMap &fieldmetadatamap)
610  {
611  TiXmlDocument doc(filename);
612  bool loadOkay = doc.LoadFile();
613 
614  std::stringstream errstr;
615  errstr << "Unable to load file: " << filename << std::endl;
616  errstr << "Reason: " << doc.ErrorDesc() << std::endl;
617  errstr << "Position: Line " << doc.ErrorRow() << ", Column " << doc.ErrorCol() << std::endl;
618  ASSERTL0(loadOkay, errstr.str());
619 
620  ImportFieldMetaData(doc,fieldmetadatamap);
621  }
622 
623 
624  void FieldIO::ImportFieldMetaData(TiXmlDocument &doc,
625  FieldMetaDataMap &fieldmetadatamap)
626  {
627 
628  TiXmlHandle docHandle(&doc);
629  TiXmlElement* master = 0; // Master tag within which all data is contained.
630  TiXmlElement* metadata = 0;
631 
632  master = doc.FirstChildElement("NEKTAR");
633  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
634  std::string strLoop = "NEKTAR";
635 
636  // Retain original metadata structure for backwards compatibility
637  // TODO: Remove old metadata format
638  metadata = master->FirstChildElement("FIELDMETADATA");
639  if(metadata)
640  {
641  TiXmlElement *param = metadata->FirstChildElement("P");
642 
643  while (param)
644  {
645  TiXmlAttribute *paramAttr = param->FirstAttribute();
646  std::string attrName(paramAttr->Name());
647  std::string paramString;
648 
649  if(attrName == "PARAM")
650  {
651  paramString.insert(0,paramAttr->Value());
652  }
653  else
654  {
655  ASSERTL0(false,"PARAM not provided as an attribute in FIELDMETADATA section");
656  }
657 
658  // Now read body of param
659  std::string paramBodyStr;
660 
661  TiXmlNode *paramBody = param->FirstChild();
662 
663  paramBodyStr += paramBody->ToText()->Value();
664 
665  fieldmetadatamap[paramString] = paramBodyStr;
666  param = param->NextSiblingElement("P");
667  }
668  }
669 
670  // New metadata format
671  metadata = master->FirstChildElement("Metadata");
672  if(metadata)
673  {
674  TiXmlElement *param = metadata->FirstChildElement();
675 
676  while (param)
677  {
678  std::string paramString = param->Value();
679  if (paramString != "Provenance")
680  {
681  // Now read body of param
682  TiXmlNode *paramBody = param->FirstChild();
683  std::string paramBodyStr = paramBody->ToText()->Value();
684 
685  fieldmetadatamap[paramString] = paramBodyStr;
686  }
687  param = param->NextSiblingElement();
688  }
689  }
690 
691  }
692 
693 
694  /**
695  * The bool decides if the FieldDefs are in <EXPANSIONS> or in <NEKTAR>.
696  */
697  void FieldIO::ImportFieldDefs(TiXmlDocument &doc, std::vector<FieldDefinitionsSharedPtr> &fielddefs,
698  bool expChild)
699  {
700  TiXmlHandle docHandle(&doc);
701  TiXmlElement* master = NULL; // Master tag within which all data is contained.
702 
703  master = doc.FirstChildElement("NEKTAR");
704  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
705  std::string strLoop = "NEKTAR";
706  TiXmlElement* loopXml = master;
707 
708  TiXmlElement *expansionTypes;
709  if(expChild)
710  {
711  expansionTypes = master->FirstChildElement("EXPANSIONS");
712  ASSERTL0(expansionTypes, "Unable to find EXPANSIONS tag in file.");
713  loopXml = expansionTypes;
714  strLoop = "EXPANSIONS";
715  }
716 
717  // Loop through all nektar tags, finding all of the element tags.
718  while (loopXml)
719  {
720  TiXmlElement* element = loopXml->FirstChildElement("ELEMENTS");
721  ASSERTL0(element, "Unable to find ELEMENTS tag within nektar tag.");
722 
723  while (element)
724  {
725  // Extract the attributes.
726  std::string idString;
727  std::string shapeString;
728  std::string basisString;
729  std::string homoLengthsString;
730  std::string homoZIDsString;
731  std::string homoYIDsString;
732  std::string numModesString;
733  std::string numPointsString;
734  std::string fieldsString;
735  std::string pointsString;
736  bool pointDef = false;
737  bool numPointDef = false;
738  TiXmlAttribute *attr = element->FirstAttribute();
739  while (attr)
740  {
741  std::string attrName(attr->Name());
742  if (attrName == "FIELDS")
743  {
744  fieldsString.insert(0, attr->Value());
745  }
746  else if (attrName == "SHAPE")
747  {
748  shapeString.insert(0, attr->Value());
749  }
750  else if (attrName == "BASIS")
751  {
752  basisString.insert(0, attr->Value());
753  }
754  else if (attrName == "HOMOGENEOUSLENGTHS")
755  {
756  homoLengthsString.insert(0,attr->Value());
757  }
758  else if (attrName == "HOMOGENEOUSZIDS")
759  {
760  homoZIDsString.insert(0,attr->Value());
761  }
762  else if (attrName == "HOMOGENEOUSYIDS")
763  {
764  homoYIDsString.insert(0,attr->Value());
765  }
766  else if (attrName == "NUMMODESPERDIR")
767  {
768  numModesString.insert(0, attr->Value());
769  }
770  else if (attrName == "ID")
771  {
772  idString.insert(0, attr->Value());
773  }
774  else if (attrName == "POINTSTYPE")
775  {
776  pointsString.insert(0, attr->Value());
777  pointDef = true;
778  }
779  else if (attrName == "NUMPOINTSPERDIR")
780  {
781  numPointsString.insert(0, attr->Value());
782  numPointDef = true;
783  }
784  else
785  {
786  std::string errstr("Unknown attribute: ");
787  errstr += attrName;
788  ASSERTL1(false, errstr.c_str());
789  }
790 
791  // Get the next attribute.
792  attr = attr->Next();
793  }
794 
795  // Check to see if homogeneous expansion and if so
796  // strip down the shapeString definition
797  int numHomoDir = 0;
798  size_t loc;
799  //---> This finds the first location of 'n'!
800  if((loc = shapeString.find_first_of("-"))!=string::npos)
801  {
802  if(shapeString.find("Exp1D")!=string::npos)
803  {
804  numHomoDir = 1;
805  }
806  else // HomogeneousExp1D
807  {
808  numHomoDir = 2;
809  }
810 
811  shapeString.erase(loc,shapeString.length());
812  }
813 
814  // Reconstruct the fielddefs.
815  std::vector<unsigned int> elementIds;
816  {
817  bool valid = ParseUtils::GenerateSeqVector(idString.c_str(), elementIds);
818  ASSERTL0(valid, "Unable to correctly parse the element ids.");
819  }
820 
821  // Get the geometrical shape
822  ShapeType shape;
823  bool valid = false;
824  for (unsigned int j = 0; j < SIZE_ShapeType; j++)
825  {
826  if (ShapeTypeMap[j] == shapeString)
827  {
828  shape = (ShapeType) j;
829  valid = true;
830  break;
831  }
832  }
833 
834  ASSERTL0(valid, std::string("Unable to correctly parse the shape type: ").append(shapeString).c_str());
835 
836  // Get the basis
837  std::vector<std::string> basisStrings;
838  std::vector<BasisType> basis;
839  valid = ParseUtils::GenerateOrderedStringVector(basisString.c_str(), basisStrings);
840  ASSERTL0(valid, "Unable to correctly parse the basis types.");
841  for (std::vector<std::string>::size_type i = 0; i < basisStrings.size(); i++)
842  {
843  valid = false;
844  for (unsigned int j = 0; j < SIZE_BasisType; j++)
845  {
846  if (BasisTypeMap[j] == basisStrings[i])
847  {
848  basis.push_back((BasisType) j);
849  valid = true;
850  break;
851  }
852  }
853  ASSERTL0(valid, std::string("Unable to correctly parse the basis type: ").append(basisStrings[i]).c_str());
854  }
855 
856  // Get homoLengths
857  std::vector<NekDouble> homoLengths;
858  if(numHomoDir)
859  {
860  valid = ParseUtils::GenerateUnOrderedVector(homoLengthsString.c_str(), homoLengths);
861  ASSERTL0(valid, "Unable to correctly parse the number of homogeneous lengths.");
862  }
863 
864  // Get Homogeneous points IDs
865  std::vector<unsigned int> homoZIDs;
866  std::vector<unsigned int> homoYIDs;
867 
868  if(numHomoDir == 1)
869  {
870  valid = ParseUtils::GenerateSeqVector(homoZIDsString.c_str(), homoZIDs);
871  ASSERTL0(valid, "Unable to correctly parse homogeneous planes IDs.");
872  }
873 
874  if(numHomoDir == 2)
875  {
876  valid = ParseUtils::GenerateSeqVector(homoZIDsString.c_str(), homoZIDs);
877  ASSERTL0(valid, "Unable to correctly parse homogeneous lines IDs in z-direction.");
878  valid = ParseUtils::GenerateSeqVector(homoYIDsString.c_str(), homoYIDs);
879  ASSERTL0(valid, "Unable to correctly parse homogeneous lines IDs in y-direction.");
880  }
881 
882 
883  // Get points type
884  std::vector<PointsType> points;
885 
886  if(pointDef)
887  {
888  std::vector<std::string> pointsStrings;
889  valid = ParseUtils::GenerateOrderedStringVector(pointsString.c_str(), pointsStrings);
890  ASSERTL0(valid, "Unable to correctly parse the points types.");
891  for (std::vector<std::string>::size_type i = 0; i < pointsStrings.size(); i++)
892  {
893  valid = false;
894  for (unsigned int j = 0; j < SIZE_PointsType; j++)
895  {
896  if (kPointsTypeStr[j] == pointsStrings[i])
897  {
898  points.push_back((PointsType) j);
899  valid = true;
900  break;
901  }
902  }
903 
904  ASSERTL0(valid, std::string("Unable to correctly parse the points type: ").append(pointsStrings[i]).c_str());
905  }
906  }
907 
908  // Get numModes
909  std::vector<unsigned int> numModes;
910  bool UniOrder = false;
911 
912  if(strstr(numModesString.c_str(),"UNIORDER:"))
913  {
914  UniOrder = true;
915  }
916 
917  valid = ParseUtils::GenerateOrderedVector(numModesString.c_str()+9, numModes);
918  ASSERTL0(valid, "Unable to correctly parse the number of modes.");
919 
920  // Get numPoints
921  std::vector<unsigned int> numPoints;
922  if(numPointDef)
923  {
924  valid = ParseUtils::GenerateOrderedVector(numPointsString.c_str(), numPoints);
925  ASSERTL0(valid, "Unable to correctly parse the number of points.");
926  }
927 
928  // Get fields names
929  std::vector<std::string> Fields;
930  valid = ParseUtils::GenerateOrderedStringVector(fieldsString.c_str(), Fields);
931  ASSERTL0(valid, "Unable to correctly parse the number of fields.");
932 
933  FieldDefinitionsSharedPtr fielddef = MemoryManager<FieldDefinitions>::AllocateSharedPtr(shape, elementIds, basis, UniOrder, numModes, Fields, numHomoDir, homoLengths, homoZIDs, homoYIDs, points, pointDef, numPoints, numPointDef);
934 
935  fielddefs.push_back(fielddef);
936 
937  element = element->NextSiblingElement("ELEMENTS");
938  }
939  loopXml = loopXml->NextSiblingElement(strLoop);
940  }
941  }
942 
943 
944  /**
945  *
946  */
947  void FieldIO::ImportFieldData(TiXmlDocument &doc, const std::vector<FieldDefinitionsSharedPtr> &fielddefs, std::vector<std::vector<NekDouble> > &fielddata)
948  {
949  int cntdumps = 0;
950 
951  TiXmlHandle docHandle(&doc);
952  TiXmlElement* master = NULL; // Master tag within which all data is contained.
953 
954  master = doc.FirstChildElement("NEKTAR");
955  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
956 
957  // Loop through all nektar tags, finding all of the element tags.
958  while (master)
959  {
960  TiXmlElement* element = master->FirstChildElement("ELEMENTS");
961  ASSERTL0(element, "Unable to find ELEMENTS tag within nektar tag.");
962  while (element)
963  {
964  // Extract the body, which the "data".
965  TiXmlNode* elementChild = element->FirstChild();
966  ASSERTL0(elementChild, "Unable to extract the data from the element tag.");
967  std::string elementStr;
968  while(elementChild)
969  {
970  if (elementChild->Type() == TiXmlNode::TEXT)
971  {
972  elementStr += elementChild->ToText()->ValueStr();
973  }
974  elementChild = elementChild->NextSibling();
975  }
976 
977  // Convert from base64 to binary.
978  typedef boost::archive::iterators::transform_width<
979  boost::archive::iterators::binary_from_base64<
980  std::string::const_iterator>, 8, 6 > binary_t;
981  std::string vCompressed(binary_t(elementStr.begin()),
982  binary_t(elementStr.end()));
983 
984  std::vector<NekDouble> elementFieldData;
985  ASSERTL0(Z_OK == Inflate(vCompressed, elementFieldData),
986  "Failed to decompress field data.");
987 
988  fielddata.push_back(elementFieldData);
989 
990  int datasize = CheckFieldDefinition(fielddefs[cntdumps]);
991  ASSERTL0(fielddata[cntdumps].size() == datasize*fielddefs[cntdumps]->m_fields.size(),"Input data is not the same length as header infoarmation");
992 
993  cntdumps++;
994 
995  element = element->NextSiblingElement("ELEMENTS");
996  }
997  master = master->NextSiblingElement("NEKTAR");
998  }
999  }
1000 
1001 
1002  /**
1003  * \brief add information about provenance and fieldmetadata
1004  */
1005  void FieldIO::AddInfoTag(TiXmlElement * root,
1006  const FieldMetaDataMap &fieldmetadatamap)
1007  {
1008  FieldMetaDataMap ProvenanceMap;
1009 
1010  // Nektar++ release version from VERSION file
1011  ProvenanceMap["NektarVersion"] = string(NEKTAR_VERSION);
1012 
1013  // Date/time stamp
1014  ptime::time_facet *facet = new ptime::time_facet("%d-%b-%Y %H:%M:%S");
1015  std::stringstream wss;
1016  wss.imbue(locale(wss.getloc(), facet));
1017  wss << ptime::second_clock::local_time();
1018  ProvenanceMap["Timestamp"] = wss.str();
1019 
1020  // Hostname
1021  boost::system::error_code ec;
1022  ProvenanceMap["Hostname"] = ip::host_name(ec);
1023 
1024  // Git information
1025  // If built from a distributed package, do not include this
1026  if (NekConstants::kGitSha1 != "GITDIR-NOTFOUND")
1027  {
1028  ProvenanceMap["GitSHA1"] = NekConstants::kGitSha1;
1029  ProvenanceMap["GitBranch"] = NekConstants::kGitBranch;
1030  }
1031 
1032  TiXmlElement * infoTag = new TiXmlElement("Metadata");
1033  root->LinkEndChild(infoTag);
1034 
1035  TiXmlElement * v;
1036  FieldMetaDataMap::const_iterator infoit;
1037 
1038  TiXmlElement * provTag = new TiXmlElement("Provenance");
1039  infoTag->LinkEndChild(provTag);
1040  for (infoit = ProvenanceMap.begin(); infoit != ProvenanceMap.end(); ++infoit)
1041  {
1042  v = new TiXmlElement( (infoit->first).c_str() );
1043  v->LinkEndChild(new TiXmlText((infoit->second).c_str()));
1044  provTag->LinkEndChild(v);
1045  }
1046 
1047  //---------------------------------------------
1048  // write field info section
1049  if(fieldmetadatamap != NullFieldMetaDataMap)
1050  {
1051  for(infoit = fieldmetadatamap.begin(); infoit != fieldmetadatamap.end(); ++infoit)
1052  {
1053  v = new TiXmlElement( (infoit->first).c_str() );
1054  v->LinkEndChild(new TiXmlText((infoit->second).c_str()));
1055  infoTag->LinkEndChild(v);
1056  }
1057  }
1058  }
1059 
1060 
1061  /**
1062  *
1063  */
1064  void FieldIO::GenerateSeqString(const std::vector<unsigned int> &elmtids,
1065  std::string &idString)
1066  {
1067  std::stringstream idStringStream;
1068  bool setdash = true;
1069  unsigned int endval;
1070 
1071  idStringStream << elmtids[0];
1072  for (int i = 1; i < elmtids.size(); ++i)
1073  {
1074  if(elmtids[i] == elmtids[i-1]+1)
1075  {
1076  if(setdash)
1077  {
1078  idStringStream << "-";
1079  setdash = false;
1080  }
1081 
1082  if(i == elmtids.size()-1) // last element
1083  {
1084  idStringStream << elmtids[i];
1085  }
1086  else
1087  {
1088  endval = elmtids[i];
1089  }
1090  }
1091  else
1092  {
1093  if(setdash == false) // finish off previous dash sequence
1094  {
1095  idStringStream << endval;
1096  setdash = true;
1097  }
1098 
1099 
1100  idStringStream << "," << elmtids[i];
1101  }
1102  }
1103  idString = idStringStream.str();
1104  }
1105 
1106 
1107  /**
1108  *
1109  */
1110  std::string FieldIO::SetUpOutput(const std::string outname,
1111  const std::vector<FieldDefinitionsSharedPtr>& fielddefs,
1112  const FieldMetaDataMap &fieldmetadatamap)
1113  {
1114  ASSERTL0(!outname.empty(), "Empty path given to SetUpOutput()");
1115 
1116  int nprocs = m_comm->GetSize();
1117  int rank = m_comm->GetRank();
1118 
1119  // Directory name if in parallel, regular filename if in serial
1120  fs::path specPath (outname);
1121 
1122  // Remove any existing file which is in the way
1123  if(m_comm->RemoveExistingFiles())
1124  {
1125  try
1126  {
1127  fs::remove_all(specPath);
1128  }
1129  catch (fs::filesystem_error& e)
1130  {
1131  ASSERTL0(e.code().value() == berrc::no_such_file_or_directory,
1132  "Filesystem error: " + string(e.what()));
1133  }
1134  }
1135 
1136  // serial processing just add ending.
1137  if(nprocs == 1)
1138  {
1139  cout << "Writing: " << specPath << endl;
1140  return LibUtilities::PortablePath(specPath);
1141  }
1142 
1143  // Compute number of elements on this process and share with other
1144  // processes. Also construct list of elements on this process from
1145  // available vector of field definitions.
1146  std::vector<unsigned int> elmtnums(nprocs,0);
1147  std::vector<unsigned int> idlist;
1148  int i;
1149  for (i = 0; i < fielddefs.size(); ++i)
1150  {
1151  elmtnums[rank] += fielddefs[i]->m_elementIDs.size();
1152  idlist.insert(idlist.end(), fielddefs[i]->m_elementIDs.begin(),
1153  fielddefs[i]->m_elementIDs.end());
1154  }
1155  m_comm->AllReduce(elmtnums,LibUtilities::ReduceMax);
1156 
1157  // Create the destination directory
1158  try
1159  {
1160  fs::create_directory(specPath);
1161  }
1162  catch (fs::filesystem_error& e)
1163  {
1164  ASSERTL0(false, "Filesystem error: " + string(e.what()));
1165  }
1166 
1167  // Collate per-process element lists on root process to generate
1168  // the info file.
1169  if (rank == 0)
1170  {
1171  std::vector<std::vector<unsigned int> > ElementIDs(nprocs);
1172 
1173  // Populate the list of element ID lists from all processes
1174  ElementIDs[0] = idlist;
1175  for (i = 1; i < nprocs; ++i)
1176  {
1177  std::vector<unsigned int> tmp(elmtnums[i]);
1178  m_comm->Recv(i, tmp);
1179  ElementIDs[i] = tmp;
1180  }
1181 
1182  // Set up output names
1183  std::vector<std::string> filenames;
1184  for(int i = 0; i < nprocs; ++i)
1185  {
1186  boost::format pad("P%1$07d.fld");
1187  pad % i;
1188  filenames.push_back(pad.str());
1189  }
1190 
1191  // Write the Info.xml file
1192  string infofile = LibUtilities::PortablePath(
1193  specPath / fs::path("Info.xml"));
1194 
1195  cout << "Writing: " << specPath << endl;
1196  WriteMultiFldFileIDs(infofile, filenames, ElementIDs,
1197  fieldmetadatamap);
1198  }
1199  else
1200  {
1201  // Send this process's ID list to the root process
1202  m_comm->Send(0, idlist);
1203  }
1204 
1205  // Pad rank to 8char filenames, e.g. P0000000.fld
1206  boost::format pad("P%1$07d.fld");
1207  pad % m_comm->GetRank();
1208 
1209  // Generate full path name
1210  fs::path poutfile(pad.str());
1211  fs::path fulloutname = specPath / poutfile;
1212 
1213  // Return the full path to the partition for this process
1214  return LibUtilities::PortablePath(fulloutname);
1215  }
1216 
1217 
1218  /**
1219  * Compress a vector of NekDouble values into a string using zlib.
1220  */
1221  int FieldIO::Deflate(std::vector<NekDouble>& in,
1222  string& out)
1223  {
1224  int ret;
1225  unsigned have;
1226  z_stream strm;
1227  unsigned char* input = (unsigned char*)(&in[0]);
1228  string buffer;
1229  buffer.resize(CHUNK);
1230 
1231  /* allocate deflate state */
1232  strm.zalloc = Z_NULL;
1233  strm.zfree = Z_NULL;
1234  strm.opaque = Z_NULL;
1235  ret = deflateInit(&strm, Z_DEFAULT_COMPRESSION);
1236 
1237  ASSERTL0(ret == Z_OK, "Error initializing Zlib.");
1238 
1239  strm.avail_in = in.size() * sizeof(NekDouble) / sizeof(char);
1240  strm.next_in = input;
1241 
1242  // Deflate input until output buffer is no longer full.
1243  do {
1244  strm.avail_out = CHUNK;
1245  strm.next_out = (unsigned char*)(&buffer[0]);
1246 
1247  ret = deflate(&strm, Z_FINISH);
1248 
1249  // Deflate can return Z_OK, Z_STREAM_ERROR, Z_BUF_ERROR or
1250  // Z_STREAM_END. All, except Z_STREAM_ERROR are ok.
1251  ASSERTL0(ret != Z_STREAM_ERROR, "Zlib stream error");
1252 
1253  have = CHUNK - strm.avail_out;
1254  out += buffer.substr(0, have);
1255 
1256  } while (strm.avail_out == 0);
1257 
1258  // Check all input was processed.
1259  ASSERTL0(strm.avail_in == 0, "Not all input was used.");
1260 
1261  // Check stream is complete.
1262  ASSERTL0(ret == Z_STREAM_END, "Stream not finished");
1263 
1264  // Clean-up and return
1265  (void)deflateEnd(&strm);
1266  return Z_OK;
1267  }
1268 
1269 
1270  /**
1271  * Decompress a zlib-compressed string into a vector of NekDouble
1272  * values.
1273  */
1274  int FieldIO::Inflate(std::string& in,
1275  std::vector<NekDouble>& out)
1276  {
1277  int ret;
1278  unsigned have;
1279  z_stream strm;
1280  string buffer;
1281  buffer.resize(CHUNK);
1282  string output;
1283 
1284  strm.zalloc = Z_NULL;
1285  strm.zfree = Z_NULL;
1286  strm.opaque = Z_NULL;
1287  strm.avail_in = 0;
1288  strm.next_in = Z_NULL;
1289  ret = inflateInit(&strm);
1290  ASSERTL0(ret == Z_OK, "Error initializing zlib decompression.");
1291 
1292  strm.avail_in = in.size();
1293  strm.next_in = (unsigned char*)(&in[0]);
1294 
1295  do {
1296  strm.avail_out = CHUNK;
1297  strm.next_out = (unsigned char*)(&buffer[0]);
1298 
1299  ret = inflate(&strm, Z_NO_FLUSH);
1300 
1301  ASSERTL0(ret != Z_STREAM_ERROR, "Stream error occured.");
1302 
1303  switch (ret) {
1304  case Z_NEED_DICT:
1305  ret = Z_DATA_ERROR; /* and fall through */
1306  case Z_DATA_ERROR:
1307  case Z_MEM_ERROR:
1308  (void)inflateEnd(&strm);
1309  return ret;
1310  }
1311 
1312  have = CHUNK - strm.avail_out;
1313  output += buffer.substr(0, have);
1314 
1315  } while (strm.avail_out == 0);
1316 
1317  (void)inflateEnd(&strm);
1318 
1319  if (ret == Z_STREAM_END)
1320  {
1321  NekDouble* readFieldData = (NekDouble*) output.c_str();
1322  unsigned int len = output.size() * sizeof(*output.c_str())
1323  / sizeof(NekDouble);
1324  out.assign( readFieldData, readFieldData + len);
1325  return Z_OK;
1326  }
1327  else
1328  {
1329  return Z_DATA_ERROR;
1330  }
1331  }
1332 
1333  /**
1334  *
1335  */
1336  int FieldIO::CheckFieldDefinition(const FieldDefinitionsSharedPtr &fielddefs)
1337  {
1338  int i;
1339 
1340  if(fielddefs->m_elementIDs.size() == 0) // empty partition
1341  {
1342  return 0;
1343  }
1344  //ASSERTL0(fielddefs->m_elementIDs.size() > 0, "Fielddefs vector must contain at least one element of data .");
1345 
1346  unsigned int numbasis = 0;
1347 
1348  // Determine nummodes vector lists are correct length
1349  switch(fielddefs->m_shapeType)
1350  {
1351  case eSegment:
1352  numbasis = 1;
1353  if(fielddefs->m_numHomogeneousDir)
1354  {
1355  numbasis += fielddefs->m_numHomogeneousDir;
1356  }
1357 
1358  break;
1359  case eTriangle:
1360  case eQuadrilateral:
1361  if(fielddefs->m_numHomogeneousDir)
1362  {
1363  numbasis = 3;
1364  }
1365  else
1366  {
1367  numbasis = 2;
1368  }
1369  break;
1370  case eTetrahedron:
1371  case ePyramid:
1372  case ePrism:
1373  case eHexahedron:
1374  numbasis = 3;
1375  break;
1376  default:
1377  ASSERTL0(false, "Unsupported shape type.");
1378  break;
1379  }
1380 
1381  unsigned int datasize = 0;
1382 
1383  ASSERTL0(fielddefs->m_basis.size() == numbasis, "Length of basis vector is incorrect");
1384 
1385  if(fielddefs->m_uniOrder == true)
1386  {
1387  unsigned int cnt = 0;
1388  // calculate datasize
1389  switch(fielddefs->m_shapeType)
1390  {
1391  case eSegment:
1392  {
1393  int l = fielddefs->m_numModes[cnt++];
1394  if(fielddefs->m_numHomogeneousDir == 1)
1395  {
1396  datasize += l*fielddefs->m_numModes[cnt++];
1397  }
1398  else if(fielddefs->m_numHomogeneousDir == 2)
1399  {
1400  int m = fielddefs->m_numModes[cnt++];
1401  datasize += l*m*fielddefs->m_numModes[cnt++];
1402  }
1403  else
1404  {
1405  datasize += l;
1406  }
1407  }
1408  break;
1409  case eTriangle:
1410  {
1411  int l = fielddefs->m_numModes[cnt++];
1412  int m = fielddefs->m_numModes[cnt++];
1413 
1414  if(fielddefs->m_numHomogeneousDir == 1)
1415  {
1416  datasize += StdTriData::getNumberOfCoefficients(l,m)*
1417  fielddefs->m_homogeneousZIDs.size();
1418  }
1419  else
1420  {
1421  datasize += StdTriData::getNumberOfCoefficients(l,m);
1422  }
1423  }
1424  break;
1425  case eQuadrilateral:
1426  {
1427  int l = fielddefs->m_numModes[cnt++];
1428  int m = fielddefs->m_numModes[cnt++];
1429  if(fielddefs->m_numHomogeneousDir == 1)
1430  {
1431  datasize += l*m*fielddefs->m_homogeneousZIDs.size();
1432  }
1433  else
1434  {
1435  datasize += l*m;
1436  }
1437  }
1438  break;
1439  case eTetrahedron:
1440  {
1441  int l = fielddefs->m_numModes[cnt++];
1442  int m = fielddefs->m_numModes[cnt++];
1443  int n = fielddefs->m_numModes[cnt++];
1444  datasize += StdTetData::getNumberOfCoefficients(l,m,n);
1445  }
1446  break;
1447  case ePyramid:
1448  {
1449  int l = fielddefs->m_numModes[cnt++];
1450  int m = fielddefs->m_numModes[cnt++];
1451  int n = fielddefs->m_numModes[cnt++];
1452  datasize += StdPyrData::getNumberOfCoefficients(l,m,n);
1453  }
1454  break;
1455  case ePrism:
1456  {
1457  int l = fielddefs->m_numModes[cnt++];
1458  int m = fielddefs->m_numModes[cnt++];
1459  int n = fielddefs->m_numModes[cnt++];
1460  datasize += StdPrismData::getNumberOfCoefficients(l,m,n);
1461  }
1462  break;
1463  case eHexahedron:
1464  {
1465  int l = fielddefs->m_numModes[cnt++];
1466  int m = fielddefs->m_numModes[cnt++];
1467  int n = fielddefs->m_numModes[cnt++];
1468  datasize += l*m*n;
1469  }
1470  break;
1471  default:
1472  ASSERTL0(false, "Unsupported shape type.");
1473  break;
1474  }
1475 
1476  datasize *= fielddefs->m_elementIDs.size();
1477  }
1478  else
1479  {
1480  unsigned int cnt = 0;
1481  // calculate data length
1482  for(i = 0; i < fielddefs->m_elementIDs.size(); ++i)
1483  {
1484  switch(fielddefs->m_shapeType)
1485  {
1486  case eSegment:
1487  datasize += fielddefs->m_numModes[cnt++];
1488  break;
1489  case eTriangle:
1490  {
1491  int l = fielddefs->m_numModes[cnt++];
1492  int m = fielddefs->m_numModes[cnt++];
1493  datasize += StdTriData::getNumberOfCoefficients(l,m);
1494  }
1495  break;
1496  case eQuadrilateral:
1497  {
1498  int l = fielddefs->m_numModes[cnt++];
1499  int m = fielddefs->m_numModes[cnt++];
1500  datasize += l*m;
1501  }
1502  break;
1503  case eTetrahedron:
1504  {
1505  int l = fielddefs->m_numModes[cnt++];
1506  int m = fielddefs->m_numModes[cnt++];
1507  int n = fielddefs->m_numModes[cnt++];
1508  datasize += StdTetData::getNumberOfCoefficients(l,m,n);
1509  }
1510  break;
1511  case ePyramid:
1512  {
1513  int l = fielddefs->m_numModes[cnt++];
1514  int m = fielddefs->m_numModes[cnt++];
1515  int n = fielddefs->m_numModes[cnt++];
1516  datasize += StdPyrData::getNumberOfCoefficients(l,m,n);
1517  }
1518  break;
1519  case ePrism:
1520  {
1521  int l = fielddefs->m_numModes[cnt++];
1522  int m = fielddefs->m_numModes[cnt++];
1523  int n = fielddefs->m_numModes[cnt++];
1524  datasize += StdPrismData::getNumberOfCoefficients(l,m,n);
1525  }
1526  break;
1527  case eHexahedron:
1528  {
1529  int l = fielddefs->m_numModes[cnt++];
1530  int m = fielddefs->m_numModes[cnt++];
1531  int n = fielddefs->m_numModes[cnt++];
1532  datasize += l*m*n;
1533  }
1534  break;
1535  default:
1536  ASSERTL0(false, "Unsupported shape type.");
1537  break;
1538  }
1539  }
1540  }
1541 
1542  return datasize;
1543  }
1544  }
1545 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2