35#include <boost/algorithm/string.hpp>
46#define TIME_RESULT(verb, msg, timer) \
49 std::cout << " - " << msg << ": " << timer.TimePerTest(1) << "\n" \
97std::pair<size_t, size_t>
SplitWork(
size_t vecsize,
int rank,
int nprocs)
99 size_t div = vecsize / nprocs;
100 size_t rem = vecsize % nprocs;
103 return std::make_pair(rank * (div + 1), div + 1);
107 return std::make_pair((rank - rem) * div + rem * (div + 1), div);
111template <class T, typename std::enable_if<T::kDim == 0, int>::type = 0>
113 [[maybe_unused]] std::map<
int, std::shared_ptr<T>> &geomMap)
118template <class T, typename std::enable_if<T::kDim == 1, int>::type = 0>
120 [[maybe_unused]] std::map<
int, std::shared_ptr<T>> &geomMap)
125template <class T, typename std::enable_if<T::kDim == 2, int>::type = 0>
127 [[maybe_unused]] std::map<
int, std::shared_ptr<T>> &geomMap)
132template <class T, typename std::enable_if<T::kDim == 3, int>::type = 0>
134 [[maybe_unused]] std::map<
int, std::shared_ptr<T>> &geomMap)
140inline void UniqueValues([[maybe_unused]] std::unordered_set<int> &unique)
146 const std::vector<int> &input, T &...args)
158 "use-hdf5-node-comm",
"",
159 "Use a per-node communicator for HDF5 partitioning.");
172 const bool isRoot = comm->TreatAsRankZero();
179 session->InitSession();
184 m_xmlGeom = session->GetElement(
"NEKTAR/GEOMETRY");
185 TiXmlAttribute *attr =
m_xmlGeom->FirstAttribute();
186 int meshDimension = 3;
187 int spaceDimension = 3;
191 std::string attrName(attr->Name());
192 if (attrName ==
"DIM")
194 err = attr->QueryIntValue(&meshDimension);
195 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read mesh dimension.");
197 else if (attrName ==
"SPACE")
199 err = attr->QueryIntValue(&spaceDimension);
200 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read space dimension.");
202 else if (attrName ==
"PARTITION")
205 "PARTITION parameter should only be used in XML meshes");
207 else if (attrName ==
"HDF5FILE")
211 else if (attrName ==
"PARTITIONED")
214 "PARTITIONED parameter should only be used in XML meshes");
218 std::string errstr(
"Unknown attribute: ");
227 ASSERTL0(meshDimension <= spaceDimension,
228 "Mesh dimension greater than space dimension.");
237 if (commMesh->GetSize() > 1)
240 parallelProps = H5::PList::FileAccess();
241 parallelProps->SetMpio(commMesh);
243 m_readPL = H5::PList::DatasetXfer();
249 auto root =
m_file->OpenGroup(
"NEKTAR");
250 ASSERTL0(root,
"Cannot find NEKTAR group in HDF5 file.");
252 auto root2 = root->OpenGroup(
"GEOMETRY");
253 ASSERTL0(root2,
"Cannot find NEKTAR/GEOMETRY group in HDF5 file.");
258 for (; attrIt != attrEnd; ++attrIt)
260 if (*attrIt ==
"FORMAT_VERSION")
266 "Unable to determine Nektar++ geometry HDF5 file version.");
271 " is higher than supported in "
272 "this version of Nektar++");
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.");
289 std::vector<std::tuple<std::string, int, LibUtilities::ShapeType>>>
300 const bool verbRoot = isRoot && session->DefinesCmdLineArgument(
"verbose");
304 std::cout <<
"Reading HDF5 geometry..." << std::endl;
312 int innerRank = 0, innerSize = 1,
313 interRank = interComm->GetRowComm()->GetRank(),
314 interSize = interComm->GetRowComm()->GetSize();
316 if (session->DefinesCmdLineArgument(
"use-hdf5-node-comm"))
318 ASSERTL0(comm->GetSize() == commMesh->GetSize(),
319 "--use-hdf5-node-comm not available with Parallel-in-Time")
321 auto splitComm = comm->SplitCommNode();
322 innerComm = splitComm.first;
323 interComm = splitComm.second;
324 innerRank = innerComm->GetRank();
325 innerSize = innerComm->GetSize();
329 interRank = interComm->GetRank();
330 interSize = interComm->GetSize();
336 std::unordered_set<int> toRead;
345 std::cout <<
" - beginning partitioning" << std::endl;
355 const bool verbRoot2 =
356 isRoot && session->DefinesCmdLineArgument(
"verbose");
359 std::vector<int> ids;
363 std::vector<MeshEntity> elmts;
364 std::unordered_map<int, int> row2id, id2row;
374 auto parallelProps = H5::PList::FileAccess();
375 parallelProps->SetMpio(interComm);
378 readPL = H5::PList::DatasetXfer();
379 readPL->SetDxMpioCollective();
380 file = H5::File::Open(
m_hdf5Name, H5F_ACC_RDONLY, parallelProps);
382 auto root = file->OpenGroup(
"NEKTAR");
383 auto root2 = root->OpenGroup(
"GEOMETRY");
384 mesh = root2->OpenGroup(
"MESH");
385 maps = root2->OpenGroup(
"MAPS");
389 for (
auto &it : dataSets[meshDimension])
391 std::string ds = std::get<0>(it);
393 if (!mesh->ContainsDataSet(ds))
401 std::vector<hsize_t> dims = space->GetDims();
405 std::vector<hsize_t> mdims = mspace->GetDims();
410 std::vector<int> tmpElmts, tmpIds;
411 mdata->Read(tmpIds, mspace, readPL);
412 data->Read(tmpElmts, space, readPL);
414 const int nGeomData = std::get<1>(it);
416 for (
int i = 0, cnt = 0; i < tmpIds.size(); ++i, ++rowCount)
419 row2id[rowCount] = tmpIds[i];
420 id2row[tmpIds[i]] = row2id[rowCount];
424 e.
list = std::vector<unsigned int>(&tmpElmts[cnt],
425 &tmpElmts[cnt + nGeomData]);
431 interComm->GetRowComm()->Block();
438 size_t numElmt = elmts.size();
439 ASSERTL0(commMesh->GetSize() <= numElmt,
440 "This mesh has more processors than elements!");
442 auto elRange =
SplitWork(numElmt, interRank, interSize);
445 std::map<int, MeshEntity> partElmts;
446 std::unordered_set<int> facetIDs;
450 for (
int el = elRange.first; el < elRange.first + elRange.second;
455 partElmts[el] = elmt;
457 for (
auto &facet : elmt.
list)
459 facetIDs.insert(facet);
466 for (
int i = 0; i < numElmt; ++i)
469 if (i >= elRange.first && i < elRange.first + elRange.second)
478 for (
auto &eId : elmt.
list)
480 if (facetIDs.find(eId) != facetIDs.end())
490 partElmts[elmt.
id] = elmt;
497 std::string partitionerName =
498 commMesh->GetSize() > 1 ?
"PtScotch" :
"Scotch";
501 partitionerName =
"ParMetis";
503 if (session->DefinesCmdLineArgument(
"use-parmetis"))
505 partitionerName =
"ParMetis";
507 if (session->DefinesCmdLineArgument(
"use-ptscotch"))
509 partitionerName =
"PtScotch";
514 partitionerName, session, interComm, meshDimension, partElmts,
518 TIME_RESULT(verbRoot2,
" - partitioner setup", t2);
521 partitioner->PartitionMesh(interSize,
true,
false, nLocal);
528 std::vector<unsigned int> nodeElmts;
529 partitioner->GetElementIDs(interRank, nodeElmts);
534 std::map<int, MeshEntity> partElmts;
535 std::unordered_map<int, int> row2elmtid, elmtid2row;
542 for (
auto &elmtRow : nodeElmts)
544 row2elmtid[vcnt] = elmts[elmtRow].origId;
545 elmtid2row[elmts[elmtRow].origId] = vcnt;
548 partElmts[vcnt++] = elmt;
557 "Scotch", session, tmpComm, meshDimension, partElmts,
561 TIME_RESULT(verbRoot2,
" - inner partition setup", t2);
564 partitioner->PartitionMesh(innerSize,
true,
false, 0);
567 TIME_RESULT(verbRoot2,
" - inner partitioning", t2);
571 for (
int i = 1; i < innerSize; ++i)
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)
578 tmp[j] = row2elmtid[tmp[j]];
580 innerComm->Send(i, tmpsize);
581 innerComm->Send(i, tmp);
585 TIME_RESULT(verbRoot2,
" - inner partition scatter", t2);
587 std::vector<unsigned int> tmp;
588 partitioner->GetElementIDs(0, tmp);
590 for (
auto &tmpId : tmp)
592 toRead.insert(row2elmtid[tmpId]);
597 for (
auto &tmpId : nodeElmts)
599 toRead.insert(row2id[tmpId]);
608 innerComm->Recv(0, tmpSize);
609 std::vector<unsigned int> tmp(tmpSize);
610 innerComm->Recv(0, tmp);
612 for (
auto &tmpId : tmp)
614 toRead.insert(tmpId);
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;
636 auto &prismGeoms =
m_meshGraph->GetAllPrismGeoms();
641 if (meshDimension == 3)
651 UniqueValues(toRead, hexData, pyrData, prismData, tetData);
656 if (meshDimension >= 2)
669 if (meshDimension >= 1)
692 if (meshDimension >= 1)
696 for (
auto &edge : segIDs)
703 FillGeomMap(segGeoms, curvedEdges, segIDs, segData);
708 if (meshDimension >= 2)
712 for (
auto &face : triIDs)
716 for (
auto &face : quadIDs)
723 FillGeomMap(triGeoms, curvedFaces, triIDs, triData);
724 FillGeomMap(quadGeoms, curvedFaces, quadIDs, quadData);
729 if (meshDimension >= 3)
741 if (session->DefinesElement(
"NEKTAR/CONDITIONS"))
743 std::set<int> vBndRegionIdList;
744 TiXmlElement *vConditions =
745 new TiXmlElement(*session->GetElement(
"Nektar/Conditions"));
746 TiXmlElement *vBndRegions =
747 vConditions->FirstChildElement(
"BOUNDARYREGIONS");
754 auto &graph_bndRegOrder =
m_meshGraph->GetBndRegionOrdering();
755 vItem = vBndRegions->FirstChildElement();
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);
763 std::vector<unsigned int> vSeq;
766 int p = atoi(vItem->Attribute(
"ID"));
768 graph_bndRegOrder[
p] = vSeq;
769 vItem = vItem->NextSiblingElement();
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,
788 std::map<
int, std::shared_ptr<PointGeom>> &geomMap,
int id,
NekDouble *data,
792 m_meshGraph->GetSpaceDimension(),
id, data[0], data[1], data[2]);
797 std::map<
int, std::shared_ptr<SegGeom>> &geomMap,
int id,
int *data,
808 std::map<
int, std::shared_ptr<TriGeom>> &geomMap,
int id,
int *data,
819 std::map<
int, std::shared_ptr<QuadGeom>> &geomMap,
int id,
int *data,
830 std::map<
int, std::shared_ptr<TetGeom>> &geomMap,
int id,
int *data,
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]))};
841 geomMap[id] = tetGeom;
846 std::map<
int, std::shared_ptr<PyrGeom>> &geomMap,
int id,
int *data,
857 geomMap[id] = pyrGeom;
862 std::map<
int, std::shared_ptr<PrismGeom>> &geomMap,
int id,
int *data,
873 geomMap[id] = prismGeom;
878 std::map<
int, std::shared_ptr<HexGeom>> &geomMap,
int id,
int *data,
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>(
892 geomMap[id] = hexGeom;
895template <
class T,
typename DataType>
898 std::vector<int> &ids,
899 std::vector<DataType> &geomData)
902 const int nRows = geomData.size() / nGeomData;
906 if (curveMap.size() > 0)
908 for (
int i = 0, cnt = 0; i < nRows; i++, cnt += nGeomData)
910 auto cIt = curveMap.find(ids[i]);
912 cIt == curveMap.end() ? empty : cIt->second);
917 for (
int i = 0, cnt = 0; i < nRows; i++, cnt += nGeomData)
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)
930 if (!
m_mesh->ContainsDataSet(dataSet))
938 std::vector<hsize_t> dims = space->GetDims();
943 std::vector<hsize_t> mdims = mspace->GetDims();
945 ASSERTL0(mdims[0] == dims[0],
"map and data set lengths do not match");
950 std::vector<int> allIds;
951 mdata->Read(allIds, mspace);
958 std::vector<hsize_t> coords;
959 for (
auto &
id : allIds)
961 if (readIds.find(
id) != readIds.end())
963 for (
int j = 0; j < nGeomData; ++j)
973 space->SetSelection(coords.size() / 2, coords);
976 data->Read(geomData, space,
m_readPL);
980 const std::unordered_set<int> &readIds)
983 if (!
m_mesh->ContainsDataSet(dsName))
997 std::vector<int> ids, newIds;
998 idData->Read(ids, idSpace);
999 curveSpace->ClearRange();
1002 std::vector<hsize_t> curveSel;
1005 for (
auto &
id : ids)
1007 if (readIds.find(
id) != readIds.end())
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);
1022 auto toRead = newIds.size();
1023 m_session->GetComm()->GetRowComm()->AllReduce(toRead,
1032 std::vector<int> curveInfo;
1033 curveSpace->SetSelection(curveSel.size() / 2, curveSel);
1034 curveData->Read(curveInfo, curveSpace,
m_readPL);
1038 std::unordered_map<int, int> curvePtOffset;
1041 for (
int i = 0, cnt = 0, cnt2 = 0; i < curveInfo.size() / 3; ++i, cnt += 3)
1046 curve->m_points.resize(curveInfo[cnt]);
1048 const int ptOffset = curveInfo[cnt + 2];
1050 for (
int j = 0; j < curveInfo[cnt]; ++j)
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);
1064 curvePtOffset[newIds[i]] = 3 * cnt2;
1065 cnt2 += curveInfo[cnt];
1067 curveMap[newIds[i]] = curve;
1076 nodeSpace->ClearRange();
1077 nodeSpace->SetSelection(curveSel.size() / 2, curveSel);
1079 std::vector<NekDouble> nodeRawData;
1080 nodeData->Read(nodeRawData, nodeSpace,
m_readPL);
1083 for (
auto &cIt : curvePtOffset)
1088 int cnt = cIt.second;
1089 for (
int i = 0; i < curve->m_points.size(); ++i, cnt += 3)
1092 0,
m_meshGraph->GetSpaceDimension(), nodeRawData[cnt],
1093 nodeRawData[cnt + 1], nodeRawData[cnt + 2]);
1104 std::map<int, CompositeSharedPtr> fullDomain;
1108 std::vector<std::string> data;
1109 dst->ReadVectorString(data, space,
m_readPL);
1110 m_meshGraph->GetCompositeList(data[0], fullDomain);
1111 domain[0] = fullDomain;
1116 std::vector<CompositeMap> fullDomain;
1120 std::vector<std::string> data;
1121 dst->ReadVectorString(data, space,
m_readPL);
1122 for (
auto &dIt : data)
1125 m_meshGraph->GetCompositeList(dIt, fullDomain.back());
1131 std::vector<int> ids;
1132 mdata->Read(ids, mspace);
1134 for (
int i = 0; i < ids.size(); ++i)
1136 domain[ids[i]] = fullDomain[i];
1152 std::string nm =
"COMPOSITE";
1156 std::vector<hsize_t> dims = space->GetDims();
1158 std::vector<std::string> comps;
1159 data->ReadVectorString(comps, space);
1163 std::vector<hsize_t> mdims = mspace->GetDims();
1165 std::vector<int> ids;
1166 mdata->Read(ids, mspace);
1168 auto &graph_compOrder =
m_meshGraph->GetCompositeOrdering();
1169 for (
int i = 0; i < dims[0]; i++)
1171 std::string compStr = comps[i];
1174 std::istringstream strm(compStr);
1180 std::string::size_type indxBeg = compStr.find_first_of(
'[') + 1;
1181 std::string::size_type indxEnd = compStr.find_last_of(
']') - 1;
1183 std::string indxStr = compStr.substr(indxBeg, indxEnd - indxBeg + 1);
1184 std::vector<unsigned int> seqVector;
1188 graph_compOrder[ids[i]] = seqVector;
1193 for (
auto &i : seqVector)
1195 auto it = vertSet.find(i);
1196 if (it != vertSet.end())
1198 comp->m_geomVec.push_back(it->second);
1204 for (
auto &i : seqVector)
1206 auto it = segGeoms.find(i);
1207 if (it != segGeoms.end())
1209 comp->m_geomVec.push_back(it->second);
1214 for (
auto &i : seqVector)
1216 auto it = quadGeoms.find(i);
1217 if (it != quadGeoms.end())
1221 comp->m_geomVec.push_back(it->second);
1227 for (
auto &i : seqVector)
1229 auto it = triGeoms.find(i);
1230 if (it != triGeoms.end())
1234 comp->m_geomVec.push_back(it->second);
1240 for (
auto &i : seqVector)
1242 auto it1 = quadGeoms.find(i);
1243 if (it1 != quadGeoms.end())
1247 comp->m_geomVec.push_back(it1->second);
1250 auto it2 = triGeoms.find(i);
1251 if (it2 != triGeoms.end())
1255 comp->m_geomVec.push_back(it2->second);
1261 for (
auto &i : seqVector)
1263 auto it = tetGeoms.find(i);
1264 if (it != tetGeoms.end())
1268 comp->m_geomVec.push_back(it->second);
1274 for (
auto &i : seqVector)
1276 auto it = pyrGeoms.find(i);
1277 if (it != pyrGeoms.end())
1281 comp->m_geomVec.push_back(it->second);
1287 for (
auto &i : seqVector)
1289 auto it = prismGeoms.find(i);
1290 if (it != prismGeoms.end())
1294 comp->m_geomVec.push_back(it->second);
1300 for (
auto &i : seqVector)
1302 auto it = hexGeoms.find(i);
1303 if (it != hexGeoms.end())
1307 comp->m_geomVec.push_back(it->second);
1314 if (comp->m_geomVec.size() > 0)
1316 meshComposites[ids[i]] = comp;
1322 std::unordered_map<int, int> &id2row)
1326 std::string nm =
"COMPOSITE";
1330 std::vector<hsize_t> dims = space->GetDims();
1332 std::vector<std::string> comps;
1333 data->ReadVectorString(comps, space);
1337 std::vector<hsize_t> mdims = mspace->GetDims();
1339 std::vector<int> ids;
1340 mdata->Read(ids, mspace);
1342 for (
int i = 0; i < dims[0]; i++)
1344 std::string compStr = comps[i];
1347 std::istringstream strm(compStr);
1351 std::string::size_type indxBeg = compStr.find_first_of(
'[') + 1;
1352 std::string::size_type indxEnd = compStr.find_last_of(
']') - 1;
1354 std::string indxStr = compStr.substr(indxBeg, indxEnd - indxBeg + 1);
1355 std::vector<unsigned int> seqVector;
1395 std::vector<int> filteredVector;
1396 for (
auto &compElmt : seqVector)
1398 if (id2row.find(compElmt) == id2row.end())
1403 filteredVector.push_back(compElmt);
1406 if (filteredVector.size() == 0)
1411 ret[ids[i]] = std::make_pair(shapeType, filteredVector);
1417template <class T, typename std::enable_if<T::kDim == 0, int>::type = 0>
1423template <class T, typename std::enable_if<T::kDim == 1, int>::type = 0>
1426 return geom->GetVid(i);
1429template <class T, typename std::enable_if<T::kDim == 2, int>::type = 0>
1430inline int GetGeomData(std::shared_ptr<T> &geom,
int i)
1432 return geom->GetEid(i);
1435template <class T, typename std::enable_if<T::kDim == 3, int>::type = 0>
1436inline int GetGeomData(std::shared_ptr<T> &geom,
int i)
1438 return geom->GetFid(i);
1443 std::map<
int, std::shared_ptr<T>> &geomMap, std::string datasetName)
1445 typedef typename std::conditional<std::is_same_v<T, PointGeom>,
NekDouble,
1446 int>::type DataType;
1449 const size_t nGeom = geomMap.size();
1457 std::vector<int> idMap(nGeom);
1458 std::vector<DataType> data(nGeom * nGeomData);
1460 int cnt1 = 0, cnt2 = 0;
1461 for (
auto &it : geomMap)
1463 idMap[cnt1++] = it.first;
1465 for (
int j = 0; j < nGeomData; ++j)
1473 std::vector<hsize_t> dims = {
static_cast<hsize_t
>(nGeom),
1474 static_cast<hsize_t
>(nGeomData)};
1479 dst->Write(data, ds);
1481 tp = H5::DataType::OfObject(idMap[0]);
1483 ds = std::shared_ptr<H5::DataSpace>(
new H5::DataSpace(dims));
1484 dst =
m_maps->CreateDataSet(datasetName, tp, ds);
1485 dst->Write(idMap, ds);
1492 std::vector<int> data, map;
1495 for (
auto &c : curves)
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);
1502 ptOffset += c.second->m_points.size();
1504 for (
auto &pt : c.second->m_points)
1508 pt->GetCoords(v.
x, v.
y, v.
z);
1509 curvedPts.
pts.push_back(v);
1510 curvedPts.
index.push_back(newIdx++);
1515 std::vector<hsize_t> dims = {data.size() / 3, 3};
1520 dst->Write(data, ds);
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);
1531 std::vector<double> vertData(curvedPts.
pts.size() * 3);
1534 for (
auto &pt : curvedPts.
pts)
1536 vertData[cnt++] = pt.x;
1537 vertData[cnt++] = pt.y;
1538 vertData[cnt++] = pt.z;
1541 std::vector<hsize_t> dims = {curvedPts.
pts.size(), 3};
1546 dst->Write(vertData, ds);
1551 std::vector<std::string> comps;
1556 std::vector<int> c_map;
1558 for (
auto &cIt : composites)
1560 if (cIt.second->m_geomVec.size() == 0)
1566 c_map.push_back(cIt.first);
1572 dst->WriteVectorString(comps, ds, tp);
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);
1585 std::vector<int> d_map;
1586 std::vector<std::vector<unsigned int>> idxList;
1589 for (
auto &dIt : domain)
1591 idxList.push_back(std::vector<unsigned int>());
1592 for (
auto cIt = dIt.second.begin(); cIt != dIt.second.end(); ++cIt)
1594 idxList[cnt].push_back(cIt->first);
1598 d_map.push_back(dIt.first);
1601 std::stringstream domString;
1602 std::vector<std::string> doms;
1603 for (
auto &cIt : idxList)
1611 dst->WriteVectorString(doms, ds, tp);
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);
1620 const std::string &outfilename,
bool defaultExp,
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";
1648 TiXmlDocument *doc =
new TiXmlDocument;
1650 TiXmlElement *geomTag;
1652 if (fs::exists(filenameXml.c_str()))
1654 std::ifstream file(filenameXml.c_str());
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");
1664 TiXmlDeclaration *decl =
new TiXmlDeclaration(
"1.0",
"utf-8",
"");
1665 doc->LinkEndChild(decl);
1666 root =
new TiXmlElement(
"NEKTAR");
1667 doc->LinkEndChild(root);
1669 geomTag =
new TiXmlElement(
"GEOMETRY");
1670 root->LinkEndChild(geomTag);
1674 geomTag->SetAttribute(
"DIM",
m_meshGraph->GetMeshDimension());
1675 geomTag->SetAttribute(
"SPACE",
m_meshGraph->GetSpaceDimension());
1676 geomTag->SetAttribute(
"HDF5FILE", filenameHdf5);
1682 TiXmlElement *expTag =
new TiXmlElement(
"EXPANSIONS");
1684 for (
auto it = meshComposites.begin(); it != meshComposites.end(); it++)
1686 if (it->second->m_geomVec[0]->GetShapeDim() ==
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");
1696 expTag->LinkEndChild(exp);
1699 root->LinkEndChild(expTag);
1705 movement->WriteMovement(root);
1708 doc->SaveFile(filenameXml);
1715 m_file = H5::File::Create(filenameHdf5, H5F_ACC_TRUNC);
1716 auto hdfRoot =
m_file->CreateGroup(
"NEKTAR");
1717 auto hdfRoot2 = hdfRoot->CreateGroup(
"GEOMETRY");
1723 m_mesh = hdfRoot2->CreateGroup(
"MESH");
1724 m_maps = hdfRoot2->CreateGroup(
"MAPS");
1742 int ptOffset = 0, newIdx = 0;
1744 WriteCurveMap(curvedEdges,
"CURVE_EDGE", curvePts, ptOffset, newIdx);
1745 WriteCurveMap(curvedFaces,
"CURVE_FACE", curvePts, ptOffset, newIdx);
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define TIME_RESULT(verb, msg, timer)
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.
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.
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)
static std::string cmdSwitch
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)
unsigned int m_inFormatVersion
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)
static std::string className
void WriteComposites(CompositeMap &comps)
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()
BndRegionOrdering m_bndRegOrder
LibUtilities::SessionReaderSharedPtr m_session
std::string GetCompositeString(CompositeSharedPtr comp)
Returns a string representation of a composite.
CompositeDescriptor CreateCompositeDescriptor()
MeshGraphSharedPtr m_meshGraph
CompositeOrdering m_compOrder
std::shared_ptr< DataSpace > DataSpaceSharedPtr
std::shared_ptr< PList > PListSharedPtr
std::shared_ptr< File > FileSharedPtr
std::shared_ptr< DataType > DataTypeSharedPtr
std::shared_ptr< Group > GroupSharedPtr
std::shared_ptr< DataSet > DataSetSharedPtr
std::map< std::string, std::string > FieldMetaDataMap
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< DomainRange > DomainRangeShPtr
CommFactory & GetCommFactory()
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
std::map< int, TriGeomSharedPtr > TriGeomMap
NekDouble GetGeomData(std::shared_ptr< T > &geom, int i)
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
std::map< int, PyrGeomSharedPtr > PyrGeomMap
int GetGeomDataDim(std::map< int, std::shared_ptr< T > > &geomMap)
std::map< int, QuadGeomSharedPtr > QuadGeomMap
std::map< int, std::pair< LibUtilities::ShapeType, std::vector< int > > > CompositeDescriptor
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
std::map< int, SegGeomSharedPtr > SegGeomMap
std::shared_ptr< Curve > CurveSharedPtr
std::unordered_map< int, CurveSharedPtr > CurveMap
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::map< int, TetGeomSharedPtr > TetGeomMap
MeshGraphIOFactory & GetMeshGraphIOFactory()
std::map< int, PrismGeomSharedPtr > PrismGeomMap
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< MeshPartition > MeshPartitionSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
std::map< int, HexGeomSharedPtr > HexGeomMap
std::map< int, PointGeomSharedPtr > PointGeomMap
void UniqueValues(std::unordered_set< int > &unique)
std::map< int, CompositeSharedPtr > CompositeMap
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