47#include <boost/format.hpp>
64 const bool isRoot = comm->TreatAsRankZero();
78 int isPartitioned = 0;
81 if (
m_session->DefinesElement(
"Nektar/Geometry"))
83 if (
m_session->GetElement(
"Nektar/Geometry")
84 ->Attribute(
"PARTITION"))
86 std::cout <<
"Using pre-partitioned mesh." << std::endl;
90 ->QueryIntAttribute(
"DIM", &meshDimension);
93 comm->Bcast(isPartitioned, 0);
94 comm->Bcast(meshDimension, 0);
111 std::string partitionerName =
"Scotch";
112 if (!haveScotch && haveMetis)
114 partitionerName =
"Metis";
118 if (session->DefinesCmdLineArgument(
"use-metis"))
120 partitionerName =
"Metis";
122 if (session->DefinesCmdLineArgument(
"use-scotch"))
124 partitionerName =
"Scotch";
130 if (session->DefinesCmdLineArgument(
"part-only") ||
131 session->DefinesCmdLineArgument(
"part-only-overlapping"))
136 "The 'part-only' option should be used in serial.");
147 partitionerName, session, session->GetComm(), meshDimension,
150 if (session->DefinesCmdLineArgument(
"part-only"))
152 nParts = session->GetCmdLineArgument<
int>(
"part-only");
153 partitioner->PartitionMesh(nParts,
true);
158 session->GetCmdLineArgument<
int>(
"part-only-overlapping");
159 partitioner->PartitionMesh(nParts,
true,
true);
162 std::vector<std::set<unsigned int>> elmtIDs;
163 std::vector<unsigned int> parts(nParts);
164 for (
int i = 0; i < nParts; ++i)
166 std::vector<unsigned int> elIDs;
167 std::set<unsigned int> tmp;
168 partitioner->GetElementIDs(i, elIDs);
169 tmp.insert(elIDs.begin(), elIDs.end());
170 elmtIDs.push_back(tmp);
176 if (isRoot && session->DefinesCmdLineArgument(
"part-info"))
178 partitioner->PrintPartInfo(std::cout);
185 if (commMesh->GetSize() > 1)
188 "Valid partitioner not found! Either Scotch or METIS "
191 int nParts = commMesh->GetSize();
193 if (session->GetSharedFilesystem())
195 std::vector<unsigned int> keys, vals;
211 partitionerName, session, session->GetComm(),
214 partitioner->PartitionMesh(nParts,
true);
216 std::vector<std::set<unsigned int>> elmtIDs;
217 std::vector<unsigned int> parts(nParts);
218 for (i = 0; i < nParts; ++i)
220 std::vector<unsigned int> elIDs;
221 std::set<unsigned int> tmp;
222 partitioner->GetElementIDs(i, elIDs);
223 tmp.insert(elIDs.begin(), elIDs.end());
224 elmtIDs.push_back(tmp);
240 comm->Bcast(keys, 0);
251 vals[i++] = cIt.second.size();
255 comm->Bcast(keys, 0);
256 comm->Bcast(vals, 0);
259 comm->Bcast(cIt.second, 0);
272 vals[i++] = bIt.second.size();
278 comm->Bcast(keys, 0);
282 comm->Bcast(vals, 0);
286 comm->Bcast(bIt.second, 0);
290 if (session->DefinesCmdLineArgument(
"part-info"))
292 partitioner->PrintPartInfo(std::cout);
298 comm->Bcast(keys, 0);
300 int cmpSize = keys[0];
301 int bndSize = keys[1];
303 keys.resize(cmpSize);
304 vals.resize(cmpSize);
305 comm->Bcast(keys, 0);
306 comm->Bcast(vals, 0);
308 for (
int i = 0; i < keys.size(); ++i)
310 std::vector<unsigned int> tmp(vals[i]);
316 keys.resize(bndSize);
317 vals.resize(bndSize);
320 comm->Bcast(keys, 0);
324 comm->Bcast(vals, 0);
326 for (
int i = 0; i < keys.size(); ++i)
328 std::vector<unsigned int> tmp(vals[i]);
349 partitionerName, session, session->GetComm(),
352 partitioner->PartitionMesh(nParts,
false);
354 std::vector<unsigned int> parts(1), tmp;
355 parts[0] = commMesh->GetRank();
356 std::vector<std::set<unsigned int>> elIDs(1);
357 partitioner->GetElementIDs(parts[0], tmp);
358 elIDs[0].insert(tmp.begin(), tmp.end());
367 if (
m_session->DefinesCmdLineArgument(
"part-info") && isRoot)
369 partitioner->PrintPartInfo(std::cout);
376 std::string dirname =
m_session->GetSessionName() +
"_xml";
377 fs::path pdirname(dirname);
379 pad % comm->GetRowComm()->GetRank();
380 fs::path pFilename(pad.str());
381 fs::path fullpath = pdirname / pFilename;
383 std::vector<std::string> filenames = {
408 TiXmlAttribute *attr =
m_xmlGeom->FirstAttribute();
413 int meshDimension =
m_meshGraph->GetMeshDimension();
414 int spaceDimension =
m_meshGraph->GetSpaceDimension();
421 std::string attrName(attr->Name());
422 if (attrName ==
"DIM")
424 err = attr->QueryIntValue(&meshDimension);
425 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read mesh dimension.");
427 else if (attrName ==
"SPACE")
429 err = attr->QueryIntValue(&spaceDimension);
430 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read space dimension.");
432 else if (attrName ==
"PARTITION")
435 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read partition.");
441 std::string errstr(
"Unknown attribute: ");
453 ASSERTL0(meshDimension <= spaceDimension,
454 "Mesh dimension greater than space dimension");
458 if (meshDimension >= 2)
461 if (meshDimension == 3)
479 int spaceDimension =
m_meshGraph->GetSpaceDimension();
482 TiXmlElement *element =
m_xmlGeom->FirstChildElement(
"VERTEX");
483 ASSERTL0(element,
"Unable to find mesh VERTEX tag in file.");
490 const char *xscal = element->Attribute(
"XSCALE");
497 std::string xscalstr = xscal;
499 xscale = expEvaluator.
Evaluate(expr_id);
502 const char *yscal = element->Attribute(
"YSCALE");
509 std::string yscalstr = yscal;
511 yscale = expEvaluator.
Evaluate(expr_id);
514 const char *zscal = element->Attribute(
"ZSCALE");
521 std::string zscalstr = zscal;
523 zscale = expEvaluator.
Evaluate(expr_id);
531 const char *xmov = element->Attribute(
"XMOVE");
538 std::string xmovstr = xmov;
540 xmove = expEvaluator.
Evaluate(expr_id);
543 const char *ymov = element->Attribute(
"YMOVE");
550 std::string ymovstr = ymov;
552 ymove = expEvaluator.
Evaluate(expr_id);
555 const char *zmov = element->Attribute(
"ZMOVE");
562 std::string zmovstr = zmov;
564 zmove = expEvaluator.
Evaluate(expr_id);
567 TiXmlElement *vertex = element->FirstChildElement(
"V");
573 TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
574 std::string attrName(vertexAttr->Name());
577 (std::string(
"Unknown attribute name: ") + attrName).c_str());
579 int err = vertexAttr->QueryIntValue(&indx);
580 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
583 std::string vertexBodyStr;
585 TiXmlNode *vertexBody = vertex->FirstChild();
590 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
592 vertexBodyStr += vertexBody->ToText()->Value();
593 vertexBodyStr +=
" ";
596 vertexBody = vertexBody->NextSibling();
600 "Vertex definitions must contain vertex data.");
604 std::istringstream vertexDataStrm(vertexBodyStr.c_str());
608 while (!vertexDataStrm.fail())
610 vertexDataStrm >> xval >> yval >> zval;
612 xval = xval * xscale + xmove;
613 yval = yval * yscale + ymove;
614 zval = zval * zscale + zmove;
619 if (!vertexDataStrm.fail())
623 spaceDimension, indx, xval, yval, zval));
624 vertSet[indx] = vert;
630 ASSERTL0(
false,
"Unable to read VERTEX data.");
633 vertex = vertex->NextSiblingElement(
"V");
641 int meshDimension =
m_meshGraph->GetMeshDimension();
645 TiXmlElement *element =
m_xmlGeom->FirstChildElement(
"VERTEX");
646 ASSERTL0(element,
"Unable to find mesh VERTEX tag in file.");
651 const char *xscal = element->Attribute(
"XSCALE");
658 std::string xscalstr = xscal;
660 xscale = expEvaluator.
Evaluate(expr_id);
663 const char *yscal = element->Attribute(
"YSCALE");
670 std::string yscalstr = yscal;
672 yscale = expEvaluator.
Evaluate(expr_id);
675 const char *zscal = element->Attribute(
"ZSCALE");
682 std::string zscalstr = zscal;
684 zscale = expEvaluator.
Evaluate(expr_id);
692 const char *xmov = element->Attribute(
"XMOVE");
699 std::string xmovstr = xmov;
701 xmove = expEvaluator.
Evaluate(expr_id);
704 const char *ymov = element->Attribute(
"YMOVE");
711 std::string ymovstr = ymov;
713 ymove = expEvaluator.
Evaluate(expr_id);
716 const char *zmov = element->Attribute(
"ZMOVE");
723 std::string zmovstr = zmov;
725 zmove = expEvaluator.
Evaluate(expr_id);
742 TiXmlElement *edgelement =
field->FirstChildElement(
"E");
744 int edgeindx, edgeid;
748 std::string edge(edgelement->ValueStr());
750 (std::string(
"Unknown 3D curve type:") + edge).c_str());
753 err = edgelement->QueryIntAttribute(
"ID", &edgeindx);
754 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read curve attribute ID.");
757 err = edgelement->QueryIntAttribute(
"EDGEID", &edgeid);
759 "Unable to read curve attribute EDGEID.");
762 std::string elementStr;
763 TiXmlNode *elementChild = edgelement->FirstChild();
768 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
770 elementStr += elementChild->ToText()->ValueStr();
773 elementChild = elementChild->NextSibling();
776 ASSERTL0(!elementStr.empty(),
"Unable to read curve description body.");
784 std::string typeStr = edgelement->Attribute(
"TYPE");
785 ASSERTL0(!typeStr.empty(),
"TYPE must be specified in "
786 "points definition");
790 const std::string *endStr =
792 const std::string *ptsStr =
std::find(begStr, endStr, typeStr);
794 ASSERTL0(ptsStr != endStr,
"Invalid points type.");
798 err = edgelement->QueryIntAttribute(
"NUMPOINTS", &numPts);
800 "Unable to read curve attribute NUMPOINTS.");
806 std::istringstream elementDataStrm(elementStr.c_str());
809 while (!elementDataStrm.fail())
811 elementDataStrm >> xval >> yval >> zval;
813 xval = xval * xscale + xmove;
814 yval = yval * yscale + ymove;
815 zval = zval * zscale + zmove;
820 if (!elementDataStrm.fail())
824 meshDimension, edgeindx, xval, yval, zval));
826 curve->m_points.push_back(vert);
833 (std::string(
"Unable to read curve data for EDGE: ") +
838 ASSERTL0(curve->m_points.size() == numPts,
839 "Number of points specificed by attribute "
840 "NUMPOINTS is different from number of points "
841 "in list (edgeid = " +
842 std::to_string(edgeid));
844 curvedEdges[edgeid] = curve;
846 edgelement = edgelement->NextSiblingElement(
"E");
852 TiXmlElement *facelement =
field->FirstChildElement(
"F");
853 int faceindx, faceid;
857 std::string face(facelement->ValueStr());
859 (std::string(
"Unknown 3D curve type: ") + face).c_str());
862 err = facelement->QueryIntAttribute(
"ID", &faceindx);
863 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read curve attribute ID.");
866 err = facelement->QueryIntAttribute(
"FACEID", &faceid);
868 "Unable to read curve attribute FACEID.");
871 std::string elementStr;
872 TiXmlNode *elementChild = facelement->FirstChild();
877 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
879 elementStr += elementChild->ToText()->ValueStr();
882 elementChild = elementChild->NextSibling();
885 ASSERTL0(!elementStr.empty(),
"Unable to read curve description body.");
891 std::string typeStr = facelement->Attribute(
"TYPE");
892 ASSERTL0(!typeStr.empty(),
"TYPE must be specified in "
893 "points definition");
896 const std::string *endStr =
898 const std::string *ptsStr =
std::find(begStr, endStr, typeStr);
900 ASSERTL0(ptsStr != endStr,
"Invalid points type.");
903 std::string numptsStr = facelement->Attribute(
"NUMPOINTS");
905 "NUMPOINTS must be specified in points definition");
914 ASSERTL0(numPts >= 3,
"NUMPOINTS for face must be greater than 2");
918 ASSERTL0(ptsStr != endStr,
"Invalid points type.");
923 std::istringstream elementDataStrm(elementStr.c_str());
926 while (!elementDataStrm.fail())
928 elementDataStrm >> xval >> yval >> zval;
934 if (!elementDataStrm.fail())
938 meshDimension, faceindx, xval, yval, zval));
939 curve->m_points.push_back(vert);
946 (std::string(
"Unable to read curve data for FACE: ") +
950 curvedFaces[faceid] = curve;
952 facelement = facelement->NextSiblingElement(
"F");
961 TiXmlElement *element =
nullptr;
963 element =
m_xmlGeom->FirstChildElement(
"DOMAIN");
965 ASSERTL0(element,
"Unable to find DOMAIN tag in file.");
969 TiXmlElement *multidomains = element->FirstChildElement(
"D");
976 int err = multidomains->QueryIntAttribute(
"ID", &indx);
978 "Unable to read attribute ID in Domain.");
980 TiXmlNode *elementChild = multidomains->FirstChild();
981 while (elementChild &&
982 elementChild->Type() != TiXmlNode::TINYXML_TEXT)
984 elementChild = elementChild->NextSibling();
987 ASSERTL0(elementChild,
"Unable to read DOMAIN body.");
988 std::string elementStr = elementChild->ToText()->ValueStr();
990 elementStr = elementStr.substr(elementStr.find_first_not_of(
" "));
992 std::string::size_type indxBeg = elementStr.find_first_of(
'[') + 1;
993 std::string::size_type indxEnd = elementStr.find_last_of(
']') - 1;
994 std::string indxStr =
995 elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
999 "Unable to read domain's composite index (index missing?).");
1003 std::map<int, CompositeSharedPtr> unrollDomain;
1004 m_meshGraph->GetCompositeList(indxStr, unrollDomain);
1005 domain[indx] = unrollDomain;
1009 "Unable to obtain domain's referenced composite: ") +
1014 multidomains = multidomains->NextSiblingElement(
"D");
1021 TiXmlNode *elementChild = element->FirstChild();
1022 while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1024 elementChild = elementChild->NextSibling();
1027 ASSERTL0(elementChild,
"Unable to read DOMAIN body.");
1028 std::string elementStr = elementChild->ToText()->ValueStr();
1030 elementStr = elementStr.substr(elementStr.find_first_not_of(
" "));
1032 std::string::size_type indxBeg = elementStr.find_first_of(
'[') + 1;
1033 std::string::size_type indxEnd = elementStr.find_last_of(
']') - 1;
1034 std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
1037 "Unable to read domain's composite index (index missing?).");
1041 std::map<int, CompositeSharedPtr> fullDomain;
1042 m_meshGraph->GetCompositeList(indxStr, fullDomain);
1043 domain[0] = fullDomain;
1047 (std::string(
"Unable to obtain domain's referenced composite: ") +
1056 auto &curvedEdges =
m_meshGraph->GetCurvedEdges();
1057 int spaceDimension =
m_meshGraph->GetSpaceDimension();
1059 CurveMap::iterator it;
1069 TiXmlElement *edgeElement =
field->FirstChildElement(
"E");
1077 std::string edgeStr;
1082 int err = edgeElement->QueryIntAttribute(
"ID", &indx);
1083 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read edge attribute ID.");
1085 TiXmlNode *child = edgeElement->FirstChild();
1087 if (child->Type() == TiXmlNode::TINYXML_TEXT)
1089 edgeStr += child->ToText()->ValueStr();
1093 int vertex1, vertex2;
1094 std::istringstream edgeDataStrm(edgeStr.c_str());
1098 while (!edgeDataStrm.fail())
1100 edgeDataStrm >> vertex1 >> vertex2;
1107 if (!edgeDataStrm.fail())
1112 it = curvedEdges.find(indx);
1114 if (it == curvedEdges.end())
1118 indx, spaceDimension, vertices);
1124 indx, spaceDimension, vertices, it->second);
1133 (std::string(
"Unable to read edge data: ") + edgeStr).c_str());
1136 edgeElement = edgeElement->NextSiblingElement(
"E");
1142 auto &curvedFaces =
m_meshGraph->GetCurvedFaces();
1155 TiXmlElement *element =
field->FirstChildElement();
1156 CurveMap::iterator it;
1160 std::string elementType(element->ValueStr());
1162 ASSERTL0(elementType ==
"Q" || elementType ==
"T",
1163 (std::string(
"Unknown 3D face type: ") + elementType).c_str());
1167 int err = element->QueryIntAttribute(
"ID", &indx);
1168 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read face attribute ID.");
1171 it = curvedFaces.find(indx);
1174 TiXmlNode *elementChild = element->FirstChild();
1175 std::string elementStr;
1176 while (elementChild)
1178 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1180 elementStr += elementChild->ToText()->ValueStr();
1182 elementChild = elementChild->NextSibling();
1185 ASSERTL0(!elementStr.empty(),
"Unable to read face description body.");
1189 if (elementType ==
"T")
1192 int edge1, edge2, edge3;
1193 std::istringstream elementDataStrm(elementStr.c_str());
1197 elementDataStrm >> edge1;
1198 elementDataStrm >> edge2;
1199 elementDataStrm >> edge3;
1202 !elementDataStrm.fail(),
1203 (std::string(
"Unable to read face data for TRIANGLE: ") +
1213 if (it == curvedFaces.end())
1221 indx, edges, it->second);
1223 triGeoms[indx]->SetGlobalID(indx);
1229 (std::string(
"Unable to read face data for TRIANGLE: ") +
1234 else if (elementType ==
"Q")
1237 int edge1, edge2, edge3, edge4;
1238 std::istringstream elementDataStrm(elementStr.c_str());
1242 elementDataStrm >> edge1;
1243 elementDataStrm >> edge2;
1244 elementDataStrm >> edge3;
1245 elementDataStrm >> edge4;
1248 (std::string(
"Unable to read face data for QUAD: ") +
1261 if (it == curvedFaces.end())
1272 quadGeoms[indx]->SetGlobalID(indx);
1277 (std::string(
"Unable to read face data for QUAD: ") +
1283 element = element->NextSiblingElement();
1289 int meshDimension =
m_meshGraph->GetMeshDimension();
1291 switch (meshDimension)
1307 auto &curvedEdges =
m_meshGraph->GetCurvedEdges();
1309 int spaceDimension =
m_meshGraph->GetSpaceDimension();
1311 TiXmlElement *
field =
nullptr;
1321 TiXmlElement *segment =
field->FirstChildElement(
"S");
1322 CurveMap::iterator it;
1327 int err = segment->QueryIntAttribute(
"ID", &indx);
1328 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read element attribute ID.");
1330 TiXmlNode *elementChild = segment->FirstChild();
1331 while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1333 elementChild = elementChild->NextSibling();
1336 ASSERTL0(elementChild,
"Unable to read element description body.");
1337 std::string elementStr = elementChild->ToText()->ValueStr();
1342 int vertex1, vertex2;
1343 std::istringstream elementDataStrm(elementStr.c_str());
1347 elementDataStrm >> vertex1;
1348 elementDataStrm >> vertex2;
1351 (std::string(
"Unable to read element data for SEGMENT: ") +
1357 it = curvedEdges.find(indx);
1359 if (it == curvedEdges.end())
1362 indx, spaceDimension, vertices);
1367 indx, spaceDimension, vertices, it->second);
1369 segGeoms[indx]->SetGlobalID(indx);
1374 (std::string(
"Unable to read element data for segment: ") +
1379 segment = segment->NextSiblingElement(
"S");
1385 auto &curvedFaces =
m_meshGraph->GetCurvedFaces();
1395 CurveMap::iterator it;
1400 TiXmlElement *element =
field->FirstChildElement();
1404 std::string elementType(element->ValueStr());
1407 elementType ==
"Q" || elementType ==
"T",
1408 (std::string(
"Unknown 2D element type: ") + elementType).c_str());
1412 int err = element->QueryIntAttribute(
"ID", &indx);
1413 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read element attribute ID.");
1415 it = curvedFaces.find(indx);
1418 TiXmlNode *elementChild = element->FirstChild();
1419 std::string elementStr;
1420 while (elementChild)
1422 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1424 elementStr += elementChild->ToText()->ValueStr();
1426 elementChild = elementChild->NextSibling();
1430 "Unable to read element description body.");
1434 if (elementType ==
"T")
1437 int edge1, edge2, edge3;
1438 std::istringstream elementDataStrm(elementStr.c_str());
1442 elementDataStrm >> edge1;
1443 elementDataStrm >> edge2;
1444 elementDataStrm >> edge3;
1447 !elementDataStrm.fail(),
1448 (std::string(
"Unable to read element data for TRIANGLE: ") +
1458 if (it == curvedFaces.end())
1466 indx, edges, it->second);
1468 triGeoms[indx]->SetGlobalID(indx);
1474 (std::string(
"Unable to read element data for TRIANGLE: ") +
1479 else if (elementType ==
"Q")
1482 int edge1, edge2, edge3, edge4;
1483 std::istringstream elementDataStrm(elementStr.c_str());
1487 elementDataStrm >> edge1;
1488 elementDataStrm >> edge2;
1489 elementDataStrm >> edge3;
1490 elementDataStrm >> edge4;
1493 !elementDataStrm.fail(),
1494 (std::string(
"Unable to read element data for QUAD: ") +
1506 if (it == curvedFaces.end())
1517 quadGeoms[indx]->SetGlobalID(indx);
1523 (std::string(
"Unable to read element data for QUAD: ") +
1529 element = element->NextSiblingElement();
1537 auto &prismGeoms =
m_meshGraph->GetAllPrismGeoms();
1548 TiXmlElement *element =
field->FirstChildElement();
1552 std::string elementType(element->ValueStr());
1556 elementType ==
"A" || elementType ==
"P" || elementType ==
"R" ||
1558 (std::string(
"Unknown 3D element type: ") + elementType).c_str());
1562 int err = element->QueryIntAttribute(
"ID", &indx);
1563 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read element attribute ID.");
1566 TiXmlNode *elementChild = element->FirstChild();
1567 std::string elementStr;
1568 while (elementChild)
1570 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1572 elementStr += elementChild->ToText()->ValueStr();
1574 elementChild = elementChild->NextSibling();
1578 "Unable to read element description body.");
1580 std::istringstream elementDataStrm(elementStr.c_str());
1586 if (elementType ==
"A")
1600 std::stringstream errorstring;
1601 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1602 <<
" triangle face(s), and " << kNqfaces
1603 <<
" quadrilateral face(s).";
1604 for (
int i = 0; i < kNfaces; i++)
1607 elementDataStrm >> faceID;
1614 std::stringstream errorstring;
1615 errorstring <<
"Element " << indx
1616 <<
" has invalid face: " << faceID;
1617 ASSERTL0(
false, errorstring.str().c_str());
1621 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1623 std::static_pointer_cast<TriGeom>(face);
1625 else if (face->GetShapeType() ==
1628 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1636 "Unable to read element data for TETRAHEDRON: ") +
1639 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1640 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1645 tetGeoms[indx] = tetgeom;
1646 m_meshGraph->PopulateFaceToElMap(tetgeom, kNfaces);
1652 "Unable to read element data for TETRAHEDRON: ") +
1658 else if (elementType ==
"P")
1673 std::stringstream errorstring;
1674 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1675 <<
" triangle face(s), and " << kNqfaces
1676 <<
" quadrilateral face(s).";
1677 for (
int i = 0; i < kNfaces; i++)
1680 elementDataStrm >> faceID;
1687 std::stringstream errorstring;
1688 errorstring <<
"Element " << indx
1689 <<
" has invalid face: " << faceID;
1690 ASSERTL0(
false, errorstring.str().c_str());
1694 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1696 std::static_pointer_cast<TriGeom>(face);
1699 else if (face->GetShapeType() ==
1702 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1704 std::static_pointer_cast<QuadGeom>(face);
1712 !elementDataStrm.fail(),
1713 (std::string(
"Unable to read element data for PYRAMID: ") +
1716 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1717 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1722 pyrGeoms[indx] = pyrgeom;
1723 m_meshGraph->PopulateFaceToElMap(pyrgeom, kNfaces);
1729 (std::string(
"Unable to read element data for PYRAMID: ") +
1735 else if (elementType ==
"R")
1750 std::stringstream errorstring;
1751 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1752 <<
" triangle face(s), and " << kNqfaces
1753 <<
" quadrilateral face(s).";
1755 for (
int i = 0; i < kNfaces; i++)
1758 elementDataStrm >> faceID;
1765 std::stringstream errorstring;
1766 errorstring <<
"Element " << indx
1767 <<
" has invalid face: " << faceID;
1768 ASSERTL0(
false, errorstring.str().c_str());
1772 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1774 std::static_pointer_cast<TriGeom>(face);
1777 else if (face->GetShapeType() ==
1780 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1782 std::static_pointer_cast<QuadGeom>(face);
1790 !elementDataStrm.fail(),
1791 (std::string(
"Unable to read element data for PRISM: ") +
1794 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1795 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1800 prismGeoms[indx] = prismgeom;
1801 m_meshGraph->PopulateFaceToElMap(prismgeom, kNfaces);
1807 (std::string(
"Unable to read element data for PRISM: ") +
1813 else if (elementType ==
"H")
1828 std::stringstream errorstring;
1829 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1830 <<
" triangle face(s), and " << kNqfaces
1831 <<
" quadrilateral face(s).";
1832 for (
int i = 0; i < kNfaces; i++)
1835 elementDataStrm >> faceID;
1842 std::stringstream errorstring;
1843 errorstring <<
"Element " << indx
1844 <<
" has invalid face: " << faceID;
1845 ASSERTL0(
false, errorstring.str().c_str());
1849 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1851 else if (face->GetShapeType() ==
1854 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1856 std::static_pointer_cast<QuadGeom>(face);
1864 "Unable to read element data for HEXAHEDRAL: ") +
1867 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1868 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1873 hexGeoms[indx] = hexgeom;
1874 m_meshGraph->PopulateFaceToElMap(hexgeom, kNfaces);
1880 "Unable to read element data for HEXAHEDRAL: ") +
1886 element = element->NextSiblingElement();
1892 auto &meshComposites =
m_meshGraph->GetComposites();
1893 auto &compositesLabels =
m_meshGraph->GetCompositesLabels();
1895 TiXmlElement *
field =
nullptr;
1902 TiXmlElement *node =
field->FirstChildElement(
"C");
1905 int nextCompositeNumber = -1;
1912 nextCompositeNumber++;
1915 int err = node->QueryIntAttribute(
"ID", &indx);
1916 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
1920 TiXmlNode *compositeChild = node->FirstChild();
1925 while (compositeChild &&
1926 compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
1928 compositeChild = compositeChild->NextSibling();
1931 ASSERTL0(compositeChild,
"Unable to read composite definition body.");
1932 std::string compositeStr = compositeChild->ToText()->ValueStr();
1936 std::istringstream compositeDataStrm(compositeStr.c_str());
1941 std::string prevCompositeElementStr;
1943 while (!compositeDataStrm.fail())
1945 std::string compositeElementStr;
1946 compositeDataStrm >> compositeElementStr;
1948 if (!compositeDataStrm.fail())
1953 meshComposites[indx] =
1957 if (compositeElementStr.length() > 0)
1960 compositeElementStr,
1961 meshComposites[indx]);
1963 prevCompositeElementStr = compositeElementStr;
1971 (std::string(
"Unable to read COMPOSITE data for composite: ") +
1978 err = node->QueryStringAttribute(
"NAME", &
name);
1979 if (err == TIXML_SUCCESS)
1981 compositesLabels[indx] =
name;
1985 node = node->NextSiblingElement(
"C");
1989 "At least one composite must be specified.");
1993 const std::string &token,
1996 int meshDimension =
m_meshGraph->GetMeshDimension();
1998 switch (meshDimension)
2013 const std::string &token,
2021 std::istringstream tokenStream(token);
2022 std::istringstream prevTokenStream(prevToken);
2027 tokenStream >> type;
2029 std::string::size_type indxBeg = token.find_first_of(
'[') + 1;
2030 std::string::size_type indxEnd = token.find_last_of(
']') - 1;
2034 (std::string(
"Error reading index definition:") + token).c_str());
2036 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2038 typedef std::vector<unsigned int> SeqVectorType;
2039 SeqVectorType seqVector;
2044 (std::string(
"Ill-formed sequence definition: ") + indxStr)
2048 prevTokenStream >> prevType;
2051 bool validSequence =
2052 (prevToken.empty() ||
2053 (type ==
'V' && prevType ==
'V') ||
2054 (type ==
'S' && prevType ==
'S'));
2057 std::string(
"Invalid combination of composite items: ") +
2058 type +
" and " + prevType +
".");
2063 for (SeqVectorType::iterator iter = seqVector.begin();
2064 iter != seqVector.end(); ++iter)
2066 if (vertSet.find(*iter) == vertSet.end())
2069 "Unknown vertex index: " +
2070 std::to_string(*iter));
2074 composite->m_geomVec.push_back(vertSet[*iter]);
2080 for (SeqVectorType::iterator iter = seqVector.begin();
2081 iter != seqVector.end(); ++iter)
2083 if (segGeoms.find(*iter) == segGeoms.end())
2086 "Unknown segment index: " +
2087 std::to_string(*iter));
2091 composite->m_geomVec.push_back(segGeoms[*iter]);
2098 "Unrecognized composite token: " + token);
2104 "Problem processing composite token: " + token);
2111 const std::string &token,
2121 std::istringstream tokenStream(token);
2122 std::istringstream prevTokenStream(prevToken);
2127 tokenStream >> type;
2129 std::string::size_type indxBeg = token.find_first_of(
'[') + 1;
2130 std::string::size_type indxEnd = token.find_last_of(
']') - 1;
2134 (std::string(
"Error reading index definition:") + token).c_str());
2136 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2137 std::vector<unsigned int> seqVector;
2138 std::vector<unsigned int>::iterator seqIter;
2143 (std::string(
"Error reading composite elements: ") + indxStr)
2146 prevTokenStream >> prevType;
2149 bool validSequence =
2150 (prevToken.empty() ||
2151 (type ==
'V' && prevType ==
'V') ||
2152 (type ==
'E' && prevType ==
'E') ||
2153 ((type ==
'T' || type ==
'Q') &&
2154 (prevType ==
'T' || prevType ==
'Q')));
2157 std::string(
"Invalid combination of composite items: ") +
2158 type +
" and " + prevType +
".");
2163 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2166 if (segGeoms.find(*seqIter) == segGeoms.end())
2169 "Unknown edge index: " +
2170 std::to_string(*seqIter));
2174 composite->m_geomVec.push_back(segGeoms[*seqIter]);
2181 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2184 if (triGeoms.count(*seqIter) == 0)
2187 "Unknown triangle index: " +
2188 std::to_string(*seqIter));
2194 composite->m_geomVec.push_back(triGeoms[*seqIter]);
2203 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2206 if (quadGeoms.count(*seqIter) == 0)
2210 "Unknown quad index: " + std::to_string(*seqIter) +
2211 " in Composite section");
2215 if (
m_meshGraph->CheckRange(*quadGeoms[*seqIter]))
2217 composite->m_geomVec.push_back(quadGeoms[*seqIter]);
2225 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2228 if (*seqIter >= vertSet.size())
2231 "Unknown vertex index: " +
2232 std::to_string(*seqIter));
2236 composite->m_geomVec.push_back(vertSet[*seqIter]);
2243 "Unrecognized composite token: " + token);
2249 "Problem processing composite token: " + token);
2256 const std::string &token,
2270 std::istringstream tokenStream(token);
2271 std::istringstream prevTokenStream(prevToken);
2276 tokenStream >> type;
2278 std::string::size_type indxBeg = token.find_first_of(
'[') + 1;
2279 std::string::size_type indxEnd = token.find_last_of(
']') - 1;
2282 "Error reading index definition: " + token);
2284 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2286 std::vector<unsigned int> seqVector;
2287 std::vector<unsigned int>::iterator seqIter;
2291 ASSERTL0(err,
"Error reading composite elements: " + indxStr);
2293 prevTokenStream >> prevType;
2297 std::map<char, int> typeMap;
2308 bool validSequence =
2309 (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
2312 std::string(
"Invalid combination of composite items: ") +
2313 type +
" and " + prevType +
".");
2318 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2321 if (vertSet.find(*seqIter) == vertSet.end())
2324 "Unknown vertex index: " +
2325 std::to_string(*seqIter));
2329 composite->m_geomVec.push_back(vertSet[*seqIter]);
2335 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2338 if (segGeoms.find(*seqIter) == segGeoms.end())
2341 "Unknown edge index: " +
2342 std::to_string(*seqIter));
2346 composite->m_geomVec.push_back(segGeoms[*seqIter]);
2352 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2360 "Unknown face index: " +
2361 std::to_string(*seqIter));
2367 composite->m_geomVec.push_back(face);
2374 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2377 if (triGeoms.find(*seqIter) == triGeoms.end())
2380 "Unknown triangle index: " +
2381 std::to_string(*seqIter));
2387 composite->m_geomVec.push_back(triGeoms[*seqIter]);
2394 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2397 if (quadGeoms.find(*seqIter) == quadGeoms.end())
2400 "Unknown quad index: " +
2401 std::to_string(*seqIter));
2405 if (
m_meshGraph->CheckRange(*quadGeoms[*seqIter]))
2407 composite->m_geomVec.push_back(quadGeoms[*seqIter]);
2415 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2418 if (tetGeoms.find(*seqIter) == tetGeoms.end())
2421 "Unknown tet index: " +
2422 std::to_string(*seqIter));
2428 composite->m_geomVec.push_back(tetGeoms[*seqIter]);
2436 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2439 if (pyrGeoms.find(*seqIter) == pyrGeoms.end())
2442 "Unknown pyramid index: " +
2443 std::to_string(*seqIter));
2449 composite->m_geomVec.push_back(pyrGeoms[*seqIter]);
2457 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2460 if (prismGeoms.find(*seqIter) == prismGeoms.end())
2463 "Unknown prism index: " +
2464 std::to_string(*seqIter));
2468 if (
m_meshGraph->CheckRange(*prismGeoms[*seqIter]))
2470 composite->m_geomVec.push_back(
2471 prismGeoms[*seqIter]);
2479 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2482 if (hexGeoms.find(*seqIter) == hexGeoms.end())
2485 "Unknown hex index: " +
2486 std::to_string(*seqIter));
2492 composite->m_geomVec.push_back(hexGeoms[*seqIter]);
2500 "Unrecognized composite token: " + token);
2506 "Problem processing composite token: " + token);
2514 TiXmlElement *vertTag =
new TiXmlElement(
"VERTEX");
2516 for (
auto &i : verts)
2518 std::stringstream s;
2519 s << std::scientific << std::setprecision(8) << (*i.second)(0) <<
" "
2520 << (*i.second)(1) <<
" " << (*i.second)(2);
2521 TiXmlElement *v =
new TiXmlElement(
"V");
2522 v->SetAttribute(
"ID", i.second->GetGlobalID());
2523 v->LinkEndChild(
new TiXmlText(s.str()));
2524 vertTag->LinkEndChild(v);
2527 geomTag->LinkEndChild(vertTag);
2532 int meshDimension =
m_meshGraph->GetMeshDimension();
2534 TiXmlElement *edgeTag =
2535 new TiXmlElement(meshDimension == 1 ?
"ELEMENT" :
"EDGE");
2536 std::string tag = meshDimension == 1 ?
"S" :
"E";
2538 for (
auto &i : edges)
2540 std::stringstream s;
2542 s << seg->GetVid(0) <<
" " << seg->GetVid(1);
2543 TiXmlElement *e =
new TiXmlElement(tag);
2544 e->SetAttribute(
"ID", i.first);
2545 e->LinkEndChild(
new TiXmlText(s.str()));
2546 edgeTag->LinkEndChild(e);
2549 geomTag->LinkEndChild(edgeTag);
2554 std::string tag =
"T";
2556 for (
auto &i : tris)
2558 std::stringstream s;
2560 s << tri->GetEid(0) <<
" " << tri->GetEid(1) <<
" " << tri->GetEid(2);
2561 TiXmlElement *t =
new TiXmlElement(tag);
2562 t->SetAttribute(
"ID", i.first);
2563 t->LinkEndChild(
new TiXmlText(s.str()));
2564 faceTag->LinkEndChild(t);
2570 std::string tag =
"Q";
2572 for (
auto &i : quads)
2574 std::stringstream s;
2576 s << quad->GetEid(0) <<
" " << quad->GetEid(1) <<
" " << quad->GetEid(2)
2577 <<
" " << quad->GetEid(3);
2578 TiXmlElement *
q =
new TiXmlElement(tag);
2579 q->SetAttribute(
"ID", i.first);
2580 q->LinkEndChild(
new TiXmlText(s.str()));
2581 faceTag->LinkEndChild(
q);
2587 std::string tag =
"H";
2589 for (
auto &i : hexs)
2591 std::stringstream s;
2593 s << hex->GetFid(0) <<
" " << hex->GetFid(1) <<
" " << hex->GetFid(2)
2594 <<
" " << hex->GetFid(3) <<
" " << hex->GetFid(4) <<
" "
2595 << hex->GetFid(5) <<
" ";
2596 TiXmlElement *h =
new TiXmlElement(tag);
2597 h->SetAttribute(
"ID", i.first);
2598 h->LinkEndChild(
new TiXmlText(s.str()));
2599 elmtTag->LinkEndChild(h);
2605 std::string tag =
"R";
2607 for (
auto &i : pris)
2609 std::stringstream s;
2611 s << prism->GetFid(0) <<
" " << prism->GetFid(1) <<
" "
2612 << prism->GetFid(2) <<
" " << prism->GetFid(3) <<
" "
2613 << prism->GetFid(4) <<
" ";
2614 TiXmlElement *
p =
new TiXmlElement(tag);
2615 p->SetAttribute(
"ID", i.first);
2616 p->LinkEndChild(
new TiXmlText(s.str()));
2617 elmtTag->LinkEndChild(
p);
2623 std::string tag =
"P";
2625 for (
auto &i : pyrs)
2627 std::stringstream s;
2629 s << pyr->GetFid(0) <<
" " << pyr->GetFid(1) <<
" " << pyr->GetFid(2)
2630 <<
" " << pyr->GetFid(3) <<
" " << pyr->GetFid(4) <<
" ";
2631 TiXmlElement *
p =
new TiXmlElement(tag);
2632 p->SetAttribute(
"ID", i.first);
2633 p->LinkEndChild(
new TiXmlText(s.str()));
2634 elmtTag->LinkEndChild(
p);
2640 std::string tag =
"A";
2642 for (
auto &i : tets)
2644 std::stringstream s;
2646 s << tet->GetFid(0) <<
" " << tet->GetFid(1) <<
" " << tet->GetFid(2)
2647 <<
" " << tet->GetFid(3) <<
" ";
2648 TiXmlElement *t =
new TiXmlElement(tag);
2649 t->SetAttribute(
"ID", i.first);
2650 t->LinkEndChild(
new TiXmlText(s.str()));
2651 elmtTag->LinkEndChild(t);
2658 TiXmlElement *curveTag =
new TiXmlElement(
"CURVED");
2659 CurveMap::iterator curveIt;
2662 for (curveIt = edges.begin(); curveIt != edges.end(); ++curveIt)
2665 TiXmlElement *c =
new TiXmlElement(
"E");
2666 std::stringstream s;
2669 for (
int j = 0; j < curve->m_points.size(); ++j)
2672 s << std::scientific << (*p)(0) <<
" " << (*
p)(1) <<
" " << (*
p)(2)
2676 c->SetAttribute(
"ID", curveId++);
2677 c->SetAttribute(
"EDGEID", curve->m_curveID);
2678 c->SetAttribute(
"NUMPOINTS", curve->m_points.size());
2680 c->LinkEndChild(
new TiXmlText(s.str()));
2681 curveTag->LinkEndChild(c);
2684 for (curveIt = faces.begin(); curveIt != faces.end(); ++curveIt)
2687 TiXmlElement *c =
new TiXmlElement(
"F");
2688 std::stringstream s;
2691 for (
int j = 0; j < curve->m_points.size(); ++j)
2694 s << std::scientific << (*p)(0) <<
" " << (*
p)(1) <<
" " << (*
p)(2)
2698 c->SetAttribute(
"ID", curveId++);
2699 c->SetAttribute(
"FACEID", curve->m_curveID);
2700 c->SetAttribute(
"NUMPOINTS", curve->m_points.size());
2702 c->LinkEndChild(
new TiXmlText(s.str()));
2703 curveTag->LinkEndChild(c);
2706 geomTag->LinkEndChild(curveTag);
2710 std::map<int, std::string> &compLabels)
2712 auto &compositesLabels =
m_meshGraph->GetCompositesLabels();
2713 TiXmlElement *compTag =
new TiXmlElement(
"COMPOSITE");
2715 for (
auto &cIt : comps)
2717 if (cIt.second->m_geomVec.size() == 0)
2722 TiXmlElement *c =
new TiXmlElement(
"C");
2723 c->SetAttribute(
"ID", cIt.first);
2724 if (!compositesLabels[cIt.first].empty())
2726 c->SetAttribute(
"NAME", compLabels[cIt.first]);
2729 compTag->LinkEndChild(c);
2732 geomTag->LinkEndChild(compTag);
2736 std::map<int, CompositeMap> &domainMap)
2738 TiXmlElement *domTag =
new TiXmlElement(
"DOMAIN");
2740 std::vector<unsigned int> idxList;
2741 for (
auto &domain : domainMap)
2743 TiXmlElement *c =
new TiXmlElement(
"D");
2745 std::stringstream s;
2750 for (
const auto &elem : domain.second)
2752 idxList.push_back(elem.first);
2756 c->SetAttribute(
"ID", domain.first);
2757 c->LinkEndChild(
new TiXmlText(s.str()));
2758 domTag->LinkEndChild(c);
2761 geomTag->LinkEndChild(domTag);
2766 int meshDimension =
m_meshGraph->GetMeshDimension();
2767 auto &meshComposites =
m_meshGraph->GetComposites();
2769 TiXmlElement *expTag =
new TiXmlElement(
"EXPANSIONS");
2771 for (
auto it = meshComposites.begin(); it != meshComposites.end(); it++)
2773 if (it->second->m_geomVec[0]->GetShapeDim() == meshDimension)
2775 TiXmlElement *exp =
new TiXmlElement(
"E");
2776 exp->SetAttribute(
"COMPOSITE",
2777 "C[" + std::to_string(it->first) +
"]");
2778 exp->SetAttribute(
"NUMMODES", 4);
2779 exp->SetAttribute(
"TYPE",
"MODIFIED");
2780 exp->SetAttribute(
"FIELDS",
"u");
2782 expTag->LinkEndChild(exp);
2785 root->LinkEndChild(expTag);
2793 const std::string &outfilename,
bool defaultExp,
2796 int meshDimension =
m_meshGraph->GetMeshDimension();
2797 int spaceDimension =
m_meshGraph->GetSpaceDimension();
2801 TiXmlDeclaration *decl =
new TiXmlDeclaration(
"1.0",
"utf-8",
"");
2802 doc.LinkEndChild(decl);
2804 TiXmlElement *root =
new TiXmlElement(
"NEKTAR");
2805 doc.LinkEndChild(root);
2806 TiXmlElement *geomTag =
new TiXmlElement(
"GEOMETRY");
2807 root->LinkEndChild(geomTag);
2815 geomTag->SetAttribute(
"DIM", meshDimension);
2816 geomTag->SetAttribute(
"SPACE", spaceDimension);
2820 geomTag->SetAttribute(
"PARTITION",
m_session->GetComm()->GetRank());
2829 if (meshDimension > 1)
2831 TiXmlElement *faceTag =
2832 new TiXmlElement(meshDimension == 2 ?
"ELEMENT" :
"FACE");
2836 geomTag->LinkEndChild(faceTag);
2838 if (meshDimension > 2)
2840 TiXmlElement *elmtTag =
new TiXmlElement(
"ELEMENT");
2847 geomTag->LinkEndChild(elmtTag);
2863 movement->WriteMovement(root);
2867 doc.SaveFile(outfilename);
2871 std::string outname, std::vector<std::set<unsigned int>> elements,
2872 std::vector<unsigned int> partitions)
2883 std::string dirname = outname +
"_xml";
2884 fs::path pdirname(dirname);
2886 if (!fs::is_directory(dirname))
2888 fs::create_directory(dirname);
2891 ASSERTL0(elements.size() == partitions.size(),
2892 "error in partitioned information");
2894 for (
int i = 0; i < partitions.size(); i++)
2897 TiXmlDeclaration *decl =
new TiXmlDeclaration(
"1.0",
"utf-8",
"");
2898 doc.LinkEndChild(decl);
2900 TiXmlElement *root = doc.FirstChildElement(
"NEKTAR");
2901 TiXmlElement *geomTag;
2906 root =
new TiXmlElement(
"NEKTAR");
2907 doc.LinkEndChild(root);
2909 geomTag =
new TiXmlElement(
"GEOMETRY");
2910 root->LinkEndChild(geomTag);
2915 geomTag = root->FirstChildElement(
"GEOMETRY");
2919 geomTag =
new TiXmlElement(
"GEOMETRY");
2920 root->LinkEndChild(geomTag);
2924 int meshDimension =
m_meshGraph->GetMeshDimension();
2925 int spaceDimension =
m_meshGraph->GetSpaceDimension();
2927 geomTag->SetAttribute(
"DIM", meshDimension);
2928 geomTag->SetAttribute(
"SPACE", spaceDimension);
2929 geomTag->SetAttribute(
"PARTITION", partitions[i]);
2939 auto &prismGeoms =
m_meshGraph->GetAllPrismGeoms();
2941 auto &curvedEdges =
m_meshGraph->GetCurvedEdges();
2942 auto &curvedFaces =
m_meshGraph->GetCurvedFaces();
2943 auto &meshComposites =
m_meshGraph->GetComposites();
2957 std::vector<std::set<unsigned int>> entityIds(4);
2958 entityIds[meshDimension] = elements[i];
2960 switch (meshDimension)
2964 for (
auto &j : entityIds[3])
2967 if (hexGeoms.count(j))
2970 localHex[j] = hexGeoms[j];
2972 else if (pyrGeoms.count(j))
2975 localPyr[j] = pyrGeoms[j];
2977 else if (prismGeoms.count(j))
2980 localPrism[j] = prismGeoms[j];
2982 else if (tetGeoms.count(j))
2985 localTet[j] = tetGeoms[j];
2989 ASSERTL0(
false,
"element in partition not found");
2992 for (
int k = 0; k < g->GetNumFaces(); k++)
2994 entityIds[2].insert(g->GetFid(k));
2996 for (
int k = 0; k < g->GetNumEdges(); k++)
2998 entityIds[1].insert(g->GetEid(k));
3000 for (
int k = 0; k < g->GetNumVerts(); k++)
3002 entityIds[0].insert(g->GetVid(k));
3009 for (
auto &j : entityIds[2])
3012 if (triGeoms.count(j))
3015 localTri[j] = triGeoms[j];
3017 else if (quadGeoms.count(j))
3020 localQuad[j] = quadGeoms[j];
3024 ASSERTL0(
false,
"element in partition not found");
3027 for (
int k = 0; k < g->GetNumEdges(); k++)
3029 entityIds[1].insert(g->GetEid(k));
3031 for (
int k = 0; k < g->GetNumVerts(); k++)
3033 entityIds[0].insert(g->GetVid(k));
3040 for (
auto &j : entityIds[1])
3043 if (segGeoms.count(j))
3046 localEdge[j] = segGeoms[j];
3050 ASSERTL0(
false,
"element in partition not found");
3053 for (
int k = 0; k < g->GetNumVerts(); k++)
3055 entityIds[0].insert(g->GetVid(k));
3061 if (meshDimension > 2)
3063 for (
auto &j : entityIds[2])
3065 if (triGeoms.count(j))
3067 localTri[j] = triGeoms[j];
3069 else if (quadGeoms.count(j))
3071 localQuad[j] = quadGeoms[j];
3080 if (meshDimension > 1)
3082 for (
auto &j : entityIds[1])
3084 if (segGeoms.count(j))
3086 localEdge[j] = segGeoms[j];
3095 for (
auto &j : entityIds[0])
3097 if (vertSet.count(j))
3099 localVert[j] = vertSet[j];
3109 if (meshDimension > 1)
3111 TiXmlElement *faceTag =
3112 new TiXmlElement(meshDimension == 2 ?
"ELEMENT" :
"FACE");
3116 geomTag->LinkEndChild(faceTag);
3118 if (meshDimension > 2)
3120 TiXmlElement *elmtTag =
new TiXmlElement(
"ELEMENT");
3127 geomTag->LinkEndChild(elmtTag);
3130 for (
auto &j : localTri)
3132 if (curvedFaces.count(j.first))
3134 localCurveFace[j.first] = curvedFaces[j.first];
3137 for (
auto &j : localQuad)
3139 if (curvedFaces.count(j.first))
3141 localCurveFace[j.first] = curvedFaces[j.first];
3144 for (
auto &j : localEdge)
3146 if (curvedEdges.count(j.first))
3148 localCurveEdge[j.first] = curvedEdges[j.first];
3155 std::map<int, std::string> localCompLabels;
3157 for (
auto &j : meshComposites)
3160 int dim = j.second->m_geomVec[0]->GetShapeDim();
3162 for (
int k = 0; k < j.second->m_geomVec.size(); k++)
3164 if (entityIds[dim].count(j.second->m_geomVec[k]->GetGlobalID()))
3166 comp->m_geomVec.push_back(j.second->m_geomVec[k]);
3170 if (comp->m_geomVec.size())
3172 localComp[j.first] = comp;
3173 auto compositesLabels =
m_meshGraph->GetCompositesLabels();
3174 if (!compositesLabels[j.first].empty())
3176 localCompLabels[j.first] = compositesLabels[j.first];
3183 std::map<int, CompositeMap> domain;
3184 for (
auto &j : localComp)
3186 for (
auto &dom : globalDomain)
3188 for (
auto &dIt : dom.second)
3190 if (j.first == dIt.first)
3192 domain[dom.first][j.first] = j.second;
3201 if (
m_session->DefinesElement(
"NEKTAR/CONDITIONS"))
3203 std::set<int> vBndRegionIdList;
3204 TiXmlElement *vConditions =
3205 new TiXmlElement(*
m_session->GetElement(
"Nektar/Conditions"));
3206 TiXmlElement *vBndRegions =
3207 vConditions->FirstChildElement(
"BOUNDARYREGIONS");
3210 TiXmlElement *vBndConditions =
3211 vConditions->FirstChildElement(
"BOUNDARYCONDITIONS");
3215 TiXmlElement *vItem;
3221 vConditions->FirstChildElement(
"BOUNDARYREGIONS")
3222 ->FirstChildElement(
"TIMELEVEL") !=
nullptr;
3224 TiXmlElement *vNewBndRegions =
3225 multiLevel ?
new TiXmlElement(
"TIMELEVEL")
3226 :
new TiXmlElement(
"BOUNDARYREGIONS");
3227 vItem = vBndRegions->FirstChildElement();
3228 auto graph_bndRegOrder =
m_meshGraph->GetBndRegionOrdering();
3231 std::string vSeqStr =
3232 vItem->FirstChild()->ToText()->Value();
3233 std::string::size_type indxBeg =
3234 vSeqStr.find_first_of(
'[') + 1;
3235 std::string::size_type indxEnd =
3236 vSeqStr.find_last_of(
']') - 1;
3237 vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
3238 std::vector<unsigned int> vSeq;
3241 std::vector<unsigned int> idxList;
3243 for (
unsigned int i = 0; i < vSeq.size(); ++i)
3245 if (localComp.find(vSeq[i]) != localComp.end())
3247 idxList.push_back(vSeq[i]);
3250 int p = atoi(vItem->Attribute(
"ID"));
3252 std::string vListStr =
3255 if (vListStr.length() == 0)
3257 TiXmlElement *tmp = vItem;
3258 vItem = vItem->NextSiblingElement();
3259 vBndRegions->RemoveChild(tmp);
3263 vListStr =
"C[" + vListStr +
"]";
3264 TiXmlText *vList =
new TiXmlText(vListStr);
3265 TiXmlElement *vNewElement =
new TiXmlElement(
"B");
3266 vNewElement->SetAttribute(
"ID",
p);
3267 vNewElement->LinkEndChild(vList);
3268 vNewBndRegions->LinkEndChild(vNewElement);
3269 vBndRegionIdList.insert(
p);
3270 vItem = vItem->NextSiblingElement();
3275 graph_bndRegOrder[
p] = vSeq;
3280 size_t timeLevel = 0;
3283 vNewBndRegions->SetAttribute(
"VALUE", timeLevel);
3284 vConditions->FirstChildElement(
"BOUNDARYREGIONS")
3285 ->ReplaceChild(vBndRegions, *vNewBndRegions);
3286 vBndRegions = vBndRegions->NextSiblingElement();
3292 vConditions->ReplaceChild(vBndRegions, *vNewBndRegions);
3298 vItem = vBndConditions->FirstChildElement();
3301 std::set<int>::iterator x;
3302 if ((x = vBndRegionIdList.find(atoi(vItem->Attribute(
3303 "REF")))) != vBndRegionIdList.end())
3305 vItem->SetAttribute(
"REF", *x);
3306 vItem = vItem->NextSiblingElement();
3310 TiXmlElement *tmp = vItem;
3311 vItem = vItem->NextSiblingElement();
3312 vBndConditions->RemoveChild(tmp);
3316 root->LinkEndChild(vConditions);
3320 TiXmlElement *vSrc =
3321 m_session->GetElement(
"Nektar")->FirstChildElement();
3324 std::string vName = boost::to_upper_copy(vSrc->ValueStr());
3325 if (vName !=
"GEOMETRY" && vName !=
"CONDITIONS")
3327 root->LinkEndChild(
new TiXmlElement(*vSrc));
3329 vSrc = vSrc->NextSiblingElement();
3335 pad % partitions[i];
3336 fs::path pFilename(pad.str());
3338 fs::path fullpath = pdirname / pFilename;
3345 auto &meshComposites =
m_meshGraph->GetComposites();
3350 for (
auto &c : meshComposites)
3352 bool fillComp =
true;
3353 for (
auto &
d : domain[0])
3355 if (c.second ==
d.second)
3362 std::vector<unsigned int> ids;
3363 for (
auto &elmt : c.second->m_geomVec)
3365 ids.push_back(elmt->GetGlobalID());
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
static void AddInfoTag(TagWriterSharedPtr root, const FieldMetaDataMap &fieldmetadatamap)
Add provenance information to the field metadata map.
Interpreter class for the evaluation of mathematical expressions.
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
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.
bool ModuleExists(tKey idKey)
Checks if a particular module is available.
static void GetXMLElementTimeLevel(TiXmlElement *&element, const size_t timeLevel, const bool enableCheck=true)
Get XML elment time level (Parallel-in-Time)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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.
static const int kNqfaces
static const int kNtfaces
BndRegionOrdering m_bndRegOrder
LibUtilities::SessionReaderSharedPtr m_session
void ReadGeometry(LibUtilities::DomainRangeShPtr rng, bool fillGraph)
std::string GetCompositeString(CompositeSharedPtr comp)
Returns a string representation of a composite.
CompositeDescriptor CreateCompositeDescriptor()
MeshGraphSharedPtr m_meshGraph
CompositeOrdering m_compOrder
std::map< int, MeshEntity > CreateMeshEntities()
Create mesh entities for this graph.
virtual void v_ReadElements1D()
virtual void v_WriteHexs(TiXmlElement *elmtTag, HexGeomMap &hexs)
void ResolveGeomRef2D(const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
void ResolveGeomRef3D(const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
void v_WriteGeometry(const std::string &outfilename, bool defaultExp=false, const LibUtilities::FieldMetaDataMap &metadata=LibUtilities::NullFieldMetaDataMap) override
Write out an XML file containing the GEOMETRY block representing this MeshGraph instance inside a NEK...
virtual void v_ReadElements3D()
virtual void v_WriteCurves(TiXmlElement *geomTag, CurveMap &edges, CurveMap &faces)
void ResolveGeomRef(const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
virtual void v_WritePyrs(TiXmlElement *elmtTag, PyrGeomMap &pyrs)
virtual void v_ReadEdges()
virtual void v_ReadFaces()
void v_PartitionMesh(LibUtilities::SessionReaderSharedPtr session) override
virtual void v_ReadCurves()
void WriteDefaultExpansion(TiXmlElement *root)
void v_ReadGeometry(LibUtilities::DomainRangeShPtr rng, bool fillGraph) override
virtual void v_ReadElements2D()
virtual void v_WriteTets(TiXmlElement *elmtTag, TetGeomMap &tets)
void ResolveGeomRef1D(const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
virtual void v_WriteEdges(TiXmlElement *geomTag, SegGeomMap &edges)
virtual void v_ReadVertices()
virtual void v_WritePrisms(TiXmlElement *elmtTag, PrismGeomMap &pris)
void WriteComposites(TiXmlElement *geomTag, CompositeMap &comps, std::map< int, std::string > &compLabels)
static std::string className
void WriteXMLGeometry(std::string outname, std::vector< std::set< unsigned int > > elements, std::vector< unsigned int > partitions)
virtual void v_WriteVertices(TiXmlElement *geomTag, PointGeomMap &verts)
CompositeOrdering CreateCompositeOrdering()
virtual void v_WriteQuads(TiXmlElement *faceTag, QuadGeomMap &quads)
virtual void v_WriteTris(TiXmlElement *faceTag, TriGeomMap &tris)
void WriteDomain(TiXmlElement *geomTag, std::map< int, CompositeMap > &domain)
static MeshGraphIOSharedPtr create()
static const int kNqfaces
static const int kNtfaces
static const int kNtfaces
static const int kNqfaces
static const int kNqfaces
static const int kNtfaces
static const int kNedges
Get the orientation of face1.
std::map< std::string, std::string > FieldMetaDataMap
const std::string kPointsTypeStr[]
static std::string PortablePath(const fs::path &path)
create portable path on different platforms for std::filesystem path.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< DomainRange > DomainRangeShPtr
std::shared_ptr< XmlTagWriter > XmlTagWriterSharedPtr
static DomainRangeShPtr NullDomainRangeShPtr
@ SIZE_PointsType
Length of enum list.
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
std::map< int, TriGeomSharedPtr > TriGeomMap
std::map< int, std::vector< unsigned int > > CompositeOrdering
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
std::map< int, PyrGeomSharedPtr > PyrGeomMap
std::map< int, QuadGeomSharedPtr > QuadGeomMap
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
MeshPartitionFactory & GetMeshPartitionFactory()
std::shared_ptr< Composite > CompositeSharedPtr
std::map< int, SegGeomSharedPtr > SegGeomMap
std::shared_ptr< Curve > CurveSharedPtr
std::unordered_map< int, CurveSharedPtr > CurveMap
std::shared_ptr< HexGeom > HexGeomSharedPtr
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::map< int, TetGeomSharedPtr > TetGeomMap
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
MeshGraphIOFactory & GetMeshGraphIOFactory()
std::shared_ptr< TetGeom > TetGeomSharedPtr
std::map< int, PrismGeomSharedPtr > PrismGeomMap
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< MeshPartition > MeshPartitionSharedPtr
std::shared_ptr< Geometry > GeometrySharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
std::map< int, HexGeomSharedPtr > HexGeomMap
std::map< int, PointGeomSharedPtr > PointGeomMap
std::map< int, CompositeSharedPtr > CompositeMap
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
std::vector< double > d(NPUPPER *NPUPPER)
std::vector< double > q(NPUPPER *NPUPPER)