Nektar++
MeshGraphHDF5.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MeshGraphHDF5.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: HDF5-based mesh format for Nektar++.
32 //
33 ////////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/algorithm/string.hpp>
36 #include <boost/core/ignore_unused.hpp>
37 #include <boost/filesystem.hpp>
38 #include <tinyxml.h>
39 #include <type_traits>
40 
45 
46 #define TIME_RESULT(verb, msg, timer) \
47  if (verb) \
48  { \
49  std::cout << " - " << msg << ": " << timer.TimePerTest(1) \
50  << std::endl; \
51  }
52 
53 using namespace std;
54 using namespace Nektar::LibUtilities;
55 
56 namespace Nektar
57 {
58 namespace SpatialDomains
59 {
60 
61 /// Version of the Nektar++ HDF5 geometry format, which is embedded into the
62 /// main NEKTAR/GEOMETRY group as an attribute.
63 const unsigned int MeshGraphHDF5::FORMAT_VERSION = 1;
64 
65 std::string MeshGraphHDF5::className =
66  GetMeshGraphFactory().RegisterCreatorFunction("HDF5", MeshGraphHDF5::create,
67  "IO with HDF5 geometry");
68 
69 void MeshGraphHDF5::ReadGeometry(DomainRangeShPtr rng, bool fillGraph)
70 {
71  boost::ignore_unused(rng);
72 
73  ReadComposites();
74  ReadDomain();
75 
76  // Close up shop.
77  m_mesh->Close();
78  m_maps->Close();
79  m_file->Close();
80 
81  if (fillGraph)
82  {
83  MeshGraph::FillGraph();
84  }
85 }
86 
87 /**
88  * @brief Utility function to split a vector equally amongst a number of
89  * processors.
90  *
91  * @param vecsize Size of the total amount of work
92  * @param rank Rank of this process
93  * @param nprocs Number of processors in the group
94  *
95  * @return A pair with the offset this process should occupy, along with the
96  * count for the amount of work.
97  */
98 std::pair<size_t, size_t> SplitWork(size_t vecsize, int rank, int nprocs)
99 {
100  size_t div = vecsize / nprocs;
101  size_t rem = vecsize % nprocs;
102  if (rank < rem)
103  {
104  return std::make_pair(rank * (div + 1), div + 1);
105  }
106  else
107  {
108  return std::make_pair((rank - rem) * div + rem * (div + 1), div);
109  }
110 }
111 
112 template <class T, typename std::enable_if<T::kDim == 0, int>::type = 0>
113 inline int GetGeomDataDim(std::map<int, std::shared_ptr<T>> &geomMap)
114 {
115  boost::ignore_unused(geomMap);
116  return 3;
117 }
118 
119 template <class T, typename std::enable_if<T::kDim == 1, int>::type = 0>
120 inline int GetGeomDataDim(std::map<int, std::shared_ptr<T>> &geomMap)
121 {
122  boost::ignore_unused(geomMap);
123  return T::kNverts;
124 }
125 
126 template <class T, typename std::enable_if<T::kDim == 2, int>::type = 0>
127 inline int GetGeomDataDim(std::map<int, std::shared_ptr<T>> &geomMap)
128 {
129  boost::ignore_unused(geomMap);
130  return T::kNedges;
131 }
132 
133 template <class T, typename std::enable_if<T::kDim == 3, int>::type = 0>
134 inline int GetGeomDataDim(std::map<int, std::shared_ptr<T>> &geomMap)
135 {
136  boost::ignore_unused(geomMap);
137  return T::kNfaces;
138 }
139 
140 template <class... T> inline void UniqueValues(std::unordered_set<int> &unique)
141 {
142  boost::ignore_unused(unique);
143 }
144 
145 template <class... T>
146 inline void UniqueValues(std::unordered_set<int> &unique,
147  const std::vector<int> &input, T &...args)
148 {
149  for (auto i : input)
150  {
151  unique.insert(i);
152  }
153 
154  UniqueValues(unique, args...);
155 }
156 
157 std::string MeshGraphHDF5::cmdSwitch =
158  LibUtilities::SessionReader::RegisterCmdLineFlag(
159  "use-hdf5-node-comm", "",
160  "Use a per-node communicator for HDF5 partitioning.");
161 
162 /**
163  * @brief Partition the mesh
164  */
165 void MeshGraphHDF5::PartitionMesh(LibUtilities::SessionReaderSharedPtr session)
166 {
168  all.Start();
169  int err;
170  LibUtilities::CommSharedPtr comm = session->GetComm();
171  int rank = comm->GetRank(), nproc = comm->GetSize();
172 
173  // By default, only the root process will have read the session file, which
174  // is done to avoid every process needing to read the XML file. For HDF5, we
175  // don't care about this, so just have every process parse the session file.
176  if (rank > 0)
177  {
178  session->InitSession();
179  }
180 
181  // We use the XML geometry to find information about the HDF5 file.
182  m_session = session;
183  m_xmlGeom = m_session->GetElement("NEKTAR/GEOMETRY");
184  TiXmlAttribute *attr = m_xmlGeom->FirstAttribute();
185  m_meshPartitioned = true;
186  m_meshDimension = 3;
187  m_spaceDimension = 3;
188 
189  while (attr)
190  {
191  std::string attrName(attr->Name());
192  if (attrName == "DIM")
193  {
194  err = attr->QueryIntValue(&m_meshDimension);
195  ASSERTL0(err == TIXML_SUCCESS, "Unable to read mesh dimension.");
196  }
197  else if (attrName == "SPACE")
198  {
199  err = attr->QueryIntValue(&m_spaceDimension);
200  ASSERTL0(err == TIXML_SUCCESS, "Unable to read space dimension.");
201  }
202  else if (attrName == "PARTITION")
203  {
204  ASSERTL0(false,
205  "PARTITION parameter should only be used in XML meshes");
206  }
207  else if (attrName == "HDF5FILE")
208  {
209  m_hdf5Name = attr->Value();
210  ASSERTL1(err == TIXML_SUCCESS, "Unable to read hdf5 name.");
211  }
212  else if (attrName == "PARTITIONED")
213  {
214  ASSERTL0(false,
215  "PARTITIONED parameter should only be used in XML meshes");
216  }
217  else
218  {
219  std::string errstr("Unknown attribute: ");
220  errstr += attrName;
221  ASSERTL1(false, errstr.c_str());
222  }
223  // Get the next attribute.
224  attr = attr->Next();
225  }
226 
227  ASSERTL0(m_hdf5Name.size() > 0, "Unable to obtain mesh file name.");
228  ASSERTL0(m_meshDimension <= m_spaceDimension,
229  "Mesh dimension greater than space dimension.");
230 
231  // Open handle to the HDF5 mesh
232  LibUtilities::H5::PListSharedPtr parallelProps = H5::PList::Default();
233  m_readPL = H5::PList::Default();
234 
235  if (nproc > 1)
236  {
237  // Use MPI/O to access the file
238  parallelProps = H5::PList::FileAccess();
239  parallelProps->SetMpio(comm);
240  // Use collective IO
241  m_readPL = H5::PList::DatasetXfer();
242  m_readPL->SetDxMpioCollective();
243  }
244 
245  m_file = H5::File::Open(m_hdf5Name, H5F_ACC_RDONLY, parallelProps);
246 
247  auto root = m_file->OpenGroup("NEKTAR");
248  ASSERTL0(root, "Cannot find NEKTAR group in HDF5 file.");
249 
250  auto root2 = root->OpenGroup("GEOMETRY");
251  ASSERTL0(root2, "Cannot find NEKTAR/GEOMETRY group in HDF5 file.");
252 
253  // Check format version
254  unsigned int formatVersion;
255  H5::Group::AttrIterator attrIt = root2->attr_begin();
256  H5::Group::AttrIterator attrEnd = root2->attr_end();
257  for (; attrIt != attrEnd; ++attrIt)
258  {
259  if (*attrIt == "FORMAT_VERSION")
260  {
261  break;
262  }
263  }
264  ASSERTL0(attrIt != attrEnd,
265  "Unable to determine Nektar++ geometry HDF5 file version.");
266  root2->GetAttribute("FORMAT_VERSION", formatVersion);
267 
268  ASSERTL0(formatVersion <= FORMAT_VERSION,
269  "File format in " + m_hdf5Name +
270  " is higher than supported in "
271  "this version of Nektar++");
272 
273  m_mesh = root2->OpenGroup("MESH");
274  ASSERTL0(m_mesh, "Cannot find NEKTAR/GEOMETRY/MESH group in HDF5 file.");
275  m_maps = root2->OpenGroup("MAPS");
276  ASSERTL0(m_mesh, "Cannot find NEKTAR/GEOMETRY/MAPS group in HDF5 file.");
277 
278  // Depending on dimension, read element IDs.
279  std::map<int,
280  std::vector<std::tuple<std::string, int, LibUtilities::ShapeType>>>
281  dataSets;
282 
283  dataSets[1] = {make_tuple("SEG", 2, LibUtilities::eSegment)};
284  dataSets[2] = {make_tuple("TRI", 3, LibUtilities::eTriangle),
285  make_tuple("QUAD", 4, LibUtilities::eQuadrilateral)};
286  dataSets[3] = {make_tuple("TET", 4, LibUtilities::eTetrahedron),
287  make_tuple("PYR", 5, LibUtilities::ePyramid),
288  make_tuple("PRISM", 5, LibUtilities::ePrism),
289  make_tuple("HEX", 6, LibUtilities::eHexahedron)};
290 
291  bool verbRoot = rank == 0 && m_session->DefinesCmdLineArgument("verbose");
292 
293  if (verbRoot)
294  {
295  std::cout << "Reading HDF5 geometry..." << std::endl;
296  }
297 
298  // If we want to use an multi-level communicator, then split the
299  // communicator at this point. We set by default the inter-node communicator
300  // to be the normal communicator: this way if the multi-level partitioning
301  // is disabled we proceed as 'normal'.
302  LibUtilities::CommSharedPtr innerComm, interComm = comm;
303  int innerRank = 0, innerSize = 1, interRank = rank, interSize = nproc;
304 
305  if (session->DefinesCmdLineArgument("use-hdf5-node-comm"))
306  {
307  auto splitComm = comm->SplitCommNode();
308  innerComm = splitComm.first;
309  interComm = splitComm.second;
310  innerRank = innerComm->GetRank();
311  innerSize = innerComm->GetSize();
312 
313  if (innerRank == 0)
314  {
315  interRank = interComm->GetRank();
316  interSize = interComm->GetSize();
317  }
318  }
319 
320  // Unordered set of rows of the dataset this process needs to read for the
321  // elements of dimension m_meshDimension.
322  std::unordered_set<int> toRead;
323 
324  // Calculate reasonably even distribution of processors for calling
325  // ptScotch. We'll do this work on only one process per node.
327  t.Start();
328 
329  if (verbRoot)
330  {
331  std::cout << " - beginning partitioning" << std::endl;
332  }
333 
334  // Perform initial read if either (a) we are on all ranks and multi-level
335  // partitioning is not enabled; or (b) rank 0 of all nodes if it is.
336  if (innerRank == 0)
337  {
339  t2.Start();
340 
341  bool verbRoot2 =
342  session->DefinesCmdLineArgument("verbose") && interRank == 0;
343 
344  // Read IDs for partitioning purposes
345  std::vector<int> ids;
346 
347  // Map from element ID to 'row' which is a contiguous ordering required
348  // for parallel partitioning.
349  std::vector<MeshEntity> elmts;
350  std::unordered_map<int, int> row2id, id2row;
351 
352  LibUtilities::H5::FileSharedPtr file = m_file;
353  LibUtilities::H5::PListSharedPtr readPL = m_readPL;
354  LibUtilities::H5::GroupSharedPtr mesh = m_mesh, maps = m_maps;
355 
356  if (innerComm)
357  {
358  // For per-node partitioning, create a temporary reader (otherwise
359  // communicators are inconsistent).
360  auto parallelProps = H5::PList::FileAccess();
361  parallelProps->SetMpio(interComm);
362 
363  // Use collective IO
364  readPL = H5::PList::DatasetXfer();
365  readPL->SetDxMpioCollective();
366  file = H5::File::Open(m_hdf5Name, H5F_ACC_RDONLY, parallelProps);
367 
368  auto root = file->OpenGroup("NEKTAR");
369  auto root2 = root->OpenGroup("GEOMETRY");
370  mesh = root2->OpenGroup("MESH");
371  maps = root2->OpenGroup("MAPS");
372  }
373 
374  int rowCount = 0;
375  for (auto &it : dataSets[m_meshDimension])
376  {
377  std::string ds = std::get<0>(it);
378 
379  if (!mesh->ContainsDataSet(ds))
380  {
381  continue;
382  }
383 
384  // Open metadata dataset
385  H5::DataSetSharedPtr data = mesh->OpenDataSet(ds);
386  H5::DataSpaceSharedPtr space = data->GetSpace();
387  vector<hsize_t> dims = space->GetDims();
388 
389  H5::DataSetSharedPtr mdata = maps->OpenDataSet(ds);
390  H5::DataSpaceSharedPtr mspace = mdata->GetSpace();
391  vector<hsize_t> mdims = mspace->GetDims();
392 
393  // TODO: This could perhaps be done more intelligently; reads all
394  // IDs for the top-level elements so that we can construct the dual
395  // graph of the mesh.
396  vector<int> tmpElmts, tmpIds;
397  mdata->Read(tmpIds, mspace, readPL);
398  data->Read(tmpElmts, space, readPL);
399 
400  const int nGeomData = std::get<1>(it);
401 
402  for (int i = 0, cnt = 0; i < tmpIds.size(); ++i, ++rowCount)
403  {
404  MeshEntity e;
405  row2id[rowCount] = tmpIds[i];
406  id2row[tmpIds[i]] = row2id[rowCount];
407  e.id = rowCount;
408  e.origId = tmpIds[i];
409  e.ghost = false;
410  e.list = std::vector<unsigned int>(&tmpElmts[cnt],
411  &tmpElmts[cnt + nGeomData]);
412  elmts.push_back(e);
413  cnt += nGeomData;
414  }
415  }
416 
417  interComm->Block();
418 
419  t2.Stop();
420  TIME_RESULT(verbRoot2, " - initial read", t2);
421  t2.Start();
422 
423  // Check to see we have at least as many processors as elements.
424  size_t numElmt = elmts.size();
425  ASSERTL0(nproc <= numElmt,
426  "This mesh has more processors than elements!");
427 
428  auto elRange = SplitWork(numElmt, interRank, interSize);
429 
430  // Construct map of element entities for partitioner.
431  std::map<int, MeshEntity> partElmts;
432  std::unordered_set<int> facetIDs;
433 
434  int vcnt = 0;
435 
436  for (int el = elRange.first; el < elRange.first + elRange.second;
437  ++el, ++vcnt)
438  {
439  MeshEntity elmt = elmts[el];
440  elmt.ghost = false;
441  partElmts[el] = elmt;
442 
443  for (auto &facet : elmt.list)
444  {
445  facetIDs.insert(facet);
446  }
447  }
448 
449  // Now identify ghost vertices for the graph. This could also probably
450  // be improved.
451  int nLocal = vcnt;
452  for (int i = 0; i < numElmt; ++i)
453  {
454  // Ignore anything we already read.
455  if (i >= elRange.first && i < elRange.first + elRange.second)
456  {
457  continue;
458  }
459 
460  MeshEntity elmt = elmts[i];
461  bool insert = false;
462 
463  // Check for connections to local elements.
464  for (auto &eId : elmt.list)
465  {
466  if (facetIDs.find(eId) != facetIDs.end())
467  {
468  insert = true;
469  break;
470  }
471  }
472 
473  if (insert)
474  {
475  elmt.ghost = true;
476  partElmts[elmt.id] = elmt;
477  }
478  }
479 
480  // Create partitioner. Default partitioner to use is PtScotch. Use
481  // ParMetis as default if it is installed. Override default with
482  // command-line flags if they are set.
483  string partitionerName = nproc > 1 ? "PtScotch" : "Scotch";
484  if (GetMeshPartitionFactory().ModuleExists("ParMetis"))
485  {
486  partitionerName = "ParMetis";
487  }
488  if (session->DefinesCmdLineArgument("use-parmetis"))
489  {
490  partitionerName = "ParMetis";
491  }
492  if (session->DefinesCmdLineArgument("use-ptscotch"))
493  {
494  partitionerName = "PtScotch";
495  }
496 
497  MeshPartitionSharedPtr partitioner =
499  partitionerName, session, interComm, m_meshDimension, partElmts,
500  CreateCompositeDescriptor(id2row));
501 
502  t2.Stop();
503  TIME_RESULT(verbRoot2, " - partitioner setup", t2);
504  t2.Start();
505 
506  partitioner->PartitionMesh(interSize, true, false, nLocal);
507  t2.Stop();
508  TIME_RESULT(verbRoot2, " - partitioning", t2);
509  t2.Start();
510 
511  // Now construct a second graph that is partitioned in serial by this
512  // rank.
513  std::vector<unsigned int> nodeElmts;
514  partitioner->GetElementIDs(interRank, nodeElmts);
515 
516  if (innerSize > 1)
517  {
518  // Construct map of element entities for partitioner.
519  std::map<int, MeshEntity> partElmts;
520  std::unordered_map<int, int> row2elmtid, elmtid2row;
521 
522  int vcnt = 0;
523 
524  // We need to keep track of which elements in the new partition
525  // correspond to elemental IDs for later (in a similar manner to
526  // row2id).
527  for (auto &elmtRow : nodeElmts)
528  {
529  row2elmtid[vcnt] = elmts[elmtRow].origId;
530  elmtid2row[elmts[elmtRow].origId] = vcnt;
531  MeshEntity elmt = elmts[elmtRow];
532  elmt.ghost = false;
533  partElmts[vcnt++] = elmt;
534  }
535 
536  // Create temporary serial communicator for serial partitioning.
537  auto tmpComm =
539 
540  MeshPartitionSharedPtr partitioner =
542  "Scotch", session, tmpComm, m_meshDimension, partElmts,
543  CreateCompositeDescriptor(elmtid2row));
544 
545  t2.Stop();
546  TIME_RESULT(verbRoot2, " - inner partition setup", t2);
547  t2.Start();
548 
549  partitioner->PartitionMesh(innerSize, true, false, 0);
550 
551  t2.Stop();
552  TIME_RESULT(verbRoot2, " - inner partitioning", t2);
553  t2.Start();
554 
555  // Send contributions to remaining processors.
556  for (int i = 1; i < innerSize; ++i)
557  {
558  std::vector<unsigned int> tmp;
559  partitioner->GetElementIDs(i, tmp);
560  size_t tmpsize = tmp.size();
561  for (int j = 0; j < tmpsize; ++j)
562  {
563  tmp[j] = row2elmtid[tmp[j]];
564  }
565  innerComm->Send(i, tmpsize);
566  innerComm->Send(i, tmp);
567  }
568 
569  t2.Stop();
570  TIME_RESULT(verbRoot2, " - inner partition scatter", t2);
571 
572  std::vector<unsigned int> tmp;
573  partitioner->GetElementIDs(0, tmp);
574 
575  for (auto &tmpId : tmp)
576  {
577  toRead.insert(row2elmtid[tmpId]);
578  }
579  }
580  else
581  {
582  for (auto &tmpId : nodeElmts)
583  {
584  toRead.insert(row2id[tmpId]);
585  }
586  }
587  }
588  else
589  {
590  // For multi-level partitioning, the innermost rank receives its
591  // partitions from rank 0 on each node.
592  size_t tmpSize;
593  innerComm->Recv(0, tmpSize);
594  std::vector<unsigned int> tmp(tmpSize);
595  innerComm->Recv(0, tmp);
596 
597  for (auto &tmpId : tmp)
598  {
599  toRead.insert(tmpId);
600  }
601  }
602 
603  t.Stop();
604  TIME_RESULT(verbRoot, "partitioning total", t);
605 
606  // Since objects are going to be constructed starting from vertices, we now
607  // need to recurse down the geometry facet dimensions to figure out which
608  // rows to read from each dataset.
609  std::vector<int> vertIDs, segIDs, triIDs, quadIDs;
610  std::vector<int> tetIDs, prismIDs, pyrIDs, hexIDs;
611  std::vector<int> segData, triData, quadData, tetData;
612  std::vector<int> prismData, pyrData, hexData;
613  std::vector<NekDouble> vertData;
614 
615  if (m_meshDimension == 3)
616  {
617  t.Start();
618  // Read 3D data
619  ReadGeometryData(m_hexGeoms, "HEX", toRead, hexIDs, hexData);
620  ReadGeometryData(m_pyrGeoms, "PYR", toRead, pyrIDs, pyrData);
621  ReadGeometryData(m_prismGeoms, "PRISM", toRead, prismIDs, prismData);
622  ReadGeometryData(m_tetGeoms, "TET", toRead, tetIDs, tetData);
623 
624  toRead.clear();
625  UniqueValues(toRead, hexData, pyrData, prismData, tetData);
626  t.Stop();
627  TIME_RESULT(verbRoot, "read 3D elements", t);
628  }
629 
630  if (m_meshDimension >= 2)
631  {
632  t.Start();
633  // Read 2D data
634  ReadGeometryData(m_triGeoms, "TRI", toRead, triIDs, triData);
635  ReadGeometryData(m_quadGeoms, "QUAD", toRead, quadIDs, quadData);
636 
637  toRead.clear();
638  UniqueValues(toRead, triData, quadData);
639  t.Stop();
640  TIME_RESULT(verbRoot, "read 2D elements", t);
641  }
642 
643  if (m_meshDimension >= 1)
644  {
645  t.Start();
646  // Read 2D data
647  ReadGeometryData(m_segGeoms, "SEG", toRead, segIDs, segData);
648 
649  toRead.clear();
650  UniqueValues(toRead, segData);
651  t.Stop();
652  TIME_RESULT(verbRoot, "read 1D elements", t);
653  }
654 
655  t.Start();
656  ReadGeometryData(m_vertSet, "VERT", toRead, vertIDs, vertData);
657  t.Stop();
658  TIME_RESULT(verbRoot, "read 0D elements", t);
659 
660  // Now start to construct geometry objects, starting from vertices upwards.
661  t.Start();
662  FillGeomMap(m_vertSet, CurveMap(), vertIDs, vertData);
663  t.Stop();
664  TIME_RESULT(verbRoot, "construct 0D elements", t);
665 
666  if (m_meshDimension >= 1)
667  {
668  // Read curves
669  toRead.clear();
670  for (auto &edge : segIDs)
671  {
672  toRead.insert(edge);
673  }
674  ReadCurveMap(m_curvedEdges, "CURVE_EDGE", toRead);
675 
676  t.Start();
677  FillGeomMap(m_segGeoms, m_curvedEdges, segIDs, segData);
678  t.Stop();
679  TIME_RESULT(verbRoot, "construct 1D elements", t);
680  }
681 
682  if (m_meshDimension >= 2)
683  {
684  // Read curves
685  toRead.clear();
686  for (auto &face : triIDs)
687  {
688  toRead.insert(face);
689  }
690  for (auto &face : quadIDs)
691  {
692  toRead.insert(face);
693  }
694  ReadCurveMap(m_curvedFaces, "CURVE_FACE", toRead);
695 
696  t.Start();
697  FillGeomMap(m_triGeoms, m_curvedFaces, triIDs, triData);
698  FillGeomMap(m_quadGeoms, m_curvedFaces, quadIDs, quadData);
699  t.Stop();
700  TIME_RESULT(verbRoot, "construct 2D elements", t);
701  }
702 
703  if (m_meshDimension >= 3)
704  {
705  t.Start();
706  FillGeomMap(m_hexGeoms, CurveMap(), hexIDs, hexData);
707  FillGeomMap(m_prismGeoms, CurveMap(), prismIDs, prismData);
708  FillGeomMap(m_pyrGeoms, CurveMap(), pyrIDs, pyrData);
709  FillGeomMap(m_tetGeoms, CurveMap(), tetIDs, tetData);
710  t.Stop();
711  TIME_RESULT(verbRoot, "construct 3D elements", t);
712  }
713 
714  // Populate m_bndRegOrder.
715  if (m_session->DefinesElement("NEKTAR/CONDITIONS"))
716  {
717  std::set<int> vBndRegionIdList;
718  TiXmlElement *vConditions =
719  new TiXmlElement(*m_session->GetElement("Nektar/Conditions"));
720  TiXmlElement *vBndRegions =
721  vConditions->FirstChildElement("BOUNDARYREGIONS");
722  TiXmlElement *vItem;
723 
724  if (vBndRegions)
725  {
726  vItem = vBndRegions->FirstChildElement();
727  while (vItem)
728  {
729  std::string vSeqStr = vItem->FirstChild()->ToText()->Value();
730  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
731  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
732  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
733 
734  std::vector<unsigned int> vSeq;
735  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
736 
737  int p = atoi(vItem->Attribute("ID"));
738  m_bndRegOrder[p] = vSeq;
739  vItem = vItem->NextSiblingElement();
740  }
741  }
742  }
743 
744  all.Stop();
745  TIME_RESULT(verbRoot, "total time", all);
746 }
747 
748 template <class T, typename DataType>
749 void MeshGraphHDF5::ConstructGeomObject(
750  std::map<int, std::shared_ptr<T>> &geomMap, int id, DataType *data,
751  CurveSharedPtr curve)
752 {
753  boost::ignore_unused(geomMap, id, data, curve);
754 }
755 
756 template <>
757 void MeshGraphHDF5::ConstructGeomObject(
758  std::map<int, std::shared_ptr<PointGeom>> &geomMap, int id, NekDouble *data,
759  CurveSharedPtr curve)
760 {
761  boost::ignore_unused(curve);
763  m_spaceDimension, id, data[0], data[1], data[2]);
764 }
765 
766 template <>
767 void MeshGraphHDF5::ConstructGeomObject(
768  std::map<int, std::shared_ptr<SegGeom>> &geomMap, int id, int *data,
769  CurveSharedPtr curve)
770 {
771  PointGeomSharedPtr pts[2] = {GetVertex(data[0]), GetVertex(data[1])};
773  id, m_spaceDimension, pts, curve);
774 }
775 
776 template <>
777 void MeshGraphHDF5::ConstructGeomObject(
778  std::map<int, std::shared_ptr<TriGeom>> &geomMap, int id, int *data,
779  CurveSharedPtr curve)
780 {
781  SegGeomSharedPtr segs[3] = {GetSegGeom(data[0]), GetSegGeom(data[1]),
782  GetSegGeom(data[2])};
783  geomMap[id] = MemoryManager<TriGeom>::AllocateSharedPtr(id, segs, curve);
784 }
785 
786 template <>
787 void MeshGraphHDF5::ConstructGeomObject(
788  std::map<int, std::shared_ptr<QuadGeom>> &geomMap, int id, int *data,
789  CurveSharedPtr curve)
790 {
791  SegGeomSharedPtr segs[4] = {GetSegGeom(data[0]), GetSegGeom(data[1]),
792  GetSegGeom(data[2]), GetSegGeom(data[3])};
793  geomMap[id] = MemoryManager<QuadGeom>::AllocateSharedPtr(id, segs, curve);
794 }
795 
796 template <>
797 void MeshGraphHDF5::ConstructGeomObject(
798  std::map<int, std::shared_ptr<TetGeom>> &geomMap, int id, int *data,
799  CurveSharedPtr curve)
800 {
801  boost::ignore_unused(curve);
802  TriGeomSharedPtr faces[4] = {
803  std::static_pointer_cast<TriGeom>(GetGeometry2D(data[0])),
804  std::static_pointer_cast<TriGeom>(GetGeometry2D(data[1])),
805  std::static_pointer_cast<TriGeom>(GetGeometry2D(data[2])),
806  std::static_pointer_cast<TriGeom>(GetGeometry2D(data[3]))};
807 
808  auto tetGeom = MemoryManager<TetGeom>::AllocateSharedPtr(id, faces);
809  PopulateFaceToElMap(tetGeom, TetGeom::kNfaces);
810  geomMap[id] = tetGeom;
811 }
812 
813 template <>
814 void MeshGraphHDF5::ConstructGeomObject(
815  std::map<int, std::shared_ptr<PyrGeom>> &geomMap, int id, int *data,
816  CurveSharedPtr curve)
817 {
818  boost::ignore_unused(curve);
819  Geometry2DSharedPtr faces[5] = {
820  GetGeometry2D(data[0]), GetGeometry2D(data[1]), GetGeometry2D(data[2]),
821  GetGeometry2D(data[3]), GetGeometry2D(data[4])};
822 
823  auto pyrGeom = MemoryManager<PyrGeom>::AllocateSharedPtr(id, faces);
824  PopulateFaceToElMap(pyrGeom, PyrGeom::kNfaces);
825  geomMap[id] = pyrGeom;
826 }
827 
828 template <>
829 void MeshGraphHDF5::ConstructGeomObject(
830  std::map<int, std::shared_ptr<PrismGeom>> &geomMap, int id, int *data,
831  CurveSharedPtr curve)
832 {
833  boost::ignore_unused(curve);
834  Geometry2DSharedPtr faces[5] = {
835  GetGeometry2D(data[0]), GetGeometry2D(data[1]), GetGeometry2D(data[2]),
836  GetGeometry2D(data[3]), GetGeometry2D(data[4])};
837 
838  auto prismGeom = MemoryManager<PrismGeom>::AllocateSharedPtr(id, faces);
839  PopulateFaceToElMap(prismGeom, PrismGeom::kNfaces);
840  geomMap[id] = prismGeom;
841 }
842 
843 template <>
844 void MeshGraphHDF5::ConstructGeomObject(
845  std::map<int, std::shared_ptr<HexGeom>> &geomMap, int id, int *data,
846  CurveSharedPtr curve)
847 {
848  boost::ignore_unused(curve);
849  QuadGeomSharedPtr faces[6] = {
850  std::static_pointer_cast<QuadGeom>(GetGeometry2D(data[0])),
851  std::static_pointer_cast<QuadGeom>(GetGeometry2D(data[1])),
852  std::static_pointer_cast<QuadGeom>(GetGeometry2D(data[2])),
853  std::static_pointer_cast<QuadGeom>(GetGeometry2D(data[3])),
854  std::static_pointer_cast<QuadGeom>(GetGeometry2D(data[4])),
855  std::static_pointer_cast<QuadGeom>(GetGeometry2D(data[5]))};
856 
857  auto hexGeom = MemoryManager<HexGeom>::AllocateSharedPtr(id, faces);
858  PopulateFaceToElMap(hexGeom, HexGeom::kNfaces);
859  geomMap[id] = hexGeom;
860 }
861 
862 template <class T, typename DataType>
863 void MeshGraphHDF5::FillGeomMap(std::map<int, std::shared_ptr<T>> &geomMap,
864  const CurveMap &curveMap, std::vector<int> &ids,
865  std::vector<DataType> &geomData)
866 {
867  const int nGeomData = GetGeomDataDim(geomMap);
868  const int nRows = geomData.size() / nGeomData;
869  CurveSharedPtr empty;
870 
871  // Construct geometry object.
872  if (curveMap.size() > 0)
873  {
874  for (int i = 0, cnt = 0; i < nRows; i++, cnt += nGeomData)
875  {
876  auto cIt = curveMap.find(ids[i]);
877  ConstructGeomObject(geomMap, ids[i], &geomData[cnt],
878  cIt == curveMap.end() ? empty : cIt->second);
879  }
880  }
881  else
882  {
883  for (int i = 0, cnt = 0; i < nRows; i++, cnt += nGeomData)
884  {
885  ConstructGeomObject(geomMap, ids[i], &geomData[cnt], empty);
886  }
887  }
888 }
889 
890 template <class T, typename DataType>
891 void MeshGraphHDF5::ReadGeometryData(std::map<int, std::shared_ptr<T>> &geomMap,
892  std::string dataSet,
893  const std::unordered_set<int> &readIds,
894  std::vector<int> &ids,
895  std::vector<DataType> &geomData)
896 {
897  if (!m_mesh->ContainsDataSet(dataSet))
898  {
899  return;
900  }
901 
902  // Open mesh dataset
903  H5::DataSetSharedPtr data = m_mesh->OpenDataSet(dataSet);
904  H5::DataSpaceSharedPtr space = data->GetSpace();
905  vector<hsize_t> dims = space->GetDims();
906 
907  // Open metadata dataset
908  H5::DataSetSharedPtr mdata = m_maps->OpenDataSet(dataSet);
909  H5::DataSpaceSharedPtr mspace = mdata->GetSpace();
910  vector<hsize_t> mdims = mspace->GetDims();
911 
912  ASSERTL0(mdims[0] == dims[0], "map and data set lengths do not match");
913 
914  const int nGeomData = GetGeomDataDim(geomMap);
915 
916  // Read all IDs
917  vector<int> allIds;
918  mdata->Read(allIds, mspace);
919 
920  // Selective reading; clear data space range so that we can select certain
921  // rows from the datasets.
922  space->ClearRange();
923 
924  int i = 0;
925  std::vector<hsize_t> coords;
926  for (auto &id : allIds)
927  {
928  if (readIds.find(id) != readIds.end())
929  {
930  for (int j = 0; j < nGeomData; ++j)
931  {
932  coords.push_back(i);
933  coords.push_back(j);
934  }
935  ids.push_back(id);
936  }
937  ++i;
938  }
939 
940  space->SetSelection(coords.size() / 2, coords);
941 
942  // Read selected data.
943  data->Read(geomData, space, m_readPL);
944 }
945 
946 void MeshGraphHDF5::ReadCurveMap(CurveMap &curveMap, std::string dsName,
947  const std::unordered_set<int> &readIds)
948 {
949  // If dataset does not exist, exit.
950  if (!m_mesh->ContainsDataSet(dsName))
951  {
952  return;
953  }
954 
955  // Open up curve map data.
956  H5::DataSetSharedPtr curveData = m_mesh->OpenDataSet(dsName);
957  H5::DataSpaceSharedPtr curveSpace = curveData->GetSpace();
958 
959  // Open up ID data set.
960  H5::DataSetSharedPtr idData = m_maps->OpenDataSet(dsName);
961  H5::DataSpaceSharedPtr idSpace = idData->GetSpace();
962 
963  // Read all IDs and clear data space.
964  vector<int> ids, newIds;
965  idData->Read(ids, idSpace);
966  curveSpace->ClearRange();
967 
968  // Search IDs to figure out which curves to read.
969  vector<hsize_t> curveSel;
970 
971  int cnt = 0;
972  for (auto &id : ids)
973  {
974  if (readIds.find(id) != readIds.end())
975  {
976  curveSel.push_back(cnt);
977  curveSel.push_back(0);
978  curveSel.push_back(cnt);
979  curveSel.push_back(1);
980  curveSel.push_back(cnt);
981  curveSel.push_back(2);
982  newIds.push_back(id);
983  }
984 
985  ++cnt;
986  }
987 
988  // Check to see whether any processor will read anything
989  auto toRead = newIds.size();
990  m_session->GetComm()->AllReduce(toRead, LibUtilities::ReduceSum);
991 
992  if (toRead == 0)
993  {
994  return;
995  }
996 
997  // Now read curve map and read data.
998  vector<int> curveInfo;
999  curveSpace->SetSelection(curveSel.size() / 2, curveSel);
1000  curveData->Read(curveInfo, curveSpace, m_readPL);
1001 
1002  curveSel.clear();
1003 
1004  std::unordered_map<int, int> curvePtOffset;
1005 
1006  // Construct curves. We'll populate nodes in a minute!
1007  for (int i = 0, cnt = 0, cnt2 = 0; i < curveInfo.size() / 3; ++i, cnt += 3)
1008  {
1010  newIds[i], (LibUtilities::PointsType)curveInfo[cnt + 1]);
1011 
1012  curve->m_points.resize(curveInfo[cnt]);
1013 
1014  const int ptOffset = curveInfo[cnt + 2];
1015 
1016  for (int j = 0; j < curveInfo[cnt]; ++j)
1017  {
1018  // ptoffset gives us the row, multiply by 3 for number of
1019  // coordinates.
1020  curveSel.push_back(ptOffset + j);
1021  curveSel.push_back(0);
1022  curveSel.push_back(ptOffset + j);
1023  curveSel.push_back(1);
1024  curveSel.push_back(ptOffset + j);
1025  curveSel.push_back(2);
1026  }
1027 
1028  // Store the offset so we know to come back later on to fill in these
1029  // points.
1030  curvePtOffset[newIds[i]] = 3 * cnt2;
1031  cnt2 += curveInfo[cnt];
1032 
1033  curveMap[newIds[i]] = curve;
1034  }
1035 
1036  curveInfo.clear();
1037 
1038  // Open node data spacee.
1039  H5::DataSetSharedPtr nodeData = m_mesh->OpenDataSet("CURVE_NODES");
1040  H5::DataSpaceSharedPtr nodeSpace = nodeData->GetSpace();
1041 
1042  nodeSpace->ClearRange();
1043  nodeSpace->SetSelection(curveSel.size() / 2, curveSel);
1044 
1045  vector<NekDouble> nodeRawData;
1046  nodeData->Read(nodeRawData, nodeSpace, m_readPL);
1047 
1048  // Go back and populate data from nodes.
1049  for (auto &cIt : curvePtOffset)
1050  {
1051  CurveSharedPtr curve = curveMap[cIt.first];
1052 
1053  // Create nodes.
1054  int cnt = cIt.second;
1055  for (int i = 0; i < curve->m_points.size(); ++i, cnt += 3)
1056  {
1057  curve->m_points[i] = MemoryManager<PointGeom>::AllocateSharedPtr(
1058  0, m_spaceDimension, nodeRawData[cnt], nodeRawData[cnt + 1],
1059  nodeRawData[cnt + 2]);
1060  }
1061  }
1062 }
1063 
1064 void MeshGraphHDF5::ReadDomain()
1065 {
1066  map<int, CompositeSharedPtr> fullDomain;
1067  H5::DataSetSharedPtr dst = m_mesh->OpenDataSet("DOMAIN");
1068  H5::DataSpaceSharedPtr space = dst->GetSpace();
1069 
1070  vector<string> data;
1071  dst->ReadVectorString(data, space, m_readPL);
1072  GetCompositeList(data[0], fullDomain);
1073  m_domain[0] = fullDomain;
1074 }
1075 
1076 void MeshGraphHDF5::ReadComposites()
1077 {
1078  string nm = "COMPOSITE";
1079 
1080  H5::DataSetSharedPtr data = m_mesh->OpenDataSet(nm);
1081  H5::DataSpaceSharedPtr space = data->GetSpace();
1082  vector<hsize_t> dims = space->GetDims();
1083 
1084  vector<string> comps;
1085  data->ReadVectorString(comps, space);
1086 
1087  H5::DataSetSharedPtr mdata = m_maps->OpenDataSet(nm);
1088  H5::DataSpaceSharedPtr mspace = mdata->GetSpace();
1089  vector<hsize_t> mdims = mspace->GetDims();
1090 
1091  vector<int> ids;
1092  mdata->Read(ids, mspace);
1093 
1094  for (int i = 0; i < dims[0]; i++)
1095  {
1096  string compStr = comps[i];
1097 
1098  char type;
1099  istringstream strm(compStr);
1100 
1101  strm >> type;
1102 
1104 
1105  string::size_type indxBeg = compStr.find_first_of('[') + 1;
1106  string::size_type indxEnd = compStr.find_last_of(']') - 1;
1107 
1108  string indxStr = compStr.substr(indxBeg, indxEnd - indxBeg + 1);
1109  vector<unsigned int> seqVector;
1110 
1111  ParseUtils::GenerateSeqVector(indxStr, seqVector);
1112  m_compOrder[ids[i]] = seqVector;
1113 
1114  switch (type)
1115  {
1116  case 'V':
1117  for (auto &i : seqVector)
1118  {
1119  auto it = m_vertSet.find(i);
1120  if (it != m_vertSet.end())
1121  {
1122  comp->m_geomVec.push_back(it->second);
1123  }
1124  }
1125  break;
1126  case 'S':
1127  case 'E':
1128  for (auto &i : seqVector)
1129  {
1130  auto it = m_segGeoms.find(i);
1131  if (it != m_segGeoms.end())
1132  {
1133  comp->m_geomVec.push_back(it->second);
1134  }
1135  }
1136  break;
1137  case 'Q':
1138  for (auto &i : seqVector)
1139  {
1140  auto it = m_quadGeoms.find(i);
1141  if (it != m_quadGeoms.end())
1142  {
1143  comp->m_geomVec.push_back(it->second);
1144  }
1145  }
1146  break;
1147  case 'T':
1148  for (auto &i : seqVector)
1149  {
1150  auto it = m_triGeoms.find(i);
1151  if (it != m_triGeoms.end())
1152  {
1153  comp->m_geomVec.push_back(it->second);
1154  }
1155  }
1156  break;
1157  case 'F':
1158  for (auto &i : seqVector)
1159  {
1160  auto it1 = m_quadGeoms.find(i);
1161  if (it1 != m_quadGeoms.end())
1162  {
1163  comp->m_geomVec.push_back(it1->second);
1164  continue;
1165  }
1166  auto it2 = m_triGeoms.find(i);
1167  if (it2 != m_triGeoms.end())
1168  {
1169  comp->m_geomVec.push_back(it2->second);
1170  }
1171  }
1172  break;
1173  case 'A':
1174  for (auto &i : seqVector)
1175  {
1176  auto it = m_tetGeoms.find(i);
1177  if (it != m_tetGeoms.end())
1178  {
1179  comp->m_geomVec.push_back(it->second);
1180  }
1181  }
1182  break;
1183  case 'P':
1184  for (auto &i : seqVector)
1185  {
1186  auto it = m_pyrGeoms.find(i);
1187  if (it != m_pyrGeoms.end())
1188  {
1189  comp->m_geomVec.push_back(it->second);
1190  }
1191  }
1192  break;
1193  case 'R':
1194  for (auto &i : seqVector)
1195  {
1196  auto it = m_prismGeoms.find(i);
1197  if (it != m_prismGeoms.end())
1198  {
1199  comp->m_geomVec.push_back(it->second);
1200  }
1201  }
1202  break;
1203  case 'H':
1204  for (auto &i : seqVector)
1205  {
1206  auto it = m_hexGeoms.find(i);
1207  if (it != m_hexGeoms.end())
1208  {
1209  comp->m_geomVec.push_back(it->second);
1210  }
1211  }
1212  break;
1213  }
1214 
1215  if (comp->m_geomVec.size() > 0)
1216  {
1217  m_meshComposites[ids[i]] = comp;
1218  }
1219  }
1220 }
1221 
1222 CompositeDescriptor MeshGraphHDF5::CreateCompositeDescriptor(
1223  std::unordered_map<int, int> &id2row)
1224 {
1225  CompositeDescriptor ret;
1226 
1227  string nm = "COMPOSITE";
1228 
1229  H5::DataSetSharedPtr data = m_mesh->OpenDataSet(nm);
1230  H5::DataSpaceSharedPtr space = data->GetSpace();
1231  vector<hsize_t> dims = space->GetDims();
1232 
1233  vector<string> comps;
1234  data->ReadVectorString(comps, space);
1235 
1236  H5::DataSetSharedPtr mdata = m_maps->OpenDataSet(nm);
1237  H5::DataSpaceSharedPtr mspace = mdata->GetSpace();
1238  vector<hsize_t> mdims = mspace->GetDims();
1239 
1240  vector<int> ids;
1241  mdata->Read(ids, mspace);
1242 
1243  for (int i = 0; i < dims[0]; i++)
1244  {
1245  string compStr = comps[i];
1246 
1247  char type;
1248  istringstream strm(compStr);
1249 
1250  strm >> type;
1251 
1252  string::size_type indxBeg = compStr.find_first_of('[') + 1;
1253  string::size_type indxEnd = compStr.find_last_of(']') - 1;
1254 
1255  string indxStr = compStr.substr(indxBeg, indxEnd - indxBeg + 1);
1256  vector<unsigned int> seqVector;
1257  ParseUtils::GenerateSeqVector(indxStr, seqVector);
1258 
1259  LibUtilities::ShapeType shapeType;
1260 
1261  switch (type)
1262  {
1263  case 'V':
1264  shapeType = LibUtilities::ePoint;
1265  break;
1266  case 'S':
1267  case 'E':
1268  shapeType = LibUtilities::eSegment;
1269  break;
1270  case 'Q':
1271  case 'F':
1272  // Note that for HDF5, the composite descriptor is only used for
1273  // partitioning purposes so 'F' tag is not really going to be
1274  // critical in this context.
1275  shapeType = LibUtilities::eQuadrilateral;
1276  break;
1277  case 'T':
1278  shapeType = LibUtilities::eTriangle;
1279  break;
1280  case 'A':
1281  shapeType = LibUtilities::eTetrahedron;
1282  break;
1283  case 'P':
1284  shapeType = LibUtilities::ePyramid;
1285  break;
1286  case 'R':
1287  shapeType = LibUtilities::ePrism;
1288  break;
1289  case 'H':
1290  shapeType = LibUtilities::eHexahedron;
1291  break;
1292  }
1293 
1294  std::vector<int> filteredVector;
1295  for (auto &compElmt : seqVector)
1296  {
1297  if (id2row.find(compElmt) == id2row.end())
1298  {
1299  continue;
1300  }
1301 
1302  filteredVector.push_back(compElmt);
1303  }
1304 
1305  if (filteredVector.size() == 0)
1306  {
1307  continue;
1308  }
1309 
1310  ret[ids[i]] = std::make_pair(shapeType, filteredVector);
1311  }
1312 
1313  return ret;
1314 }
1315 
1316 template <class T, typename std::enable_if<T::kDim == 0, int>::type = 0>
1317 inline NekDouble GetGeomData(std::shared_ptr<T> &geom, int i)
1318 {
1319  return (*geom)(i);
1320 }
1321 
1322 template <class T, typename std::enable_if<T::kDim == 1, int>::type = 0>
1323 inline int GetGeomData(std::shared_ptr<T> &geom, int i)
1324 {
1325  return geom->GetVid(i);
1326 }
1327 
1328 template <class T, typename std::enable_if<T::kDim == 2, int>::type = 0>
1329 inline int GetGeomData(std::shared_ptr<T> &geom, int i)
1330 {
1331  return geom->GetEid(i);
1332 }
1333 
1334 template <class T, typename std::enable_if<T::kDim == 3, int>::type = 0>
1335 inline int GetGeomData(std::shared_ptr<T> &geom, int i)
1336 {
1337  return geom->GetFid(i);
1338 }
1339 
1340 template <class T>
1341 void MeshGraphHDF5::WriteGeometryMap(std::map<int, std::shared_ptr<T>> &geomMap,
1342  std::string datasetName)
1343 {
1344  typedef typename std::conditional<std::is_same<T, PointGeom>::value,
1345  NekDouble, int>::type DataType;
1346 
1347  const int nGeomData = GetGeomDataDim(geomMap);
1348  const size_t nGeom = geomMap.size();
1349 
1350  if (nGeom == 0)
1351  {
1352  return;
1353  }
1354 
1355  // Construct a map storing IDs
1356  vector<int> idMap(nGeom);
1357  vector<DataType> data(nGeom * nGeomData);
1358 
1359  int cnt1 = 0, cnt2 = 0;
1360  for (auto &it : geomMap)
1361  {
1362  idMap[cnt1++] = it.first;
1363 
1364  for (int j = 0; j < nGeomData; ++j)
1365  {
1366  data[cnt2 + j] = GetGeomData(it.second, j);
1367  }
1368 
1369  cnt2 += nGeomData;
1370  }
1371 
1372  vector<hsize_t> dims = {static_cast<hsize_t>(nGeom),
1373  static_cast<hsize_t>(nGeomData)};
1374  H5::DataTypeSharedPtr tp = H5::DataType::OfObject(data[0]);
1376  std::shared_ptr<H5::DataSpace>(new H5::DataSpace(dims));
1377  H5::DataSetSharedPtr dst = m_mesh->CreateDataSet(datasetName, tp, ds);
1378  dst->Write(data, ds);
1379 
1380  tp = H5::DataType::OfObject(idMap[0]);
1381  dims = {nGeom};
1382  ds = std::shared_ptr<H5::DataSpace>(new H5::DataSpace(dims));
1383  dst = m_maps->CreateDataSet(datasetName, tp, ds);
1384  dst->Write(idMap, ds);
1385 }
1386 
1387 void MeshGraphHDF5::WriteCurveMap(CurveMap &curves, std::string dsName,
1388  MeshCurvedPts &curvedPts, int &ptOffset,
1389  int &newIdx)
1390 {
1391  vector<int> data, map;
1392 
1393  // Compile curve data.
1394  for (auto &c : curves)
1395  {
1396  map.push_back(c.first);
1397  data.push_back(c.second->m_points.size());
1398  data.push_back(c.second->m_ptype);
1399  data.push_back(ptOffset);
1400 
1401  ptOffset += c.second->m_points.size();
1402 
1403  for (auto &pt : c.second->m_points)
1404  {
1405  MeshVertex v;
1406  v.id = newIdx;
1407  pt->GetCoords(v.x, v.y, v.z);
1408  curvedPts.pts.push_back(v);
1409  curvedPts.index.push_back(newIdx++);
1410  }
1411  }
1412 
1413  // Write data.
1414  vector<hsize_t> dims = {data.size() / 3, 3};
1415  H5::DataTypeSharedPtr tp = H5::DataType::OfObject(data[0]);
1417  std::shared_ptr<H5::DataSpace>(new H5::DataSpace(dims));
1418  H5::DataSetSharedPtr dst = m_mesh->CreateDataSet(dsName, tp, ds);
1419  dst->Write(data, ds);
1420 
1421  tp = H5::DataType::OfObject(map[0]);
1422  dims = {map.size()};
1423  ds = std::shared_ptr<H5::DataSpace>(new H5::DataSpace(dims));
1424  dst = m_maps->CreateDataSet(dsName, tp, ds);
1425  dst->Write(map, ds);
1426 }
1427 
1428 void MeshGraphHDF5::WriteCurvePoints(MeshCurvedPts &curvedPts)
1429 {
1430  vector<double> vertData(curvedPts.pts.size() * 3);
1431 
1432  int cnt = 0;
1433  for (auto &pt : curvedPts.pts)
1434  {
1435  vertData[cnt++] = pt.x;
1436  vertData[cnt++] = pt.y;
1437  vertData[cnt++] = pt.z;
1438  }
1439 
1440  vector<hsize_t> dims = {curvedPts.pts.size(), 3};
1441  H5::DataTypeSharedPtr tp = H5::DataType::OfObject(vertData[0]);
1443  std::shared_ptr<H5::DataSpace>(new H5::DataSpace(dims));
1444  H5::DataSetSharedPtr dst = m_mesh->CreateDataSet("CURVE_NODES", tp, ds);
1445  dst->Write(vertData, ds);
1446 }
1447 
1448 void MeshGraphHDF5::WriteComposites(CompositeMap &composites)
1449 {
1450  vector<string> comps;
1451 
1452  // dont need location map only a id map
1453  // will filter the composites per parition on read, its easier
1454  // composites do not need to be written in paralell.
1455  vector<int> c_map;
1456 
1457  for (auto &cIt : composites)
1458  {
1459  if (cIt.second->m_geomVec.size() == 0)
1460  {
1461  continue;
1462  }
1463 
1464  comps.push_back(GetCompositeString(cIt.second));
1465  c_map.push_back(cIt.first);
1466  }
1467 
1468  H5::DataTypeSharedPtr tp = H5::DataType::String();
1469  H5::DataSpaceSharedPtr ds = H5::DataSpace::OneD(comps.size());
1470  H5::DataSetSharedPtr dst = m_mesh->CreateDataSet("COMPOSITE", tp, ds);
1471  dst->WriteVectorString(comps, ds, tp);
1472 
1473  tp = H5::DataType::OfObject(c_map[0]);
1474  ds = H5::DataSpace::OneD(c_map.size());
1475  dst = m_maps->CreateDataSet("COMPOSITE", tp, ds);
1476  dst->Write(c_map, ds);
1477 }
1478 
1479 void MeshGraphHDF5::WriteDomain(map<int, CompositeMap> &domain)
1480 {
1481  vector<unsigned int> idxList;
1482  for (auto cIt = domain[0].begin(); cIt != domain[0].end(); ++cIt)
1483  {
1484  idxList.push_back(cIt->first);
1485  }
1486  stringstream domString;
1487  vector<string> doms;
1488  doms.push_back(ParseUtils::GenerateSeqString(idxList));
1489 
1490  H5::DataTypeSharedPtr tp = H5::DataType::String();
1491  H5::DataSpaceSharedPtr ds = H5::DataSpace::OneD(doms.size());
1492  H5::DataSetSharedPtr dst = m_mesh->CreateDataSet("DOMAIN", tp, ds);
1493  dst->WriteVectorString(doms, ds, tp);
1494 }
1495 
1496 void MeshGraphHDF5::WriteGeometry(
1497  std::string &outfilename, bool defaultExp,
1498  const LibUtilities::FieldMetaDataMap &metadata)
1499 {
1500  boost::ignore_unused(metadata);
1501 
1502  vector<string> tmp;
1503  boost::split(tmp, outfilename, boost::is_any_of("."));
1504  string filenameXml = tmp[0] + ".xml";
1505  string filenameHdf5 = tmp[0] + ".nekg";
1506 
1507  //////////////////
1508  // XML part
1509  //////////////////
1510 
1511  // Check to see if a xml of the same name exists
1512  // if might have boundary conditions etc, we will just alter the geometry
1513  // tag if needed
1514  TiXmlDocument *doc = new TiXmlDocument;
1515  TiXmlElement *root;
1516  TiXmlElement *geomTag;
1517 
1518  if (boost::filesystem::exists(filenameXml.c_str()))
1519  {
1520  ifstream file(filenameXml.c_str());
1521  file >> (*doc);
1522  TiXmlHandle docHandle(doc);
1523  root = docHandle.FirstChildElement("NEKTAR").Element();
1524  ASSERTL0(root, "Unable to find NEKTAR tag in file.");
1525  geomTag = root->FirstChildElement("GEOMETRY");
1526  defaultExp = false;
1527  }
1528  else
1529  {
1530  TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", "");
1531  doc->LinkEndChild(decl);
1532  root = new TiXmlElement("NEKTAR");
1533  doc->LinkEndChild(root);
1534 
1535  geomTag = new TiXmlElement("GEOMETRY");
1536  root->LinkEndChild(geomTag);
1537  }
1538 
1539  // Update attributes with dimensions.
1540  geomTag->SetAttribute("DIM", m_meshDimension);
1541  geomTag->SetAttribute("SPACE", m_spaceDimension);
1542  geomTag->SetAttribute("HDF5FILE", filenameHdf5);
1543 
1544  geomTag->Clear();
1545 
1546  if (defaultExp)
1547  {
1548  TiXmlElement *expTag = new TiXmlElement("EXPANSIONS");
1549 
1550  for (auto it = m_meshComposites.begin(); it != m_meshComposites.end();
1551  it++)
1552  {
1553  if (it->second->m_geomVec[0]->GetShapeDim() == m_meshDimension)
1554  {
1555  TiXmlElement *exp = new TiXmlElement("E");
1556  exp->SetAttribute(
1557  "COMPOSITE",
1558  "C[" + boost::lexical_cast<string>(it->first) + "]");
1559  exp->SetAttribute("NUMMODES", 4);
1560  exp->SetAttribute("TYPE", "MODIFIED");
1561  exp->SetAttribute("FIELDS", "u");
1562 
1563  expTag->LinkEndChild(exp);
1564  }
1565  }
1566  root->LinkEndChild(expTag);
1567  }
1568 
1569  doc->SaveFile(filenameXml);
1570 
1571  //////////////////
1572  // HDF5 part
1573  //////////////////
1574 
1575  // This is serial IO so we will just override any existing file.
1576  m_file = H5::File::Create(filenameHdf5, H5F_ACC_TRUNC);
1577  auto hdfRoot = m_file->CreateGroup("NEKTAR");
1578  auto hdfRoot2 = hdfRoot->CreateGroup("GEOMETRY");
1579 
1580  // Write format version.
1581  hdfRoot2->SetAttribute("FORMAT_VERSION", FORMAT_VERSION);
1582 
1583  // Create main groups.
1584  m_mesh = hdfRoot2->CreateGroup("MESH");
1585  m_maps = hdfRoot2->CreateGroup("MAPS");
1586 
1587  WriteGeometryMap(m_vertSet, "VERT");
1588  WriteGeometryMap(m_segGeoms, "SEG");
1589  if (m_meshDimension > 1)
1590  {
1591  WriteGeometryMap(m_triGeoms, "TRI");
1592  WriteGeometryMap(m_quadGeoms, "QUAD");
1593  }
1594  if (m_meshDimension > 2)
1595  {
1596  WriteGeometryMap(m_tetGeoms, "TET");
1597  WriteGeometryMap(m_pyrGeoms, "PYR");
1598  WriteGeometryMap(m_prismGeoms, "PRISM");
1599  WriteGeometryMap(m_hexGeoms, "HEX");
1600  }
1601 
1602  // Write curves
1603  int ptOffset = 0, newIdx = 0;
1604  MeshCurvedPts curvePts;
1605  WriteCurveMap(m_curvedEdges, "CURVE_EDGE", curvePts, ptOffset, newIdx);
1606  WriteCurveMap(m_curvedFaces, "CURVE_FACE", curvePts, ptOffset, newIdx);
1607  WriteCurvePoints(curvePts);
1608 
1609  // Write composites and domain.
1610  WriteComposites(m_meshComposites);
1611  WriteDomain(m_domain);
1612 }
1613 
1614 } // namespace SpatialDomains
1615 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define TIME_RESULT(verb, msg, timer)
HDF5 DataSpace wrapper.
Definition: H5.h:303
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::shared_ptr< DataSpace > DataSpaceSharedPtr
Definition: H5.h:83
std::shared_ptr< PList > PListSharedPtr
Definition: H5.h:97
std::shared_ptr< File > FileSharedPtr
Definition: H5.h:93
std::shared_ptr< DataType > DataTypeSharedPtr
Definition: H5.h:79
std::shared_ptr< Group > GroupSharedPtr
Definition: FieldIOHdf5.h:49
std::shared_ptr< DataSet > DataSetSharedPtr
Definition: H5.h:95
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:52
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< DomainRange > DomainRangeShPtr
Definition: DomainRange.h:66
CommFactory & GetCommFactory()
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
std::map< int, std::pair< LibUtilities::ShapeType, std::vector< int > > > CompositeDescriptor
Definition: MeshGraph.h:63
int GetGeomDataDim(std::map< int, std::shared_ptr< T >> &geomMap)
MeshPartitionFactory & GetMeshPartitionFactory()
std::pair< size_t, size_t > SplitWork(size_t vecsize, int rank, int nprocs)
Utility function to split a vector equally amongst a number of processors.
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:60
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:61
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:138
int GetGeomData(std::shared_ptr< T > &geom, int i)
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
std::shared_ptr< MeshPartition > MeshPartitionSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
MeshGraphFactory & GetMeshGraphFactory()
Definition: MeshGraph.cpp:73
std::shared_ptr< Composite > CompositeSharedPtr
Definition: MeshGraph.h:137
void UniqueValues(std::unordered_set< int > &unique, const std::vector< int > &input, T &...args)
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
std::vector< MeshVertex > pts
mapping to access pts value.
std::vector< NekInt64 > index
id of this Point set
std::vector< unsigned int > list