46 #include <boost/format.hpp>
54 namespace SpatialDomains
57 std::string MeshGraphXml::className =
59 "IO with Xml geometry");
61 void MeshGraphXml::PartitionMesh(
67 const bool isRoot = comm->TreatAsRankZero();
79 int isPartitioned = 0;
82 if (m_session->DefinesElement(
"Nektar/Geometry"))
84 if (m_session->GetElement(
"Nektar/Geometry")
85 ->Attribute(
"PARTITION"))
87 std::cout <<
"Using pre-partitioned mesh." << std::endl;
92 comm->Bcast(isPartitioned, 0);
100 m_session->InitSession();
107 string partitionerName =
"Metis";
110 partitionerName =
"Scotch";
112 if (session->DefinesCmdLineArgument(
"use-metis"))
114 partitionerName =
"Metis";
116 if (session->DefinesCmdLineArgument(
"use-scotch"))
118 partitionerName =
"Scotch";
124 if (session->DefinesCmdLineArgument(
"part-only") ||
125 session->DefinesCmdLineArgument(
"part-only-overlapping"))
130 "The 'part-only' option should be used in serial.");
137 auto comp = CreateCompositeDescriptor();
141 partitionerName, session, session->GetComm(),
142 m_meshDimension, CreateMeshEntities(), comp);
144 if (session->DefinesCmdLineArgument(
"part-only"))
146 nParts = session->GetCmdLineArgument<
int>(
"part-only");
147 partitioner->PartitionMesh(nParts,
true);
152 session->GetCmdLineArgument<
int>(
"part-only-overlapping");
153 partitioner->PartitionMesh(nParts,
true,
true);
156 vector<set<unsigned int>> elmtIDs;
157 vector<unsigned int> parts(nParts);
158 for (
int i = 0; i < nParts; ++i)
160 vector<unsigned int> elIDs;
161 set<unsigned int> tmp;
162 partitioner->GetElementIDs(i, elIDs);
163 tmp.insert(elIDs.begin(), elIDs.end());
164 elmtIDs.push_back(tmp);
168 this->WriteXMLGeometry(m_session->GetSessionName(), elmtIDs, parts);
170 if (isRoot && session->DefinesCmdLineArgument(
"part-info"))
172 partitioner->PrintPartInfo(std::cout);
179 if (commMesh->GetSize() > 1)
181 int nParts = commMesh->GetSize();
183 if (session->GetSharedFilesystem())
185 vector<unsigned int> keys, vals;
194 m_compOrder = CreateCompositeOrdering();
195 auto comp = CreateCompositeDescriptor();
200 partitionerName, session, session->GetComm(),
201 m_meshDimension, CreateMeshEntities(), comp);
203 partitioner->PartitionMesh(nParts,
true);
205 vector<set<unsigned int>> elmtIDs;
206 vector<unsigned int> parts(nParts);
207 for (i = 0; i < nParts; ++i)
209 vector<unsigned int> elIDs;
210 set<unsigned int> tmp;
211 partitioner->GetElementIDs(i, elIDs);
212 tmp.insert(elIDs.begin(), elIDs.end());
213 elmtIDs.push_back(tmp);
219 this->WriteXMLGeometry(m_session->GetSessionName(), elmtIDs,
227 keys[0] = m_compOrder.size();
228 keys[1] = m_bndRegOrder.size();
229 comm->Bcast(keys, 0);
233 keys.resize(m_compOrder.size());
234 vals.resize(m_compOrder.size());
237 for (
auto &cIt : m_compOrder)
240 vals[i++] = cIt.second.size();
244 comm->Bcast(keys, 0);
245 comm->Bcast(vals, 0);
246 for (
auto &cIt : m_compOrder)
248 comm->Bcast(cIt.second, 0);
253 keys.resize(m_bndRegOrder.size());
254 vals.resize(m_bndRegOrder.size());
257 for (
auto &bIt : m_bndRegOrder)
260 vals[i++] = bIt.second.size();
266 comm->Bcast(keys, 0);
270 comm->Bcast(vals, 0);
272 for (
auto &bIt : m_bndRegOrder)
274 comm->Bcast(bIt.second, 0);
277 if (session->DefinesCmdLineArgument(
"part-info"))
279 partitioner->PrintPartInfo(std::cout);
285 comm->Bcast(keys, 0);
287 int cmpSize = keys[0];
288 int bndSize = keys[1];
290 keys.resize(cmpSize);
291 vals.resize(cmpSize);
292 comm->Bcast(keys, 0);
293 comm->Bcast(vals, 0);
295 for (
int i = 0; i < keys.size(); ++i)
297 vector<unsigned int> tmp(vals[i]);
299 m_compOrder[keys[i]] = tmp;
302 keys.resize(bndSize);
303 vals.resize(bndSize);
306 comm->Bcast(keys, 0);
310 comm->Bcast(vals, 0);
312 for (
int i = 0; i < keys.size(); ++i)
314 vector<unsigned int> tmp(vals[i]);
316 m_bndRegOrder[keys[i]] = tmp;
322 m_session->InitSession();
325 m_compOrder = CreateCompositeOrdering();
326 auto comp = CreateCompositeDescriptor();
333 partitionerName, session, session->GetComm(),
334 m_meshDimension, CreateMeshEntities(), comp);
336 partitioner->PartitionMesh(nParts,
false);
338 vector<unsigned int> parts(1), tmp;
339 parts[0] = commMesh->GetRank();
340 vector<set<unsigned int>> elIDs(1);
341 partitioner->GetElementIDs(parts[0], tmp);
342 elIDs[0].insert(tmp.begin(), tmp.end());
343 this->WriteXMLGeometry(session->GetSessionName(), elIDs, parts);
345 if (m_session->DefinesCmdLineArgument(
"part-info") && isRoot)
347 partitioner->PrintPartInfo(std::cout);
354 std::string dirname = m_session->GetSessionName() +
"_xml";
355 fs::path pdirname(dirname);
357 pad % comm->GetRowComm()->GetRank();
358 fs::path pFilename(pad.str());
359 fs::path fullpath = pdirname / pFilename;
361 std::vector<std::string> filenames = {
363 m_session->InitSession(filenames);
370 m_session->InitSession();
380 m_curvedEdges.clear();
381 m_curvedFaces.clear();
387 m_prismGeoms.clear();
389 m_meshComposites.clear();
390 m_compositesLabels.clear();
392 m_expansionMapShPtrMap.clear();
394 m_faceToElMap.clear();
397 m_xmlGeom = m_session->GetElement(
"NEKTAR/GEOMETRY");
401 TiXmlAttribute *attr = m_xmlGeom->FirstAttribute();
406 m_meshPartitioned =
false;
408 m_spaceDimension = 3;
412 std::string attrName(attr->Name());
413 if (attrName ==
"DIM")
415 err = attr->QueryIntValue(&m_meshDimension);
416 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read mesh dimension.");
418 else if (attrName ==
"SPACE")
420 err = attr->QueryIntValue(&m_spaceDimension);
421 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read space dimension.");
423 else if (attrName ==
"PARTITION")
425 err = attr->QueryIntValue(&m_partition);
426 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read partition.");
427 m_meshPartitioned =
true;
431 std::string errstr(
"Unknown attribute: ");
440 ASSERTL0(m_meshDimension <= m_spaceDimension,
441 "Mesh dimension greater than space dimension");
445 if (m_meshDimension >= 2)
448 if (m_meshDimension == 3)
459 MeshGraph::FillGraph();
463 void MeshGraphXml::ReadVertices()
466 TiXmlElement *element = m_xmlGeom->FirstChildElement(
"VERTEX");
467 ASSERTL0(element,
"Unable to find mesh VERTEX tag in file.");
474 const char *xscal = element->Attribute(
"XSCALE");
481 std::string xscalstr = xscal;
483 xscale = expEvaluator.
Evaluate(expr_id);
486 const char *yscal = element->Attribute(
"YSCALE");
493 std::string yscalstr = yscal;
495 yscale = expEvaluator.
Evaluate(expr_id);
498 const char *zscal = element->Attribute(
"ZSCALE");
505 std::string zscalstr = zscal;
507 zscale = expEvaluator.
Evaluate(expr_id);
515 const char *xmov = element->Attribute(
"XMOVE");
522 std::string xmovstr = xmov;
524 xmove = expEvaluator.
Evaluate(expr_id);
527 const char *ymov = element->Attribute(
"YMOVE");
534 std::string ymovstr = ymov;
536 ymove = expEvaluator.
Evaluate(expr_id);
539 const char *zmov = element->Attribute(
"ZMOVE");
546 std::string zmovstr = zmov;
548 zmove = expEvaluator.
Evaluate(expr_id);
551 TiXmlElement *vertex = element->FirstChildElement(
"V");
554 int nextVertexNumber = -1;
560 TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
561 std::string attrName(vertexAttr->Name());
564 (std::string(
"Unknown attribute name: ") + attrName).c_str());
566 int err = vertexAttr->QueryIntValue(&indx);
567 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
570 std::string vertexBodyStr;
572 TiXmlNode *vertexBody = vertex->FirstChild();
577 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
579 vertexBodyStr += vertexBody->ToText()->Value();
580 vertexBodyStr +=
" ";
583 vertexBody = vertexBody->NextSibling();
587 "Vertex definitions must contain vertex data.");
591 std::istringstream vertexDataStrm(vertexBodyStr.c_str());
595 while (!vertexDataStrm.fail())
597 vertexDataStrm >> xval >> yval >> zval;
599 xval = xval * xscale + xmove;
600 yval = yval * yscale + ymove;
601 zval = zval * zscale + zmove;
606 if (!vertexDataStrm.fail())
610 m_spaceDimension, indx, xval, yval, zval));
611 m_vertSet[indx] = vert;
617 ASSERTL0(
false,
"Unable to read VERTEX data.");
620 vertex = vertex->NextSiblingElement(
"V");
624 void MeshGraphXml::ReadCurves()
628 TiXmlElement *element = m_xmlGeom->FirstChildElement(
"VERTEX");
629 ASSERTL0(element,
"Unable to find mesh VERTEX tag in file.");
634 const char *xscal = element->Attribute(
"XSCALE");
641 std::string xscalstr = xscal;
643 xscale = expEvaluator.
Evaluate(expr_id);
646 const char *yscal = element->Attribute(
"YSCALE");
653 std::string yscalstr = yscal;
655 yscale = expEvaluator.
Evaluate(expr_id);
658 const char *zscal = element->Attribute(
"ZSCALE");
665 std::string zscalstr = zscal;
667 zscale = expEvaluator.
Evaluate(expr_id);
675 const char *xmov = element->Attribute(
"XMOVE");
682 std::string xmovstr = xmov;
684 xmove = expEvaluator.
Evaluate(expr_id);
687 const char *ymov = element->Attribute(
"YMOVE");
694 std::string ymovstr = ymov;
696 ymove = expEvaluator.
Evaluate(expr_id);
699 const char *zmov = element->Attribute(
"ZMOVE");
706 std::string zmovstr = zmov;
708 zmove = expEvaluator.
Evaluate(expr_id);
714 TiXmlElement *field = m_xmlGeom->FirstChildElement(
"CURVED");
725 TiXmlElement *edgelement = field->FirstChildElement(
"E");
727 int edgeindx, edgeid;
728 int nextEdgeNumber = -1;
735 std::string edge(edgelement->ValueStr());
737 (std::string(
"Unknown 3D curve type:") + edge).c_str());
740 err = edgelement->QueryIntAttribute(
"ID", &edgeindx);
741 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read curve attribute ID.");
744 err = edgelement->QueryIntAttribute(
"EDGEID", &edgeid);
746 "Unable to read curve attribute EDGEID.");
749 std::string elementStr;
750 TiXmlNode *elementChild = edgelement->FirstChild();
755 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
757 elementStr += elementChild->ToText()->ValueStr();
760 elementChild = elementChild->NextSibling();
763 ASSERTL0(!elementStr.empty(),
"Unable to read curve description body.");
771 std::string typeStr = edgelement->Attribute(
"TYPE");
772 ASSERTL0(!typeStr.empty(),
"TYPE must be specified in "
773 "points definition");
777 const std::string *endStr =
779 const std::string *ptsStr =
std::find(begStr, endStr, typeStr);
781 ASSERTL0(ptsStr != endStr,
"Invalid points type.");
785 err = edgelement->QueryIntAttribute(
"NUMPOINTS", &numPts);
787 "Unable to read curve attribute NUMPOINTS.");
793 std::istringstream elementDataStrm(elementStr.c_str());
796 while (!elementDataStrm.fail())
798 elementDataStrm >> xval >> yval >> zval;
800 xval = xval * xscale + xmove;
801 yval = yval * yscale + ymove;
802 zval = zval * zscale + zmove;
807 if (!elementDataStrm.fail())
811 m_meshDimension, edgeindx, xval, yval, zval));
813 curve->m_points.push_back(vert);
820 (std::string(
"Unable to read curve data for EDGE: ") +
825 ASSERTL0(curve->m_points.size() == numPts,
826 "Number of points specificed by attribute "
827 "NUMPOINTS is different from number of points "
828 "in list (edgeid = " +
829 boost::lexical_cast<string>(edgeid));
831 m_curvedEdges[edgeid] = curve;
833 edgelement = edgelement->NextSiblingElement(
"E");
839 TiXmlElement *facelement = field->FirstChildElement(
"F");
840 int faceindx, faceid;
844 std::string face(facelement->ValueStr());
846 (std::string(
"Unknown 3D curve type: ") + face).c_str());
849 err = facelement->QueryIntAttribute(
"ID", &faceindx);
850 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read curve attribute ID.");
853 err = facelement->QueryIntAttribute(
"FACEID", &faceid);
855 "Unable to read curve attribute FACEID.");
858 std::string elementStr;
859 TiXmlNode *elementChild = facelement->FirstChild();
864 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
866 elementStr += elementChild->ToText()->ValueStr();
869 elementChild = elementChild->NextSibling();
872 ASSERTL0(!elementStr.empty(),
"Unable to read curve description body.");
878 std::string typeStr = facelement->Attribute(
"TYPE");
879 ASSERTL0(!typeStr.empty(),
"TYPE must be specified in "
880 "points definition");
883 const std::string *endStr =
885 const std::string *ptsStr =
std::find(begStr, endStr, typeStr);
887 ASSERTL0(ptsStr != endStr,
"Invalid points type.");
890 std::string numptsStr = facelement->Attribute(
"NUMPOINTS");
892 "NUMPOINTS must be specified in points definition");
901 ASSERTL0(numPts >= 3,
"NUMPOINTS for face must be greater than 2");
905 ASSERTL0(ptsStr != endStr,
"Invalid points type.");
910 std::istringstream elementDataStrm(elementStr.c_str());
913 while (!elementDataStrm.fail())
915 elementDataStrm >> xval >> yval >> zval;
921 if (!elementDataStrm.fail())
925 m_meshDimension, faceindx, xval, yval, zval));
926 curve->m_points.push_back(vert);
933 (std::string(
"Unable to read curve data for FACE: ") +
937 m_curvedFaces[faceid] = curve;
939 facelement = facelement->NextSiblingElement(
"F");
944 void MeshGraphXml::ReadDomain()
946 TiXmlElement *domain = NULL;
948 domain = m_xmlGeom->FirstChildElement(
"DOMAIN");
950 ASSERTL0(domain,
"Unable to find DOMAIN tag in file.");
954 TiXmlElement *multidomains = domain->FirstChildElement(
"D");
961 int err = multidomains->QueryIntAttribute(
"ID", &indx);
963 "Unable to read attribute ID in Domain.");
965 TiXmlNode *elementChild = multidomains->FirstChild();
966 while (elementChild &&
967 elementChild->Type() != TiXmlNode::TINYXML_TEXT)
969 elementChild = elementChild->NextSibling();
972 ASSERTL0(elementChild,
"Unable to read DOMAIN body.");
973 std::string elementStr = elementChild->ToText()->ValueStr();
975 elementStr = elementStr.substr(elementStr.find_first_not_of(
" "));
977 std::string::size_type indxBeg = elementStr.find_first_of(
'[') + 1;
978 std::string::size_type indxEnd = elementStr.find_last_of(
']') - 1;
979 std::string indxStr =
980 elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
984 "Unable to read domain's composite index (index missing?).");
988 map<int, CompositeSharedPtr> unrollDomain;
989 GetCompositeList(indxStr, unrollDomain);
990 m_domain[indx] = unrollDomain;
994 "Unable to obtain domain's referenced composite: ") +
999 multidomains = multidomains->NextSiblingElement(
"D");
1006 TiXmlNode *elementChild = domain->FirstChild();
1007 while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1009 elementChild = elementChild->NextSibling();
1012 ASSERTL0(elementChild,
"Unable to read DOMAIN body.");
1013 std::string elementStr = elementChild->ToText()->ValueStr();
1015 elementStr = elementStr.substr(elementStr.find_first_not_of(
" "));
1017 std::string::size_type indxBeg = elementStr.find_first_of(
'[') + 1;
1018 std::string::size_type indxEnd = elementStr.find_last_of(
']') - 1;
1019 std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
1022 "Unable to read domain's composite index (index missing?).");
1026 map<int, CompositeSharedPtr> fullDomain;
1027 GetCompositeList(indxStr, fullDomain);
1028 m_domain[0] = fullDomain;
1031 !m_domain[0].empty(),
1032 (std::string(
"Unable to obtain domain's referenced composite: ") +
1038 void MeshGraphXml::ReadEdges()
1040 CurveMap::iterator it;
1043 TiXmlElement *field = m_xmlGeom->FirstChildElement(
"EDGE");
1045 ASSERTL0(field,
"Unable to find EDGE tag in file.");
1050 TiXmlElement *edge = field->FirstChildElement(
"E");
1058 std::string edgeStr;
1063 int err = edge->QueryIntAttribute(
"ID", &indx);
1064 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read edge attribute ID.");
1066 TiXmlNode *child = edge->FirstChild();
1068 if (child->Type() == TiXmlNode::TINYXML_TEXT)
1070 edgeStr += child->ToText()->ValueStr();
1074 int vertex1, vertex2;
1075 std::istringstream edgeDataStrm(edgeStr.c_str());
1079 while (!edgeDataStrm.fail())
1081 edgeDataStrm >> vertex1 >> vertex2;
1088 if (!edgeDataStrm.fail())
1091 GetVertex(vertex2)};
1093 it = m_curvedEdges.find(indx);
1095 if (it == m_curvedEdges.end())
1098 indx, m_spaceDimension, vertices);
1103 indx, m_spaceDimension, vertices, it->second);
1106 m_segGeoms[indx] = edge;
1114 (std::string(
"Unable to read edge data: ") + edgeStr).c_str());
1117 edge = edge->NextSiblingElement(
"E");
1121 void MeshGraphXml::ReadFaces()
1124 TiXmlElement *field = m_xmlGeom->FirstChildElement(
"FACE");
1126 ASSERTL0(field,
"Unable to find FACE tag in file.");
1132 TiXmlElement *element = field->FirstChildElement();
1133 CurveMap::iterator it;
1137 std::string elementType(element->ValueStr());
1139 ASSERTL0(elementType ==
"Q" || elementType ==
"T",
1140 (std::string(
"Unknown 3D face type: ") + elementType).c_str());
1144 int err = element->QueryIntAttribute(
"ID", &indx);
1145 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read face attribute ID.");
1148 it = m_curvedFaces.find(indx);
1151 TiXmlNode *elementChild = element->FirstChild();
1152 std::string elementStr;
1153 while (elementChild)
1155 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1157 elementStr += elementChild->ToText()->ValueStr();
1159 elementChild = elementChild->NextSibling();
1162 ASSERTL0(!elementStr.empty(),
"Unable to read face description body.");
1166 if (elementType ==
"T")
1169 int edge1, edge2, edge3;
1170 std::istringstream elementDataStrm(elementStr.c_str());
1174 elementDataStrm >> edge1;
1175 elementDataStrm >> edge2;
1176 elementDataStrm >> edge3;
1179 !elementDataStrm.fail(),
1180 (std::string(
"Unable to read face data for TRIANGLE: ") +
1186 GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3)};
1190 if (it == m_curvedFaces.end())
1198 indx, edges, it->second);
1201 trigeom->SetGlobalID(indx);
1203 m_triGeoms[indx] = trigeom;
1209 (std::string(
"Unable to read face data for TRIANGLE: ") +
1214 else if (elementType ==
"Q")
1217 int edge1, edge2, edge3, edge4;
1218 std::istringstream elementDataStrm(elementStr.c_str());
1222 elementDataStrm >> edge1;
1223 elementDataStrm >> edge2;
1224 elementDataStrm >> edge3;
1225 elementDataStrm >> edge4;
1228 (std::string(
"Unable to read face data for QUAD: ") +
1234 GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3),
1239 if (it == m_curvedFaces.end())
1247 indx, edges, it->second);
1249 quadgeom->SetGlobalID(indx);
1251 m_quadGeoms[indx] = quadgeom;
1256 (std::string(
"Unable to read face data for QUAD: ") +
1262 element = element->NextSiblingElement();
1266 void MeshGraphXml::ReadElements()
1268 switch (m_meshDimension)
1282 void MeshGraphXml::ReadElements1D()
1284 TiXmlElement *field = NULL;
1287 field = m_xmlGeom->FirstChildElement(
"ELEMENT");
1289 ASSERTL0(field,
"Unable to find ELEMENT tag in file.");
1294 TiXmlElement *segment = field->FirstChildElement(
"S");
1295 CurveMap::iterator it;
1300 int err = segment->QueryIntAttribute(
"ID", &indx);
1301 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read element attribute ID.");
1303 TiXmlNode *elementChild = segment->FirstChild();
1304 while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1306 elementChild = elementChild->NextSibling();
1309 ASSERTL0(elementChild,
"Unable to read element description body.");
1310 std::string elementStr = elementChild->ToText()->ValueStr();
1315 int vertex1, vertex2;
1316 std::istringstream elementDataStrm(elementStr.c_str());
1320 elementDataStrm >> vertex1;
1321 elementDataStrm >> vertex2;
1324 (std::string(
"Unable to read element data for SEGMENT: ") +
1329 GetVertex(vertex2)};
1331 it = m_curvedEdges.find(indx);
1333 if (it == m_curvedEdges.end())
1336 indx, m_spaceDimension, vertices);
1337 seg->SetGlobalID(indx);
1342 indx, m_spaceDimension, vertices, it->second);
1343 seg->SetGlobalID(indx);
1345 seg->SetGlobalID(indx);
1346 m_segGeoms[indx] = seg;
1351 (std::string(
"Unable to read element data for segment: ") +
1356 segment = segment->NextSiblingElement(
"S");
1360 void MeshGraphXml::ReadElements2D()
1363 TiXmlElement *field = m_xmlGeom->FirstChildElement(
"ELEMENT");
1365 ASSERTL0(field,
"Unable to find ELEMENT tag in file.");
1368 CurveMap::iterator it;
1373 TiXmlElement *element = field->FirstChildElement();
1377 std::string elementType(element->ValueStr());
1380 elementType ==
"Q" || elementType ==
"T",
1381 (std::string(
"Unknown 2D element type: ") + elementType).c_str());
1385 int err = element->QueryIntAttribute(
"ID", &indx);
1386 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read element attribute ID.");
1388 it = m_curvedFaces.find(indx);
1391 TiXmlNode *elementChild = element->FirstChild();
1392 std::string elementStr;
1393 while (elementChild)
1395 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1397 elementStr += elementChild->ToText()->ValueStr();
1399 elementChild = elementChild->NextSibling();
1403 "Unable to read element description body.");
1407 if (elementType ==
"T")
1410 int edge1, edge2, edge3;
1411 std::istringstream elementDataStrm(elementStr.c_str());
1415 elementDataStrm >> edge1;
1416 elementDataStrm >> edge2;
1417 elementDataStrm >> edge3;
1420 !elementDataStrm.fail(),
1421 (std::string(
"Unable to read element data for TRIANGLE: ") +
1427 GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3)};
1430 if (it == m_curvedFaces.end())
1438 indx, edges, it->second);
1440 trigeom->SetGlobalID(indx);
1442 m_triGeoms[indx] = trigeom;
1448 (std::string(
"Unable to read element data for TRIANGLE: ") +
1453 else if (elementType ==
"Q")
1456 int edge1, edge2, edge3, edge4;
1457 std::istringstream elementDataStrm(elementStr.c_str());
1461 elementDataStrm >> edge1;
1462 elementDataStrm >> edge2;
1463 elementDataStrm >> edge3;
1464 elementDataStrm >> edge4;
1467 !elementDataStrm.fail(),
1468 (std::string(
"Unable to read element data for QUAD: ") +
1474 GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3),
1478 if (it == m_curvedFaces.end())
1486 indx, edges, it->second);
1488 quadgeom->SetGlobalID(indx);
1490 m_quadGeoms[indx] = quadgeom;
1496 (std::string(
"Unable to read element data for QUAD: ") +
1502 element = element->NextSiblingElement();
1506 void MeshGraphXml::ReadElements3D()
1509 TiXmlElement *field = m_xmlGeom->FirstChildElement(
"ELEMENT");
1511 ASSERTL0(field,
"Unable to find ELEMENT tag in file.");
1516 TiXmlElement *element = field->FirstChildElement();
1520 std::string elementType(element->ValueStr());
1524 elementType ==
"A" || elementType ==
"P" || elementType ==
"R" ||
1526 (std::string(
"Unknown 3D element type: ") + elementType).c_str());
1530 int err = element->QueryIntAttribute(
"ID", &indx);
1531 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read element attribute ID.");
1534 TiXmlNode *elementChild = element->FirstChild();
1535 std::string elementStr;
1536 while (elementChild)
1538 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1540 elementStr += elementChild->ToText()->ValueStr();
1542 elementChild = elementChild->NextSibling();
1546 "Unable to read element description body.");
1548 std::istringstream elementDataStrm(elementStr.c_str());
1554 if (elementType ==
"A")
1559 const int kNfaces = TetGeom::kNfaces;
1560 const int kNtfaces = TetGeom::kNtfaces;
1561 const int kNqfaces = TetGeom::kNqfaces;
1568 std::stringstream errorstring;
1569 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1570 <<
" triangle face(s), and " << kNqfaces
1571 <<
" quadrilateral face(s).";
1572 for (
int i = 0; i < kNfaces; i++)
1575 elementDataStrm >> faceID;
1581 std::stringstream errorstring;
1582 errorstring <<
"Element " << indx
1583 <<
" has invalid face: " << faceID;
1584 ASSERTL0(
false, errorstring.str().c_str());
1588 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1589 tfaces[Ntfaces++] = static_pointer_cast<TriGeom>(face);
1591 else if (face->GetShapeType() ==
1594 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1602 "Unable to read element data for TETRAHEDRON: ") +
1605 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1606 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1611 m_tetGeoms[indx] = tetgeom;
1612 PopulateFaceToElMap(tetgeom, kNfaces);
1618 "Unable to read element data for TETRAHEDRON: ") +
1624 else if (elementType ==
"P")
1629 const int kNfaces = PyrGeom::kNfaces;
1630 const int kNtfaces = PyrGeom::kNtfaces;
1631 const int kNqfaces = PyrGeom::kNqfaces;
1639 std::stringstream errorstring;
1640 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1641 <<
" triangle face(s), and " << kNqfaces
1642 <<
" quadrilateral face(s).";
1643 for (
int i = 0; i < kNfaces; i++)
1646 elementDataStrm >> faceID;
1652 std::stringstream errorstring;
1653 errorstring <<
"Element " << indx
1654 <<
" has invalid face: " << faceID;
1655 ASSERTL0(
false, errorstring.str().c_str());
1659 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1660 faces[Nfaces++] = static_pointer_cast<TriGeom>(face);
1663 else if (face->GetShapeType() ==
1666 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1667 faces[Nfaces++] = static_pointer_cast<QuadGeom>(face);
1675 !elementDataStrm.fail(),
1676 (std::string(
"Unable to read element data for PYRAMID: ") +
1679 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1680 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1685 m_pyrGeoms[indx] = pyrgeom;
1686 PopulateFaceToElMap(pyrgeom, kNfaces);
1692 (std::string(
"Unable to read element data for PYRAMID: ") +
1698 else if (elementType ==
"R")
1703 const int kNfaces = PrismGeom::kNfaces;
1704 const int kNtfaces = PrismGeom::kNtfaces;
1705 const int kNqfaces = PrismGeom::kNqfaces;
1713 std::stringstream errorstring;
1714 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1715 <<
" triangle face(s), and " << kNqfaces
1716 <<
" quadrilateral face(s).";
1718 for (
int i = 0; i < kNfaces; i++)
1721 elementDataStrm >> faceID;
1727 std::stringstream errorstring;
1728 errorstring <<
"Element " << indx
1729 <<
" has invalid face: " << faceID;
1730 ASSERTL0(
false, errorstring.str().c_str());
1734 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1736 std::static_pointer_cast<TriGeom>(face);
1739 else if (face->GetShapeType() ==
1742 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1744 std::static_pointer_cast<QuadGeom>(face);
1752 !elementDataStrm.fail(),
1753 (std::string(
"Unable to read element data for PRISM: ") +
1756 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1757 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1762 m_prismGeoms[indx] = prismgeom;
1763 PopulateFaceToElMap(prismgeom, kNfaces);
1769 (std::string(
"Unable to read element data for PRISM: ") +
1775 else if (elementType ==
"H")
1780 const int kNfaces = HexGeom::kNfaces;
1781 const int kNtfaces = HexGeom::kNtfaces;
1782 const int kNqfaces = HexGeom::kNqfaces;
1790 std::stringstream errorstring;
1791 errorstring <<
"Element " << indx <<
" must have " << kNtfaces
1792 <<
" triangle face(s), and " << kNqfaces
1793 <<
" quadrilateral face(s).";
1794 for (
int i = 0; i < kNfaces; i++)
1797 elementDataStrm >> faceID;
1803 std::stringstream errorstring;
1804 errorstring <<
"Element " << indx
1805 <<
" has invalid face: " << faceID;
1806 ASSERTL0(
false, errorstring.str().c_str());
1810 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1814 else if (face->GetShapeType() ==
1817 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1819 std::static_pointer_cast<QuadGeom>(face);
1827 "Unable to read element data for HEXAHEDRAL: ") +
1830 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1831 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1836 m_hexGeoms[indx] = hexgeom;
1837 PopulateFaceToElMap(hexgeom, kNfaces);
1843 "Unable to read element data for HEXAHEDRAL: ") +
1849 element = element->NextSiblingElement();
1853 void MeshGraphXml::ReadComposites()
1855 TiXmlElement *field = NULL;
1858 field = m_xmlGeom->FirstChildElement(
"COMPOSITE");
1860 ASSERTL0(field,
"Unable to find COMPOSITE tag in file.");
1862 TiXmlElement *node = field->FirstChildElement(
"C");
1865 int nextCompositeNumber = -1;
1872 nextCompositeNumber++;
1875 int err = node->QueryIntAttribute(
"ID", &indx);
1876 ASSERTL0(err == TIXML_SUCCESS,
"Unable to read attribute ID.");
1880 TiXmlNode *compositeChild = node->FirstChild();
1885 while (compositeChild &&
1886 compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
1888 compositeChild = compositeChild->NextSibling();
1891 ASSERTL0(compositeChild,
"Unable to read composite definition body.");
1892 std::string compositeStr = compositeChild->ToText()->ValueStr();
1896 std::istringstream compositeDataStrm(compositeStr.c_str());
1901 std::string prevCompositeElementStr;
1903 while (!compositeDataStrm.fail())
1905 std::string compositeElementStr;
1906 compositeDataStrm >> compositeElementStr;
1908 if (!compositeDataStrm.fail())
1916 m_meshComposites[indx] = curVector;
1919 if (compositeElementStr.length() > 0)
1921 ResolveGeomRef(prevCompositeElementStr,
1922 compositeElementStr,
1923 m_meshComposites[indx]);
1925 prevCompositeElementStr = compositeElementStr;
1933 (std::string(
"Unable to read COMPOSITE data for composite: ") +
1940 err = node->QueryStringAttribute(
"NAME", &
name);
1941 if (err == TIXML_SUCCESS)
1943 m_compositesLabels[indx] =
name;
1947 node = node->NextSiblingElement(
"C");
1951 "At least one composite must be specified.");
1954 void MeshGraphXml::ResolveGeomRef(
const std::string &prevToken,
1955 const std::string &token,
1958 switch (m_meshDimension)
1961 ResolveGeomRef1D(prevToken, token, composite);
1964 ResolveGeomRef2D(prevToken, token, composite);
1967 ResolveGeomRef3D(prevToken, token, composite);
1972 void MeshGraphXml::ResolveGeomRef1D(
const std::string &prevToken,
1973 const std::string &token,
1978 std::istringstream tokenStream(token);
1979 std::istringstream prevTokenStream(prevToken);
1984 tokenStream >> type;
1986 std::string::size_type indxBeg = token.find_first_of(
'[') + 1;
1987 std::string::size_type indxEnd = token.find_last_of(
']') - 1;
1991 (std::string(
"Error reading index definition:") + token).c_str());
1993 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
1995 typedef vector<unsigned int> SeqVectorType;
1996 SeqVectorType seqVector;
1998 if (!ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector))
2001 (std::string(
"Ill-formed sequence definition: ") + indxStr)
2005 prevTokenStream >> prevType;
2008 bool validSequence =
2009 (prevToken.empty() ||
2010 (type ==
'V' && prevType ==
'V') ||
2011 (type ==
'S' && prevType ==
'S'));
2014 (std::string(
"Invalid combination of composite items: ") +
2015 type +
" and " + prevType +
".")
2021 for (SeqVectorType::iterator iter = seqVector.begin();
2022 iter != seqVector.end(); ++iter)
2024 if (m_vertSet.find(*iter) == m_vertSet.end())
2026 char errStr[16] =
"";
2027 ::sprintf(errStr,
"%d", *iter);
2029 ErrorUtil::ewarning,
2030 (std::string(
"Unknown vertex index: ") + errStr)
2035 composite->m_geomVec.push_back(m_vertSet[*iter]);
2041 for (SeqVectorType::iterator iter = seqVector.begin();
2042 iter != seqVector.end(); ++iter)
2044 if (m_segGeoms.find(*iter) == m_segGeoms.end())
2046 char errStr[16] =
"";
2047 ::sprintf(errStr,
"%d", *iter);
2049 ErrorUtil::ewarning,
2050 (std::string(
"Unknown segment index: ") + errStr)
2055 composite->m_geomVec.push_back(m_segGeoms[*iter]);
2062 (std::string(
"Unrecognized composite token: ") + token)
2069 (std::string(
"Problem processing composite token: ") + token)
2076 void MeshGraphXml::ResolveGeomRef2D(
const std::string &prevToken,
2077 const std::string &token,
2082 std::istringstream tokenStream(token);
2083 std::istringstream prevTokenStream(prevToken);
2088 tokenStream >> type;
2090 std::string::size_type indxBeg = token.find_first_of(
'[') + 1;
2091 std::string::size_type indxEnd = token.find_last_of(
']') - 1;
2095 (std::string(
"Error reading index definition:") + token).c_str());
2097 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2098 std::vector<unsigned int> seqVector;
2099 std::vector<unsigned int>::iterator seqIter;
2101 bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
2104 (std::string(
"Error reading composite elements: ") + indxStr)
2107 prevTokenStream >> prevType;
2110 bool validSequence =
2111 (prevToken.empty() ||
2112 (type ==
'V' && prevType ==
'V') ||
2113 (type ==
'E' && prevType ==
'E') ||
2114 ((type ==
'T' || type ==
'Q') &&
2115 (prevType ==
'T' || prevType ==
'Q')));
2118 (std::string(
"Invalid combination of composite items: ") +
2119 type +
" and " + prevType +
".")
2125 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2128 if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
2130 char errStr[16] =
"";
2131 ::sprintf(errStr,
"%d", *seqIter);
2133 (std::string(
"Unknown edge index: ") + errStr)
2138 composite->m_geomVec.push_back(m_segGeoms[*seqIter]);
2145 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2148 if (m_triGeoms.count(*seqIter) == 0)
2150 char errStr[16] =
"";
2151 ::sprintf(errStr,
"%d", *seqIter);
2153 ErrorUtil::ewarning,
2154 (std::string(
"Unknown triangle index: ") + errStr)
2159 if (CheckRange(*m_triGeoms[*seqIter]))
2161 composite->m_geomVec.push_back(
2162 m_triGeoms[*seqIter]);
2171 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2174 if (m_quadGeoms.count(*seqIter) == 0)
2176 char errStr[16] =
"";
2177 ::sprintf(errStr,
"%d", *seqIter);
2179 (std::string(
"Unknown quad index: ") + errStr +
2180 std::string(
" in Composite section"))
2185 if (CheckRange(*m_quadGeoms[*seqIter]))
2187 composite->m_geomVec.push_back(
2188 m_quadGeoms[*seqIter]);
2196 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2199 if (*seqIter >= m_vertSet.size())
2201 char errStr[16] =
"";
2202 ::sprintf(errStr,
"%d", *seqIter);
2204 ErrorUtil::ewarning,
2205 (std::string(
"Unknown vertex index: ") + errStr)
2210 composite->m_geomVec.push_back(m_vertSet[*seqIter]);
2217 (std::string(
"Unrecognized composite token: ") + token)
2224 (std::string(
"Problem processing composite token: ") + token)
2231 void MeshGraphXml::ResolveGeomRef3D(
const std::string &prevToken,
2232 const std::string &token,
2237 std::istringstream tokenStream(token);
2238 std::istringstream prevTokenStream(prevToken);
2243 tokenStream >> type;
2245 std::string::size_type indxBeg = token.find_first_of(
'[') + 1;
2246 std::string::size_type indxEnd = token.find_last_of(
']') - 1;
2250 (std::string(
"Error reading index definition:") + token).c_str());
2252 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2254 std::vector<unsigned int> seqVector;
2255 std::vector<unsigned int>::iterator seqIter;
2257 bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
2260 (std::string(
"Error reading composite elements: ") + indxStr)
2263 prevTokenStream >> prevType;
2267 map<char, int> typeMap;
2278 bool validSequence =
2279 (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
2282 (std::string(
"Invalid combination of composite items: ") +
2283 type +
" and " + prevType +
".")
2289 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2292 if (m_vertSet.find(*seqIter) == m_vertSet.end())
2294 char errStr[16] =
"";
2295 ::sprintf(errStr,
"%d", *seqIter);
2297 ErrorUtil::ewarning,
2298 (std::string(
"Unknown vertex index: ") + errStr)
2303 composite->m_geomVec.push_back(m_vertSet[*seqIter]);
2309 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2312 if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
2314 char errStr[16] =
"";
2315 ::sprintf(errStr,
"%d", *seqIter);
2317 (std::string(
"Unknown edge index: ") + errStr)
2322 composite->m_geomVec.push_back(m_segGeoms[*seqIter]);
2328 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2334 char errStr[16] =
"";
2335 ::sprintf(errStr,
"%d", *seqIter);
2337 (std::string(
"Unknown face index: ") + errStr)
2342 if (CheckRange(*face))
2344 composite->m_geomVec.push_back(face);
2351 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2354 if (m_triGeoms.find(*seqIter) == m_triGeoms.end())
2356 char errStr[16] =
"";
2357 ::sprintf(errStr,
"%d", *seqIter);
2359 ErrorUtil::ewarning,
2360 (std::string(
"Unknown triangle index: ") + errStr)
2365 if (CheckRange(*m_triGeoms[*seqIter]))
2367 composite->m_geomVec.push_back(
2368 m_triGeoms[*seqIter]);
2375 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2378 if (m_quadGeoms.find(*seqIter) == m_quadGeoms.end())
2380 char errStr[16] =
"";
2381 ::sprintf(errStr,
"%d", *seqIter);
2383 (std::string(
"Unknown quad index: ") + errStr)
2388 if (CheckRange(*m_quadGeoms[*seqIter]))
2390 composite->m_geomVec.push_back(
2391 m_quadGeoms[*seqIter]);
2399 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2402 if (m_tetGeoms.find(*seqIter) == m_tetGeoms.end())
2404 char errStr[16] =
"";
2405 ::sprintf(errStr,
"%d", *seqIter);
2407 (std::string(
"Unknown tet index: ") + errStr)
2412 if (CheckRange(*m_tetGeoms[*seqIter]))
2414 composite->m_geomVec.push_back(
2415 m_tetGeoms[*seqIter]);
2423 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2426 if (m_pyrGeoms.find(*seqIter) == m_pyrGeoms.end())
2428 char errStr[16] =
"";
2429 ::sprintf(errStr,
"%d", *seqIter);
2431 ErrorUtil::ewarning,
2432 (std::string(
"Unknown pyramid index: ") + errStr)
2437 if (CheckRange(*m_pyrGeoms[*seqIter]))
2439 composite->m_geomVec.push_back(
2440 m_pyrGeoms[*seqIter]);
2448 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2451 if (m_prismGeoms.find(*seqIter) == m_prismGeoms.end())
2453 char errStr[16] =
"";
2454 ::sprintf(errStr,
"%d", *seqIter);
2456 (std::string(
"Unknown prism index: ") + errStr)
2461 if (CheckRange(*m_prismGeoms[*seqIter]))
2463 composite->m_geomVec.push_back(
2464 m_prismGeoms[*seqIter]);
2472 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2475 if (m_hexGeoms.find(*seqIter) == m_hexGeoms.end())
2477 char errStr[16] =
"";
2478 ::sprintf(errStr,
"%d", *seqIter);
2480 (std::string(
"Unknown hex index: ") + errStr)
2485 if (CheckRange(*m_hexGeoms[*seqIter]))
2487 composite->m_geomVec.push_back(
2488 m_hexGeoms[*seqIter]);
2496 (std::string(
"Unrecognized composite token: ") + token)
2503 (std::string(
"Problem processing composite token: ") + token)
2510 void MeshGraphXml::WriteVertices(TiXmlElement *geomTag,
PointGeomMap &verts)
2512 TiXmlElement *vertTag =
new TiXmlElement(
"VERTEX");
2514 for (
auto &i : verts)
2517 s << scientific << setprecision(8) << (*i.second)(0) <<
" "
2518 << (*i.second)(1) <<
" " << (*i.second)(2);
2519 TiXmlElement *v =
new TiXmlElement(
"V");
2520 v->SetAttribute(
"ID", i.second->GetGlobalID());
2521 v->LinkEndChild(
new TiXmlText(s.str()));
2522 vertTag->LinkEndChild(v);
2525 geomTag->LinkEndChild(vertTag);
2528 void MeshGraphXml::WriteEdges(TiXmlElement *geomTag,
SegGeomMap &edges)
2530 TiXmlElement *edgeTag =
2531 new TiXmlElement(m_meshDimension == 1 ?
"ELEMENT" :
"EDGE");
2532 string tag = m_meshDimension == 1 ?
"S" :
"E";
2534 for (
auto &i : edges)
2538 s << seg->GetVid(0) <<
" " << seg->GetVid(1);
2539 TiXmlElement *e =
new TiXmlElement(tag);
2540 e->SetAttribute(
"ID", i.first);
2541 e->LinkEndChild(
new TiXmlText(s.str()));
2542 edgeTag->LinkEndChild(e);
2545 geomTag->LinkEndChild(edgeTag);
2548 void MeshGraphXml::WriteTris(TiXmlElement *faceTag,
TriGeomMap &tris)
2552 for (
auto &i : tris)
2556 s << tri->GetEid(0) <<
" " << tri->GetEid(1) <<
" " << tri->GetEid(2);
2557 TiXmlElement *t =
new TiXmlElement(tag);
2558 t->SetAttribute(
"ID", i.first);
2559 t->LinkEndChild(
new TiXmlText(s.str()));
2560 faceTag->LinkEndChild(t);
2564 void MeshGraphXml::WriteQuads(TiXmlElement *faceTag,
QuadGeomMap &quads)
2568 for (
auto &i : quads)
2572 s << quad->GetEid(0) <<
" " << quad->GetEid(1) <<
" " << quad->GetEid(2)
2573 <<
" " << quad->GetEid(3);
2574 TiXmlElement *q =
new TiXmlElement(tag);
2575 q->SetAttribute(
"ID", i.first);
2576 q->LinkEndChild(
new TiXmlText(s.str()));
2577 faceTag->LinkEndChild(q);
2581 void MeshGraphXml::WriteHexs(TiXmlElement *elmtTag,
HexGeomMap &hexs)
2585 for (
auto &i : hexs)
2589 s << hex->GetFid(0) <<
" " << hex->GetFid(1) <<
" " << hex->GetFid(2)
2590 <<
" " << hex->GetFid(3) <<
" " << hex->GetFid(4) <<
" "
2591 << hex->GetFid(5) <<
" ";
2592 TiXmlElement *h =
new TiXmlElement(tag);
2593 h->SetAttribute(
"ID", i.first);
2594 h->LinkEndChild(
new TiXmlText(s.str()));
2595 elmtTag->LinkEndChild(h);
2603 for (
auto &i : pris)
2607 s << prism->GetFid(0) <<
" " << prism->GetFid(1) <<
" "
2608 << prism->GetFid(2) <<
" " << prism->GetFid(3) <<
" "
2609 << prism->GetFid(4) <<
" ";
2610 TiXmlElement *
p =
new TiXmlElement(tag);
2611 p->SetAttribute(
"ID", i.first);
2612 p->LinkEndChild(
new TiXmlText(s.str()));
2613 elmtTag->LinkEndChild(
p);
2617 void MeshGraphXml::WritePyrs(TiXmlElement *elmtTag,
PyrGeomMap &pyrs)
2621 for (
auto &i : pyrs)
2625 s << pyr->GetFid(0) <<
" " << pyr->GetFid(1) <<
" " << pyr->GetFid(2)
2626 <<
" " << pyr->GetFid(3) <<
" " << pyr->GetFid(4) <<
" ";
2627 TiXmlElement *
p =
new TiXmlElement(tag);
2628 p->SetAttribute(
"ID", i.first);
2629 p->LinkEndChild(
new TiXmlText(s.str()));
2630 elmtTag->LinkEndChild(
p);
2634 void MeshGraphXml::WriteTets(TiXmlElement *elmtTag,
TetGeomMap &tets)
2638 for (
auto &i : tets)
2642 s << tet->GetFid(0) <<
" " << tet->GetFid(1) <<
" " << tet->GetFid(2)
2643 <<
" " << tet->GetFid(3) <<
" ";
2644 TiXmlElement *t =
new TiXmlElement(tag);
2645 t->SetAttribute(
"ID", i.first);
2646 t->LinkEndChild(
new TiXmlText(s.str()));
2647 elmtTag->LinkEndChild(t);
2651 void MeshGraphXml::WriteCurves(TiXmlElement *geomTag,
CurveMap &edges,
2654 TiXmlElement *curveTag =
new TiXmlElement(
"CURVED");
2655 CurveMap::iterator curveIt;
2658 for (curveIt = edges.begin(); curveIt != edges.end(); ++curveIt)
2661 TiXmlElement *c =
new TiXmlElement(
"E");
2665 for (
int j = 0; j < curve->m_points.size(); ++j)
2668 s << scientific << (*p)(0) <<
" " << (*
p)(1) <<
" " << (*
p)(2)
2672 c->SetAttribute(
"ID", curveId++);
2673 c->SetAttribute(
"EDGEID", curve->m_curveID);
2674 c->SetAttribute(
"NUMPOINTS", curve->m_points.size());
2676 c->LinkEndChild(
new TiXmlText(s.str()));
2677 curveTag->LinkEndChild(c);
2680 for (curveIt = faces.begin(); curveIt != faces.end(); ++curveIt)
2683 TiXmlElement *c =
new TiXmlElement(
"F");
2687 for (
int j = 0; j < curve->m_points.size(); ++j)
2690 s << scientific << (*p)(0) <<
" " << (*
p)(1) <<
" " << (*
p)(2)
2694 c->SetAttribute(
"ID", curveId++);
2695 c->SetAttribute(
"FACEID", curve->m_curveID);
2696 c->SetAttribute(
"NUMPOINTS", curve->m_points.size());
2698 c->LinkEndChild(
new TiXmlText(s.str()));
2699 curveTag->LinkEndChild(c);
2702 geomTag->LinkEndChild(curveTag);
2705 void MeshGraphXml::WriteComposites(TiXmlElement *geomTag,
CompositeMap &comps)
2707 TiXmlElement *compTag =
new TiXmlElement(
"COMPOSITE");
2709 for (
auto &cIt : comps)
2711 if (cIt.second->m_geomVec.size() == 0)
2716 TiXmlElement *c =
new TiXmlElement(
"C");
2717 c->SetAttribute(
"ID", cIt.first);
2718 c->LinkEndChild(
new TiXmlText(GetCompositeString(cIt.second)));
2719 compTag->LinkEndChild(c);
2722 geomTag->LinkEndChild(compTag);
2725 void MeshGraphXml::WriteDomain(TiXmlElement *geomTag,
2726 map<int, CompositeMap> &domain)
2728 TiXmlElement *domTag =
new TiXmlElement(
"DOMAIN");
2730 for (
auto &d : domain)
2732 stringstream domString;
2733 if (d.second.size())
2736 TiXmlElement *dtag =
new TiXmlElement(
"D");
2737 dtag->SetAttribute(
"ID", d.first);
2739 vector<unsigned int> idxList;
2740 for (
auto cIt = d.second.begin(); cIt != d.second.end(); ++cIt)
2742 idxList.push_back(cIt->first);
2744 domString <<
" C[" << ParseUtils::GenerateSeqString(idxList)
2746 dtag->LinkEndChild(
new TiXmlText(domString.str()));
2747 domTag->LinkEndChild(dtag);
2751 geomTag->LinkEndChild(domTag);
2754 void MeshGraphXml::WriteDefaultExpansion(TiXmlElement *root)
2756 TiXmlElement *expTag =
new TiXmlElement(
"EXPANSIONS");
2758 for (
auto it = m_meshComposites.begin(); it != m_meshComposites.end(); it++)
2760 if (it->second->m_geomVec[0]->GetShapeDim() == m_meshDimension)
2762 TiXmlElement *exp =
new TiXmlElement(
"E");
2763 exp->SetAttribute(
"COMPOSITE",
2764 "C[" + boost::lexical_cast<string>(it->first) +
2766 exp->SetAttribute(
"NUMMODES", 4);
2767 exp->SetAttribute(
"TYPE",
"MODIFIED");
2768 exp->SetAttribute(
"FIELDS",
"u");
2770 expTag->LinkEndChild(exp);
2773 root->LinkEndChild(expTag);
2780 void MeshGraphXml::WriteGeometry(std::string &outfilename,
bool defaultExp,
2785 TiXmlDeclaration *decl =
new TiXmlDeclaration(
"1.0",
"utf-8",
"");
2786 doc.LinkEndChild(decl);
2788 TiXmlElement *root =
new TiXmlElement(
"NEKTAR");
2789 doc.LinkEndChild(root);
2790 TiXmlElement *geomTag =
new TiXmlElement(
"GEOMETRY");
2791 root->LinkEndChild(geomTag);
2799 geomTag->SetAttribute(
"DIM", m_meshDimension);
2800 geomTag->SetAttribute(
"SPACE", m_spaceDimension);
2806 WriteVertices(geomTag, m_vertSet);
2807 WriteEdges(geomTag, m_segGeoms);
2808 if (m_meshDimension > 1)
2810 TiXmlElement *faceTag =
2811 new TiXmlElement(m_meshDimension == 2 ?
"ELEMENT" :
"FACE");
2813 WriteTris(faceTag, m_triGeoms);
2814 WriteQuads(faceTag, m_quadGeoms);
2815 geomTag->LinkEndChild(faceTag);
2817 if (m_meshDimension > 2)
2819 TiXmlElement *elmtTag =
new TiXmlElement(
"ELEMENT");
2821 WriteHexs(elmtTag, m_hexGeoms);
2822 WritePyrs(elmtTag, m_pyrGeoms);
2823 WritePrisms(elmtTag, m_prismGeoms);
2824 WriteTets(elmtTag, m_tetGeoms);
2826 geomTag->LinkEndChild(elmtTag);
2828 WriteCurves(geomTag, m_curvedEdges, m_curvedFaces);
2829 WriteComposites(geomTag, m_meshComposites);
2830 WriteDomain(geomTag, m_domain);
2834 WriteDefaultExpansion(root);
2838 doc.SaveFile(outfilename);
2841 void MeshGraphXml::WriteXMLGeometry(std::string outname,
2842 vector<set<unsigned int>> elements,
2843 vector<unsigned int> partitions)
2854 string dirname = outname +
"_xml";
2855 boost::filesystem::path pdirname(dirname);
2857 if (!boost::filesystem::is_directory(dirname))
2859 boost::filesystem::create_directory(dirname);
2862 ASSERTL0(elements.size() == partitions.size(),
2863 "error in partitioned information");
2865 for (
int i = 0; i < partitions.size(); i++)
2868 TiXmlDeclaration *decl =
new TiXmlDeclaration(
"1.0",
"utf-8",
"");
2869 doc.LinkEndChild(decl);
2871 TiXmlElement *root = doc.FirstChildElement(
"NEKTAR");
2872 TiXmlElement *geomTag;
2877 root =
new TiXmlElement(
"NEKTAR");
2878 doc.LinkEndChild(root);
2880 geomTag =
new TiXmlElement(
"GEOMETRY");
2881 root->LinkEndChild(geomTag);
2886 geomTag = root->FirstChildElement(
"GEOMETRY");
2890 geomTag =
new TiXmlElement(
"GEOMETRY");
2891 root->LinkEndChild(geomTag);
2895 geomTag->SetAttribute(
"DIM", m_meshDimension);
2896 geomTag->SetAttribute(
"SPACE", m_spaceDimension);
2897 geomTag->SetAttribute(
"PARTITION", partitions[i]);
2912 vector<set<unsigned int>> entityIds(4);
2913 entityIds[m_meshDimension] = elements[i];
2915 switch (m_meshDimension)
2919 for (
auto &j : entityIds[3])
2922 if (m_hexGeoms.count(j))
2925 localHex[j] = m_hexGeoms[j];
2927 else if (m_pyrGeoms.count(j))
2930 localPyr[j] = m_pyrGeoms[j];
2932 else if (m_prismGeoms.count(j))
2934 g = m_prismGeoms[j];
2935 localPrism[j] = m_prismGeoms[j];
2937 else if (m_tetGeoms.count(j))
2940 localTet[j] = m_tetGeoms[j];
2944 ASSERTL0(
false,
"element in partition not found");
2947 for (
int k = 0; k < g->GetNumFaces(); k++)
2949 entityIds[2].insert(g->GetFid(k));
2951 for (
int k = 0; k < g->GetNumEdges(); k++)
2953 entityIds[1].insert(g->GetEid(k));
2955 for (
int k = 0; k < g->GetNumVerts(); k++)
2957 entityIds[0].insert(g->GetVid(k));
2964 for (
auto &j : entityIds[2])
2967 if (m_triGeoms.count(j))
2970 localTri[j] = m_triGeoms[j];
2972 else if (m_quadGeoms.count(j))
2975 localQuad[j] = m_quadGeoms[j];
2979 ASSERTL0(
false,
"element in partition not found");
2982 for (
int k = 0; k < g->GetNumEdges(); k++)
2984 entityIds[1].insert(g->GetEid(k));
2986 for (
int k = 0; k < g->GetNumVerts(); k++)
2988 entityIds[0].insert(g->GetVid(k));
2995 for (
auto &j : entityIds[1])
2998 if (m_segGeoms.count(j))
3001 localEdge[j] = m_segGeoms[j];
3005 ASSERTL0(
false,
"element in partition not found");
3008 for (
int k = 0; k < g->GetNumVerts(); k++)
3010 entityIds[0].insert(g->GetVid(k));
3016 if (m_meshDimension > 2)
3018 for (
auto &j : entityIds[2])
3020 if (m_triGeoms.count(j))
3022 localTri[j] = m_triGeoms[j];
3024 else if (m_quadGeoms.count(j))
3026 localQuad[j] = m_quadGeoms[j];
3035 if (m_meshDimension > 1)
3037 for (
auto &j : entityIds[1])
3039 if (m_segGeoms.count(j))
3041 localEdge[j] = m_segGeoms[j];
3050 for (
auto &j : entityIds[0])
3052 if (m_vertSet.count(j))
3054 localVert[j] = m_vertSet[j];
3062 WriteVertices(geomTag, localVert);
3063 WriteEdges(geomTag, localEdge);
3064 if (m_meshDimension > 1)
3066 TiXmlElement *faceTag =
3067 new TiXmlElement(m_meshDimension == 2 ?
"ELEMENT" :
"FACE");
3069 WriteTris(faceTag, localTri);
3070 WriteQuads(faceTag, localQuad);
3071 geomTag->LinkEndChild(faceTag);
3073 if (m_meshDimension > 2)
3075 TiXmlElement *elmtTag =
new TiXmlElement(
"ELEMENT");
3077 WriteHexs(elmtTag, localHex);
3078 WritePyrs(elmtTag, localPyr);
3079 WritePrisms(elmtTag, localPrism);
3080 WriteTets(elmtTag, localTet);
3082 geomTag->LinkEndChild(elmtTag);
3085 for (
auto &j : localTri)
3087 if (m_curvedFaces.count(j.first))
3089 localCurveFace[j.first] = m_curvedFaces[j.first];
3092 for (
auto &j : localQuad)
3094 if (m_curvedFaces.count(j.first))
3096 localCurveFace[j.first] = m_curvedFaces[j.first];
3099 for (
auto &j : localEdge)
3101 if (m_curvedEdges.count(j.first))
3103 localCurveEdge[j.first] = m_curvedEdges[j.first];
3107 WriteCurves(geomTag, localCurveEdge, localCurveFace);
3111 for (
auto &j : m_meshComposites)
3114 int dim = j.second->m_geomVec[0]->GetShapeDim();
3116 for (
int k = 0; k < j.second->m_geomVec.size(); k++)
3118 if (entityIds[dim].count(j.second->m_geomVec[k]->GetGlobalID()))
3120 comp->m_geomVec.push_back(j.second->m_geomVec[k]);
3124 if (comp->m_geomVec.size())
3126 localComp[j.first] = comp;
3130 WriteComposites(geomTag, localComp);
3132 map<int, CompositeMap> domain;
3133 for (
auto &d : m_domain)
3136 for (
auto &j : localComp)
3138 if (d.second.count(j.first))
3140 domMap[j.first] = j.second;
3143 domain[d.first] = domMap;
3146 WriteDomain(geomTag, domain);
3148 if (m_session->DefinesElement(
"NEKTAR/CONDITIONS"))
3150 std::set<int> vBndRegionIdList;
3151 TiXmlElement *vConditions =
3152 new TiXmlElement(*m_session->GetElement(
"Nektar/Conditions"));
3153 TiXmlElement *vBndRegions =
3154 vConditions->FirstChildElement(
"BOUNDARYREGIONS");
3155 TiXmlElement *vBndConditions =
3156 vConditions->FirstChildElement(
"BOUNDARYCONDITIONS");
3157 TiXmlElement *vItem;
3161 TiXmlElement *vNewBndRegions =
3162 new TiXmlElement(
"BOUNDARYREGIONS");
3163 vItem = vBndRegions->FirstChildElement();
3166 std::string vSeqStr =
3167 vItem->FirstChild()->ToText()->Value();
3168 std::string::size_type indxBeg =
3169 vSeqStr.find_first_of(
'[') + 1;
3170 std::string::size_type indxEnd =
3171 vSeqStr.find_last_of(
']') - 1;
3172 vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
3173 std::vector<unsigned int> vSeq;
3174 ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
3176 vector<unsigned int> idxList;
3178 for (
unsigned int i = 0; i < vSeq.size(); ++i)
3180 if (localComp.find(vSeq[i]) != localComp.end())
3182 idxList.push_back(vSeq[i]);
3185 int p = atoi(vItem->Attribute(
"ID"));
3187 std::string vListStr =
3188 ParseUtils::GenerateSeqString(idxList);
3190 if (vListStr.length() == 0)
3192 TiXmlElement *tmp = vItem;
3193 vItem = vItem->NextSiblingElement();
3194 vBndRegions->RemoveChild(tmp);
3198 vListStr =
"C[" + vListStr +
"]";
3199 TiXmlText *vList =
new TiXmlText(vListStr);
3200 TiXmlElement *vNewElement =
new TiXmlElement(
"B");
3201 vNewElement->SetAttribute(
"ID",
p);
3202 vNewElement->LinkEndChild(vList);
3203 vNewBndRegions->LinkEndChild(vNewElement);
3204 vBndRegionIdList.insert(
p);
3205 vItem = vItem->NextSiblingElement();
3209 m_bndRegOrder[
p] = vSeq;
3211 vConditions->ReplaceChild(vBndRegions, *vNewBndRegions);
3216 vItem = vBndConditions->FirstChildElement();
3219 std::set<int>::iterator x;
3220 if ((x = vBndRegionIdList.find(atoi(vItem->Attribute(
3221 "REF")))) != vBndRegionIdList.end())
3223 vItem->SetAttribute(
"REF", *x);
3224 vItem = vItem->NextSiblingElement();
3228 TiXmlElement *tmp = vItem;
3229 vItem = vItem->NextSiblingElement();
3230 vBndConditions->RemoveChild(tmp);
3234 root->LinkEndChild(vConditions);
3238 TiXmlElement *vSrc =
3239 m_session->GetElement(
"Nektar")->FirstChildElement();
3242 std::string vName = boost::to_upper_copy(vSrc->ValueStr());
3243 if (vName !=
"GEOMETRY" && vName !=
"CONDITIONS")
3245 root->LinkEndChild(
new TiXmlElement(*vSrc));
3247 vSrc = vSrc->NextSiblingElement();
3253 pad % partitions[i];
3254 boost::filesystem::path pFilename(pad.str());
3256 boost::filesystem::path fullpath = pdirname / pFilename;
3265 for (
auto &c : m_meshComposites)
3267 bool fillComp =
true;
3268 for (
auto &d : m_domain[0])
3270 if (c.second == d.second)
3277 std::vector<unsigned int> ids;
3278 for (
auto &elmt : c.second->m_geomVec)
3280 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...
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.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::map< std::string, std::string > FieldMetaDataMap
const std::string kPointsTypeStr[]
std::string PortablePath(const boost::filesystem::path &path)
create portable path on different platforms for boost::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::map< int, SegGeomSharedPtr > SegGeomMap
std::shared_ptr< Curve > CurveSharedPtr
std::unordered_map< int, CurveSharedPtr > CurveMap
std::map< int, CompositeSharedPtr > CompositeMap
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::map< int, TetGeomSharedPtr > TetGeomMap
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
std::map< int, PrismGeomSharedPtr > PrismGeomMap
std::shared_ptr< HexGeom > HexGeomSharedPtr
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
MeshGraphFactory & GetMeshGraphFactory()
std::shared_ptr< Composite > CompositeSharedPtr
std::map< int, HexGeomSharedPtr > HexGeomMap
std::map< int, PointGeomSharedPtr > PointGeomMap
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
The above copyright notice and this permission notice shall be included.