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