35 #include <boost/algorithm/string.hpp>
36 #include <boost/core/ignore_unused.hpp>
37 #include <boost/filesystem.hpp>
39 #include <type_traits>
46 #define TIME_RESULT(verb, msg, timer) \
49 std::cout << " - " << msg << ": " << timer.TimePerTest(1) \
58 namespace SpatialDomains
63 const unsigned int MeshGraphHDF5::FORMAT_VERSION = 1;
65 std::string MeshGraphHDF5::className =
67 "IO with HDF5 geometry");
71 boost::ignore_unused(rng);
83 MeshGraph::FillGraph();
98 std::pair<size_t, size_t>
SplitWork(
size_t vecsize,
int rank,
int nprocs)
100 size_t div = vecsize / nprocs;
101 size_t rem = vecsize % nprocs;
104 return std::make_pair(rank * (div + 1), div + 1);
108 return std::make_pair((rank - rem) * div + rem * (div + 1), div);
112 template <class T, typename std::enable_if<T::kDim == 0, int>::type = 0>
115 boost::ignore_unused(geomMap);
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)
122 boost::ignore_unused(geomMap);
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)
129 boost::ignore_unused(geomMap);
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)
136 boost::ignore_unused(geomMap);
140 template <
class... T>
inline void UniqueValues(std::unordered_set<int> &unique)
142 boost::ignore_unused(unique);
145 template <
class... T>
147 const std::vector<int> &input, T &...args)
157 std::string MeshGraphHDF5::cmdSwitch =
158 LibUtilities::SessionReader::RegisterCmdLineFlag(
159 "use-hdf5-node-comm",
"",
160 "Use a per-node communicator for HDF5 partitioning.");
171 int rank = comm->GetRank(), nproc = comm->GetSize();
178 session->InitSession();
183 m_xmlGeom = m_session->GetElement(
"NEKTAR/GEOMETRY");
184 TiXmlAttribute *attr = m_xmlGeom->FirstAttribute();
185 m_meshPartitioned =
true;
187 m_spaceDimension = 3;
191 std::string attrName(attr->Name());
192 if (attrName ==
"DIM")
194 err = attr->QueryIntValue(&m_meshDimension);
195 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read mesh dimension.");
197 else if (attrName ==
"SPACE")
199 err = attr->QueryIntValue(&m_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")
209 m_hdf5Name = attr->Value();
210 ASSERTL1(err == TIXML_SUCCESS,
"Unable to read hdf5 name.");
212 else if (attrName ==
"PARTITIONED")
215 "PARTITIONED parameter should only be used in XML meshes");
219 std::string errstr(
"Unknown attribute: ");
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.");
233 m_readPL = H5::PList::Default();
238 parallelProps = H5::PList::FileAccess();
239 parallelProps->SetMpio(comm);
241 m_readPL = H5::PList::DatasetXfer();
242 m_readPL->SetDxMpioCollective();
245 m_file = H5::File::Open(m_hdf5Name, H5F_ACC_RDONLY, parallelProps);
247 auto root = m_file->OpenGroup(
"NEKTAR");
248 ASSERTL0(root,
"Cannot find NEKTAR group in HDF5 file.");
250 auto root2 = root->OpenGroup(
"GEOMETRY");
251 ASSERTL0(root2,
"Cannot find NEKTAR/GEOMETRY group in HDF5 file.");
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)
259 if (*attrIt ==
"FORMAT_VERSION")
265 "Unable to determine Nektar++ geometry HDF5 file version.");
266 root2->GetAttribute(
"FORMAT_VERSION", formatVersion);
268 ASSERTL0(formatVersion <= FORMAT_VERSION,
269 "File format in " + m_hdf5Name +
270 " is higher than supported in "
271 "this version of Nektar++");
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.");
280 std::vector<std::tuple<std::string, int, LibUtilities::ShapeType>>>
291 bool verbRoot = rank == 0 && m_session->DefinesCmdLineArgument(
"verbose");
295 std::cout <<
"Reading HDF5 geometry..." << std::endl;
303 int innerRank = 0, innerSize = 1, interRank = rank, interSize = nproc;
305 if (session->DefinesCmdLineArgument(
"use-hdf5-node-comm"))
307 auto splitComm = comm->SplitCommNode();
308 innerComm = splitComm.first;
309 interComm = splitComm.second;
310 innerRank = innerComm->GetRank();
311 innerSize = innerComm->GetSize();
315 interRank = interComm->GetRank();
316 interSize = interComm->GetSize();
322 std::unordered_set<int> toRead;
331 std::cout <<
" - beginning partitioning" << std::endl;
342 session->DefinesCmdLineArgument(
"verbose") && interRank == 0;
345 std::vector<int> ids;
349 std::vector<MeshEntity> elmts;
350 std::unordered_map<int, int> row2id, id2row;
360 auto parallelProps = H5::PList::FileAccess();
361 parallelProps->SetMpio(interComm);
364 readPL = H5::PList::DatasetXfer();
365 readPL->SetDxMpioCollective();
366 file = H5::File::Open(m_hdf5Name, H5F_ACC_RDONLY, parallelProps);
368 auto root = file->OpenGroup(
"NEKTAR");
369 auto root2 = root->OpenGroup(
"GEOMETRY");
370 mesh = root2->OpenGroup(
"MESH");
371 maps = root2->OpenGroup(
"MAPS");
375 for (
auto &it : dataSets[m_meshDimension])
377 std::string ds = std::get<0>(it);
379 if (!mesh->ContainsDataSet(ds))
387 vector<hsize_t> dims = space->GetDims();
391 vector<hsize_t> mdims = mspace->GetDims();
396 vector<int> tmpElmts, tmpIds;
397 mdata->Read(tmpIds, mspace, readPL);
398 data->Read(tmpElmts, space, readPL);
400 const int nGeomData = std::get<1>(it);
402 for (
int i = 0, cnt = 0; i < tmpIds.size(); ++i, ++rowCount)
405 row2id[rowCount] = tmpIds[i];
406 id2row[tmpIds[i]] = row2id[rowCount];
410 e.
list = std::vector<unsigned int>(&tmpElmts[cnt],
411 &tmpElmts[cnt + nGeomData]);
424 size_t numElmt = elmts.size();
426 "This mesh has more processors than elements!");
428 auto elRange =
SplitWork(numElmt, interRank, interSize);
431 std::map<int, MeshEntity> partElmts;
432 std::unordered_set<int> facetIDs;
436 for (
int el = elRange.first; el < elRange.first + elRange.second;
441 partElmts[el] = elmt;
443 for (
auto &facet : elmt.
list)
445 facetIDs.insert(facet);
452 for (
int i = 0; i < numElmt; ++i)
455 if (i >= elRange.first && i < elRange.first + elRange.second)
464 for (
auto &eId : elmt.
list)
466 if (facetIDs.find(eId) != facetIDs.end())
476 partElmts[elmt.
id] = elmt;
483 string partitionerName = nproc > 1 ?
"PtScotch" :
"Scotch";
486 partitionerName =
"ParMetis";
488 if (session->DefinesCmdLineArgument(
"use-parmetis"))
490 partitionerName =
"ParMetis";
492 if (session->DefinesCmdLineArgument(
"use-ptscotch"))
494 partitionerName =
"PtScotch";
499 partitionerName, session, interComm, m_meshDimension, partElmts,
500 CreateCompositeDescriptor(id2row));
503 TIME_RESULT(verbRoot2,
" - partitioner setup", t2);
506 partitioner->PartitionMesh(interSize,
true,
false, nLocal);
513 std::vector<unsigned int> nodeElmts;
514 partitioner->GetElementIDs(interRank, nodeElmts);
519 std::map<int, MeshEntity> partElmts;
520 std::unordered_map<int, int> row2elmtid, elmtid2row;
527 for (
auto &elmtRow : nodeElmts)
529 row2elmtid[vcnt] = elmts[elmtRow].origId;
530 elmtid2row[elmts[elmtRow].origId] = vcnt;
533 partElmts[vcnt++] = elmt;
542 "Scotch", session, tmpComm, m_meshDimension, partElmts,
543 CreateCompositeDescriptor(elmtid2row));
546 TIME_RESULT(verbRoot2,
" - inner partition setup", t2);
549 partitioner->PartitionMesh(innerSize,
true,
false, 0);
552 TIME_RESULT(verbRoot2,
" - inner partitioning", t2);
556 for (
int i = 1; i < innerSize; ++i)
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)
563 tmp[j] = row2elmtid[tmp[j]];
565 innerComm->Send(i, tmpsize);
566 innerComm->Send(i, tmp);
570 TIME_RESULT(verbRoot2,
" - inner partition scatter", t2);
572 std::vector<unsigned int> tmp;
573 partitioner->GetElementIDs(0, tmp);
575 for (
auto &tmpId : tmp)
577 toRead.insert(row2elmtid[tmpId]);
582 for (
auto &tmpId : nodeElmts)
584 toRead.insert(row2id[tmpId]);
593 innerComm->Recv(0, tmpSize);
594 std::vector<unsigned int> tmp(tmpSize);
595 innerComm->Recv(0, tmp);
597 for (
auto &tmpId : tmp)
599 toRead.insert(tmpId);
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;
615 if (m_meshDimension == 3)
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);
625 UniqueValues(toRead, hexData, pyrData, prismData, tetData);
630 if (m_meshDimension >= 2)
634 ReadGeometryData(m_triGeoms,
"TRI", toRead, triIDs, triData);
635 ReadGeometryData(m_quadGeoms,
"QUAD", toRead, quadIDs, quadData);
643 if (m_meshDimension >= 1)
647 ReadGeometryData(m_segGeoms,
"SEG", toRead, segIDs, segData);
656 ReadGeometryData(m_vertSet,
"VERT", toRead, vertIDs, vertData);
662 FillGeomMap(m_vertSet,
CurveMap(), vertIDs, vertData);
666 if (m_meshDimension >= 1)
670 for (
auto &edge : segIDs)
674 ReadCurveMap(m_curvedEdges,
"CURVE_EDGE", toRead);
677 FillGeomMap(m_segGeoms, m_curvedEdges, segIDs, segData);
682 if (m_meshDimension >= 2)
686 for (
auto &face : triIDs)
690 for (
auto &face : quadIDs)
694 ReadCurveMap(m_curvedFaces,
"CURVE_FACE", toRead);
697 FillGeomMap(m_triGeoms, m_curvedFaces, triIDs, triData);
698 FillGeomMap(m_quadGeoms, m_curvedFaces, quadIDs, quadData);
703 if (m_meshDimension >= 3)
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);
715 if (m_session->DefinesElement(
"NEKTAR/CONDITIONS"))
717 std::set<int> vBndRegionIdList;
718 TiXmlElement *vConditions =
719 new TiXmlElement(*m_session->GetElement(
"Nektar/Conditions"));
720 TiXmlElement *vBndRegions =
721 vConditions->FirstChildElement(
"BOUNDARYREGIONS");
726 vItem = vBndRegions->FirstChildElement();
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);
734 std::vector<unsigned int> vSeq;
735 ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
737 int p = atoi(vItem->Attribute(
"ID"));
738 m_bndRegOrder[
p] = vSeq;
739 vItem = vItem->NextSiblingElement();
748 template <
class T,
typename DataType>
749 void MeshGraphHDF5::ConstructGeomObject(
750 std::map<
int, std::shared_ptr<T>> &geomMap,
int id, DataType *data,
753 boost::ignore_unused(geomMap,
id, data, curve);
757 void MeshGraphHDF5::ConstructGeomObject(
758 std::map<
int, std::shared_ptr<PointGeom>> &geomMap,
int id,
NekDouble *data,
761 boost::ignore_unused(curve);
763 m_spaceDimension,
id, data[0], data[1], data[2]);
767 void MeshGraphHDF5::ConstructGeomObject(
768 std::map<
int, std::shared_ptr<SegGeom>> &geomMap,
int id,
int *data,
773 id, m_spaceDimension, pts, curve);
777 void MeshGraphHDF5::ConstructGeomObject(
778 std::map<
int, std::shared_ptr<TriGeom>> &geomMap,
int id,
int *data,
782 GetSegGeom(data[2])};
787 void MeshGraphHDF5::ConstructGeomObject(
788 std::map<
int, std::shared_ptr<QuadGeom>> &geomMap,
int id,
int *data,
792 GetSegGeom(data[2]), GetSegGeom(data[3])};
797 void MeshGraphHDF5::ConstructGeomObject(
798 std::map<
int, std::shared_ptr<TetGeom>> &geomMap,
int id,
int *data,
801 boost::ignore_unused(curve);
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]))};
809 PopulateFaceToElMap(tetGeom, TetGeom::kNfaces);
810 geomMap[id] = tetGeom;
814 void MeshGraphHDF5::ConstructGeomObject(
815 std::map<
int, std::shared_ptr<PyrGeom>> &geomMap,
int id,
int *data,
818 boost::ignore_unused(curve);
820 GetGeometry2D(data[0]), GetGeometry2D(data[1]), GetGeometry2D(data[2]),
821 GetGeometry2D(data[3]), GetGeometry2D(data[4])};
824 PopulateFaceToElMap(pyrGeom, PyrGeom::kNfaces);
825 geomMap[id] = pyrGeom;
829 void MeshGraphHDF5::ConstructGeomObject(
830 std::map<
int, std::shared_ptr<PrismGeom>> &geomMap,
int id,
int *data,
833 boost::ignore_unused(curve);
835 GetGeometry2D(data[0]), GetGeometry2D(data[1]), GetGeometry2D(data[2]),
836 GetGeometry2D(data[3]), GetGeometry2D(data[4])};
839 PopulateFaceToElMap(prismGeom, PrismGeom::kNfaces);
840 geomMap[id] = prismGeom;
844 void MeshGraphHDF5::ConstructGeomObject(
845 std::map<
int, std::shared_ptr<HexGeom>> &geomMap,
int id,
int *data,
848 boost::ignore_unused(curve);
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]))};
858 PopulateFaceToElMap(hexGeom, HexGeom::kNfaces);
859 geomMap[id] = hexGeom;
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)
868 const int nRows = geomData.size() / nGeomData;
872 if (curveMap.size() > 0)
874 for (
int i = 0, cnt = 0; i < nRows; i++, cnt += nGeomData)
876 auto cIt = curveMap.find(ids[i]);
877 ConstructGeomObject(geomMap, ids[i], &geomData[cnt],
878 cIt == curveMap.end() ? empty : cIt->second);
883 for (
int i = 0, cnt = 0; i < nRows; i++, cnt += nGeomData)
885 ConstructGeomObject(geomMap, ids[i], &geomData[cnt], empty);
890 template <
class T,
typename DataType>
891 void MeshGraphHDF5::ReadGeometryData(std::map<
int, std::shared_ptr<T>> &geomMap,
893 const std::unordered_set<int> &readIds,
894 std::vector<int> &ids,
895 std::vector<DataType> &geomData)
897 if (!m_mesh->ContainsDataSet(dataSet))
905 vector<hsize_t> dims = space->GetDims();
910 vector<hsize_t> mdims = mspace->GetDims();
912 ASSERTL0(mdims[0] == dims[0],
"map and data set lengths do not match");
918 mdata->Read(allIds, mspace);
925 std::vector<hsize_t> coords;
926 for (
auto &
id : allIds)
928 if (readIds.find(
id) != readIds.end())
930 for (
int j = 0; j < nGeomData; ++j)
940 space->SetSelection(coords.size() / 2, coords);
943 data->Read(geomData, space, m_readPL);
946 void MeshGraphHDF5::ReadCurveMap(
CurveMap &curveMap, std::string dsName,
947 const std::unordered_set<int> &readIds)
950 if (!m_mesh->ContainsDataSet(dsName))
964 vector<int> ids, newIds;
965 idData->Read(ids, idSpace);
966 curveSpace->ClearRange();
969 vector<hsize_t> curveSel;
974 if (readIds.find(
id) != readIds.end())
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);
989 auto toRead = newIds.size();
998 vector<int> curveInfo;
999 curveSpace->SetSelection(curveSel.size() / 2, curveSel);
1000 curveData->Read(curveInfo, curveSpace, m_readPL);
1004 std::unordered_map<int, int> curvePtOffset;
1007 for (
int i = 0, cnt = 0, cnt2 = 0; i < curveInfo.size() / 3; ++i, cnt += 3)
1012 curve->m_points.resize(curveInfo[cnt]);
1014 const int ptOffset = curveInfo[cnt + 2];
1016 for (
int j = 0; j < curveInfo[cnt]; ++j)
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);
1030 curvePtOffset[newIds[i]] = 3 * cnt2;
1031 cnt2 += curveInfo[cnt];
1033 curveMap[newIds[i]] = curve;
1042 nodeSpace->ClearRange();
1043 nodeSpace->SetSelection(curveSel.size() / 2, curveSel);
1045 vector<NekDouble> nodeRawData;
1046 nodeData->Read(nodeRawData, nodeSpace, m_readPL);
1049 for (
auto &cIt : curvePtOffset)
1054 int cnt = cIt.second;
1055 for (
int i = 0; i < curve->m_points.size(); ++i, cnt += 3)
1058 0, m_spaceDimension, nodeRawData[cnt], nodeRawData[cnt + 1],
1059 nodeRawData[cnt + 2]);
1064 void MeshGraphHDF5::ReadDomain()
1066 map<int, CompositeSharedPtr> fullDomain;
1070 vector<string> data;
1071 dst->ReadVectorString(data, space, m_readPL);
1072 GetCompositeList(data[0], fullDomain);
1073 m_domain[0] = fullDomain;
1076 void MeshGraphHDF5::ReadComposites()
1078 string nm =
"COMPOSITE";
1082 vector<hsize_t> dims = space->GetDims();
1084 vector<string> comps;
1085 data->ReadVectorString(comps, space);
1089 vector<hsize_t> mdims = mspace->GetDims();
1092 mdata->Read(ids, mspace);
1094 for (
int i = 0; i < dims[0]; i++)
1096 string compStr = comps[i];
1099 istringstream strm(compStr);
1105 string::size_type indxBeg = compStr.find_first_of(
'[') + 1;
1106 string::size_type indxEnd = compStr.find_last_of(
']') - 1;
1108 string indxStr = compStr.substr(indxBeg, indxEnd - indxBeg + 1);
1109 vector<unsigned int> seqVector;
1111 ParseUtils::GenerateSeqVector(indxStr, seqVector);
1112 m_compOrder[ids[i]] = seqVector;
1117 for (
auto &i : seqVector)
1119 auto it = m_vertSet.find(i);
1120 if (it != m_vertSet.end())
1122 comp->m_geomVec.push_back(it->second);
1128 for (
auto &i : seqVector)
1130 auto it = m_segGeoms.find(i);
1131 if (it != m_segGeoms.end())
1133 comp->m_geomVec.push_back(it->second);
1138 for (
auto &i : seqVector)
1140 auto it = m_quadGeoms.find(i);
1141 if (it != m_quadGeoms.end())
1143 comp->m_geomVec.push_back(it->second);
1148 for (
auto &i : seqVector)
1150 auto it = m_triGeoms.find(i);
1151 if (it != m_triGeoms.end())
1153 comp->m_geomVec.push_back(it->second);
1158 for (
auto &i : seqVector)
1160 auto it1 = m_quadGeoms.find(i);
1161 if (it1 != m_quadGeoms.end())
1163 comp->m_geomVec.push_back(it1->second);
1166 auto it2 = m_triGeoms.find(i);
1167 if (it2 != m_triGeoms.end())
1169 comp->m_geomVec.push_back(it2->second);
1174 for (
auto &i : seqVector)
1176 auto it = m_tetGeoms.find(i);
1177 if (it != m_tetGeoms.end())
1179 comp->m_geomVec.push_back(it->second);
1184 for (
auto &i : seqVector)
1186 auto it = m_pyrGeoms.find(i);
1187 if (it != m_pyrGeoms.end())
1189 comp->m_geomVec.push_back(it->second);
1194 for (
auto &i : seqVector)
1196 auto it = m_prismGeoms.find(i);
1197 if (it != m_prismGeoms.end())
1199 comp->m_geomVec.push_back(it->second);
1204 for (
auto &i : seqVector)
1206 auto it = m_hexGeoms.find(i);
1207 if (it != m_hexGeoms.end())
1209 comp->m_geomVec.push_back(it->second);
1215 if (comp->m_geomVec.size() > 0)
1217 m_meshComposites[ids[i]] = comp;
1223 std::unordered_map<int, int> &id2row)
1227 string nm =
"COMPOSITE";
1231 vector<hsize_t> dims = space->GetDims();
1233 vector<string> comps;
1234 data->ReadVectorString(comps, space);
1238 vector<hsize_t> mdims = mspace->GetDims();
1241 mdata->Read(ids, mspace);
1243 for (
int i = 0; i < dims[0]; i++)
1245 string compStr = comps[i];
1248 istringstream strm(compStr);
1252 string::size_type indxBeg = compStr.find_first_of(
'[') + 1;
1253 string::size_type indxEnd = compStr.find_last_of(
']') - 1;
1255 string indxStr = compStr.substr(indxBeg, indxEnd - indxBeg + 1);
1256 vector<unsigned int> seqVector;
1257 ParseUtils::GenerateSeqVector(indxStr, seqVector);
1294 std::vector<int> filteredVector;
1295 for (
auto &compElmt : seqVector)
1297 if (id2row.find(compElmt) == id2row.end())
1302 filteredVector.push_back(compElmt);
1305 if (filteredVector.size() == 0)
1310 ret[ids[i]] = std::make_pair(shapeType, filteredVector);
1316 template <class T, typename std::enable_if<T::kDim == 0, int>::type = 0>
1322 template <class T, typename std::enable_if<T::kDim == 1, int>::type = 0>
1325 return geom->GetVid(i);
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)
1331 return geom->GetEid(i);
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)
1337 return geom->GetFid(i);
1341 void MeshGraphHDF5::WriteGeometryMap(std::map<
int, std::shared_ptr<T>> &geomMap,
1342 std::string datasetName)
1344 typedef typename std::conditional<std::is_same<T, PointGeom>::value,
1348 const size_t nGeom = geomMap.size();
1356 vector<int> idMap(nGeom);
1357 vector<DataType> data(nGeom * nGeomData);
1359 int cnt1 = 0, cnt2 = 0;
1360 for (
auto &it : geomMap)
1362 idMap[cnt1++] = it.first;
1364 for (
int j = 0; j < nGeomData; ++j)
1372 vector<hsize_t> dims = {
static_cast<hsize_t
>(nGeom),
1373 static_cast<hsize_t
>(nGeomData)};
1378 dst->Write(data, ds);
1380 tp = H5::DataType::OfObject(idMap[0]);
1382 ds = std::shared_ptr<H5::DataSpace>(
new H5::DataSpace(dims));
1383 dst = m_maps->CreateDataSet(datasetName, tp, ds);
1384 dst->Write(idMap, ds);
1387 void MeshGraphHDF5::WriteCurveMap(
CurveMap &curves, std::string dsName,
1391 vector<int> data, map;
1394 for (
auto &c : curves)
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);
1401 ptOffset += c.second->m_points.size();
1403 for (
auto &pt : c.second->m_points)
1407 pt->GetCoords(v.
x, v.
y, v.
z);
1408 curvedPts.
pts.push_back(v);
1409 curvedPts.
index.push_back(newIdx++);
1414 vector<hsize_t> dims = {data.size() / 3, 3};
1419 dst->Write(data, ds);
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);
1430 vector<double> vertData(curvedPts.
pts.size() * 3);
1433 for (
auto &pt : curvedPts.
pts)
1435 vertData[cnt++] = pt.x;
1436 vertData[cnt++] = pt.y;
1437 vertData[cnt++] = pt.z;
1440 vector<hsize_t> dims = {curvedPts.
pts.size(), 3};
1445 dst->Write(vertData, ds);
1450 vector<string> comps;
1457 for (
auto &cIt : composites)
1459 if (cIt.second->m_geomVec.size() == 0)
1464 comps.push_back(GetCompositeString(cIt.second));
1465 c_map.push_back(cIt.first);
1471 dst->WriteVectorString(comps, ds, tp);
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);
1479 void MeshGraphHDF5::WriteDomain(map<int, CompositeMap> &domain)
1481 vector<unsigned int> idxList;
1482 for (
auto cIt = domain[0].begin(); cIt != domain[0].end(); ++cIt)
1484 idxList.push_back(cIt->first);
1486 stringstream domString;
1487 vector<string> doms;
1488 doms.push_back(ParseUtils::GenerateSeqString(idxList));
1493 dst->WriteVectorString(doms, ds, tp);
1496 void MeshGraphHDF5::WriteGeometry(
1497 std::string &outfilename,
bool defaultExp,
1500 boost::ignore_unused(metadata);
1503 boost::split(tmp, outfilename, boost::is_any_of(
"."));
1504 string filenameXml = tmp[0] +
".xml";
1505 string filenameHdf5 = tmp[0] +
".nekg";
1514 TiXmlDocument *doc =
new TiXmlDocument;
1516 TiXmlElement *geomTag;
1518 if (boost::filesystem::exists(filenameXml.c_str()))
1520 ifstream file(filenameXml.c_str());
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");
1530 TiXmlDeclaration *decl =
new TiXmlDeclaration(
"1.0",
"utf-8",
"");
1531 doc->LinkEndChild(decl);
1532 root =
new TiXmlElement(
"NEKTAR");
1533 doc->LinkEndChild(root);
1535 geomTag =
new TiXmlElement(
"GEOMETRY");
1536 root->LinkEndChild(geomTag);
1540 geomTag->SetAttribute(
"DIM", m_meshDimension);
1541 geomTag->SetAttribute(
"SPACE", m_spaceDimension);
1542 geomTag->SetAttribute(
"HDF5FILE", filenameHdf5);
1548 TiXmlElement *expTag =
new TiXmlElement(
"EXPANSIONS");
1550 for (
auto it = m_meshComposites.begin(); it != m_meshComposites.end();
1553 if (it->second->m_geomVec[0]->GetShapeDim() == m_meshDimension)
1555 TiXmlElement *exp =
new TiXmlElement(
"E");
1558 "C[" + boost::lexical_cast<string>(it->first) +
"]");
1559 exp->SetAttribute(
"NUMMODES", 4);
1560 exp->SetAttribute(
"TYPE",
"MODIFIED");
1561 exp->SetAttribute(
"FIELDS",
"u");
1563 expTag->LinkEndChild(exp);
1566 root->LinkEndChild(expTag);
1569 doc->SaveFile(filenameXml);
1576 m_file = H5::File::Create(filenameHdf5, H5F_ACC_TRUNC);
1577 auto hdfRoot = m_file->CreateGroup(
"NEKTAR");
1578 auto hdfRoot2 = hdfRoot->CreateGroup(
"GEOMETRY");
1581 hdfRoot2->SetAttribute(
"FORMAT_VERSION", FORMAT_VERSION);
1584 m_mesh = hdfRoot2->CreateGroup(
"MESH");
1585 m_maps = hdfRoot2->CreateGroup(
"MAPS");
1587 WriteGeometryMap(m_vertSet,
"VERT");
1588 WriteGeometryMap(m_segGeoms,
"SEG");
1589 if (m_meshDimension > 1)
1591 WriteGeometryMap(m_triGeoms,
"TRI");
1592 WriteGeometryMap(m_quadGeoms,
"QUAD");
1594 if (m_meshDimension > 2)
1596 WriteGeometryMap(m_tetGeoms,
"TET");
1597 WriteGeometryMap(m_pyrGeoms,
"PYR");
1598 WriteGeometryMap(m_prismGeoms,
"PRISM");
1599 WriteGeometryMap(m_hexGeoms,
"HEX");
1603 int ptOffset = 0, newIdx = 0;
1605 WriteCurveMap(m_curvedEdges,
"CURVE_EDGE", curvePts, ptOffset, newIdx);
1606 WriteCurveMap(m_curvedFaces,
"CURVE_FACE", curvePts, ptOffset, newIdx);
1607 WriteCurvePoints(curvePts);
1610 WriteComposites(m_meshComposites);
1611 WriteDomain(m_domain);
#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.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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::shared_ptr< QuadGeom > QuadGeomSharedPtr
std::map< int, std::pair< LibUtilities::ShapeType, std::vector< int > > > CompositeDescriptor
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
std::unordered_map< int, CurveSharedPtr > CurveMap
std::map< int, CompositeSharedPtr > CompositeMap
int GetGeomData(std::shared_ptr< T > &geom, int i)
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< MeshPartition > MeshPartitionSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
MeshGraphFactory & GetMeshGraphFactory()
std::shared_ptr< Composite > CompositeSharedPtr
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.
std::vector< MeshVertex > pts
mapping to access pts value.
std::vector< NekInt64 > index
id of this Point set
std::vector< unsigned int > list