Nektar++
MeshGraphIOXml.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: MeshGraphIOXml.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33//
34////////////////////////////////////////////////////////////////////////////////
35
36#include <iomanip>
37
41
46
47#include <boost/format.hpp>
48
49#include <tinyxml.h>
50
52{
53
54std::string MeshGraphIOXml::className =
56 "Xml", MeshGraphIOXml::create, "IO with Xml geometry");
57
60{
61 // Get row of comm, or the whole comm if not split
62 LibUtilities::CommSharedPtr comm = session->GetComm();
63 LibUtilities::CommSharedPtr commMesh = comm->GetRowComm();
64 const bool isRoot = comm->TreatAsRankZero();
65
66 m_session = session;
67
68 int meshDimension;
69
70 // Load file for root process only (since this is always needed)
71 // and determine if the provided geometry has already been
72 // partitioned. This will be the case if the user provides the
73 // directory of mesh partitions as an input. Partitioned geometries
74 // have the attribute
75 // PARTITION=X
76 // where X is the number of the partition (and should match the
77 // process rank). The result is shared with all other processes.
78 int isPartitioned = 0;
79 if (isRoot)
80 {
81 if (m_session->DefinesElement("Nektar/Geometry"))
82 {
83 if (m_session->GetElement("Nektar/Geometry")
84 ->Attribute("PARTITION"))
85 {
86 std::cout << "Using pre-partitioned mesh." << std::endl;
87 isPartitioned = 1;
88 }
89 m_session->GetElement("NEKTAR/GEOMETRY")
90 ->QueryIntAttribute("DIM", &meshDimension);
91 }
92 }
93 comm->Bcast(isPartitioned, 0);
94 comm->Bcast(meshDimension, 0);
95
96 // If the mesh is already partitioned, we are done. Remaining
97 // processes must load their partitions.
98 if (isPartitioned)
99 {
100 if (!isRoot)
101 {
102 m_session->InitSession();
103 }
104 }
105 else
106 {
107 // Default partitioner to use is Scotch, if it is installed.
108 bool haveScotch = GetMeshPartitionFactory().ModuleExists("Scotch");
109 bool haveMetis = GetMeshPartitionFactory().ModuleExists("Metis");
110
111 std::string partitionerName = "Scotch";
112 if (!haveScotch && haveMetis)
113 {
114 partitionerName = "Metis";
115 }
116
117 // Override default with command-line flags if they are set.
118 if (session->DefinesCmdLineArgument("use-metis"))
119 {
120 partitionerName = "Metis";
121 }
122 if (session->DefinesCmdLineArgument("use-scotch"))
123 {
124 partitionerName = "Scotch";
125 }
126
127 // Mesh has not been partitioned so do partitioning if required. Note
128 // in the serial case nothing is done as we have already loaded the
129 // mesh.
130 if (session->DefinesCmdLineArgument("part-only") ||
131 session->DefinesCmdLineArgument("part-only-overlapping"))
132 {
133 // Perform partitioning of the mesh only. For this we insist the
134 // code is run in serial (parallel execution is pointless).
135 ASSERTL0(comm->GetSize() == 1,
136 "The 'part-only' option should be used in serial.");
137
138 // Read 'lite' geometry information
140
141 // Number of partitions is specified by the parameter.
142 int nParts;
143 auto comp = CreateCompositeDescriptor();
144
145 MeshPartitionSharedPtr partitioner =
147 partitionerName, session, session->GetComm(), meshDimension,
148 CreateMeshEntities(), comp);
149
150 if (session->DefinesCmdLineArgument("part-only"))
151 {
152 nParts = session->GetCmdLineArgument<int>("part-only");
153 partitioner->PartitionMesh(nParts, true);
154 }
155 else
156 {
157 nParts =
158 session->GetCmdLineArgument<int>("part-only-overlapping");
159 partitioner->PartitionMesh(nParts, true, true);
160 }
161
162 std::vector<std::set<unsigned int>> elmtIDs;
163 std::vector<unsigned int> parts(nParts);
164 for (int i = 0; i < nParts; ++i)
165 {
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);
171 parts[i] = i;
172 }
173
174 this->WriteXMLGeometry(m_session->GetSessionName(), elmtIDs, parts);
175
176 if (isRoot && session->DefinesCmdLineArgument("part-info"))
177 {
178 partitioner->PrintPartInfo(std::cout);
179 }
180
181 session->Finalise();
182 exit(0);
183 }
184
185 if (commMesh->GetSize() > 1)
186 {
187 ASSERTL0(haveScotch || haveMetis,
188 "Valid partitioner not found! Either Scotch or METIS "
189 "should be used.");
190
191 int nParts = commMesh->GetSize();
192
193 if (session->GetSharedFilesystem())
194 {
195 std::vector<unsigned int> keys, vals;
196 int i;
197
198 if (isRoot)
199 {
200 // Read 'lite' geometry information
202
203 // Store composite ordering and boundary information.
205 m_meshGraph->SetCompositeOrdering(m_compOrder);
206 auto comp = CreateCompositeDescriptor();
207
208 // Create mesh partitioner.
209 MeshPartitionSharedPtr partitioner =
211 partitionerName, session, session->GetComm(),
212 meshDimension, CreateMeshEntities(), comp);
213
214 partitioner->PartitionMesh(nParts, true);
215
216 std::vector<std::set<unsigned int>> elmtIDs;
217 std::vector<unsigned int> parts(nParts);
218 for (i = 0; i < nParts; ++i)
219 {
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);
225 parts[i] = i;
226 }
227
228 // Call WriteGeometry to write out partition files. This
229 // will populate m_bndRegOrder.
230 this->WriteXMLGeometry(m_session->GetSessionName(), elmtIDs,
231 parts);
232
233 // Communicate orderings to the other processors.
234
235 // First send sizes of the orderings and boundary
236 // regions to allocate storage on the remote end.
237 keys.resize(2);
238 keys[0] = m_compOrder.size();
239 keys[1] = m_bndRegOrder.size();
240 comm->Bcast(keys, 0);
241
242 // Construct the keys and sizes of values for composite
243 // ordering
244 keys.resize(m_compOrder.size());
245 vals.resize(m_compOrder.size());
246
247 i = 0;
248 for (auto &cIt : m_compOrder)
249 {
250 keys[i] = cIt.first;
251 vals[i++] = cIt.second.size();
252 }
253
254 // Send across data.
255 comm->Bcast(keys, 0);
256 comm->Bcast(vals, 0);
257 for (auto &cIt : m_compOrder)
258 {
259 comm->Bcast(cIt.second, 0);
260 }
261 m_meshGraph->SetCompositeOrdering(m_compOrder);
262
263 // Construct the keys and sizes of values for composite
264 // ordering
265 keys.resize(m_bndRegOrder.size());
266 vals.resize(m_bndRegOrder.size());
267
268 i = 0;
269 for (auto &bIt : m_bndRegOrder)
270 {
271 keys[i] = bIt.first;
272 vals[i++] = bIt.second.size();
273 }
274
275 // Send across data.
276 if (!keys.empty())
277 {
278 comm->Bcast(keys, 0);
279 }
280 if (!vals.empty())
281 {
282 comm->Bcast(vals, 0);
283 }
284 for (auto &bIt : m_bndRegOrder)
285 {
286 comm->Bcast(bIt.second, 0);
287 }
288 m_meshGraph->SetBndRegionOrdering(m_bndRegOrder);
289
290 if (session->DefinesCmdLineArgument("part-info"))
291 {
292 partitioner->PrintPartInfo(std::cout);
293 }
294 }
295 else
296 {
297 keys.resize(2);
298 comm->Bcast(keys, 0);
299
300 int cmpSize = keys[0];
301 int bndSize = keys[1];
302
303 keys.resize(cmpSize);
304 vals.resize(cmpSize);
305 comm->Bcast(keys, 0);
306 comm->Bcast(vals, 0);
307
308 for (int i = 0; i < keys.size(); ++i)
309 {
310 std::vector<unsigned int> tmp(vals[i]);
311 comm->Bcast(tmp, 0);
312 m_compOrder[keys[i]] = tmp;
313 }
314 m_meshGraph->SetCompositeOrdering(m_compOrder);
315
316 keys.resize(bndSize);
317 vals.resize(bndSize);
318 if (!keys.empty())
319 {
320 comm->Bcast(keys, 0);
321 }
322 if (!vals.empty())
323 {
324 comm->Bcast(vals, 0);
325 }
326 for (int i = 0; i < keys.size(); ++i)
327 {
328 std::vector<unsigned int> tmp(vals[i]);
329 comm->Bcast(tmp, 0);
330 m_bndRegOrder[keys[i]] = tmp;
331 }
332 m_meshGraph->SetBndRegionOrdering(m_bndRegOrder);
333 }
334 }
335 else
336 {
337 m_session->InitSession();
339
341 m_meshGraph->SetCompositeOrdering(m_compOrder);
342 auto comp = CreateCompositeDescriptor();
343
344 // Partitioner now operates in parallel. Each process receives
345 // partitioning over interconnect and writes its own session
346 // file to the working directory.
347 MeshPartitionSharedPtr partitioner =
349 partitionerName, session, session->GetComm(),
350 meshDimension, CreateMeshEntities(), comp);
351
352 partitioner->PartitionMesh(nParts, false);
353
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());
359
360 // if (comm->GetTimeComm()->GetRank() == 0) // FIXME
361 // (OpenMPI 3.1.3)
362 {
363 this->WriteXMLGeometry(session->GetSessionName(), elIDs,
364 parts);
365 }
366
367 if (m_session->DefinesCmdLineArgument("part-info") && isRoot)
368 {
369 partitioner->PrintPartInfo(std::cout);
370 }
371 }
372
373 // Wait for all processors to finish their writing activities.
374 comm->Block();
375
376 std::string dirname = m_session->GetSessionName() + "_xml";
377 fs::path pdirname(dirname);
378 boost::format pad("P%1$07d.xml");
379 pad % comm->GetRowComm()->GetRank();
380 fs::path pFilename(pad.str());
381 fs::path fullpath = pdirname / pFilename;
382
383 std::vector<std::string> filenames = {
385 m_session->InitSession(filenames);
386 }
387 else if (!isRoot)
388 {
389 // No partitioning, non-root processors need to read the session
390 // file -- handles case where --npz is the same as number of
391 // processors.
392 m_session->InitSession();
393 }
394 }
395}
396
398 bool fillGraph)
399{
400 // Reset member variables.
401 m_meshGraph->Clear();
402
403 m_meshGraph->SetDomainRange(rng);
404 m_xmlGeom = m_session->GetElement("NEKTAR/GEOMETRY");
405
406 int err; /// Error value returned by TinyXML.
407
408 TiXmlAttribute *attr = m_xmlGeom->FirstAttribute();
409
410 // Initialize the mesh and space dimensions to 3 dimensions.
411 // We want to do this each time we read a file, so it should
412 // be done here and not just during class initialization.
413 int meshDimension = m_meshGraph->GetMeshDimension();
414 int spaceDimension = m_meshGraph->GetSpaceDimension();
415 m_meshPartitioned = false;
416 meshDimension = 3;
417 spaceDimension = 3;
418
419 while (attr)
420 {
421 std::string attrName(attr->Name());
422 if (attrName == "DIM")
423 {
424 err = attr->QueryIntValue(&meshDimension);
425 ASSERTL0(err == TIXML_SUCCESS, "Unable to read mesh dimension.");
426 }
427 else if (attrName == "SPACE")
428 {
429 err = attr->QueryIntValue(&spaceDimension);
430 ASSERTL0(err == TIXML_SUCCESS, "Unable to read space dimension.");
431 }
432 else if (attrName == "PARTITION")
433 {
434 err = attr->QueryIntValue(&m_partition);
435 ASSERTL0(err == TIXML_SUCCESS, "Unable to read partition.");
436 m_meshPartitioned = true;
437 m_meshGraph->SetMeshPartitioned(true);
438 }
439 else
440 {
441 std::string errstr("Unknown attribute: ");
442 errstr += attrName;
443 ASSERTL0(false, errstr.c_str());
444 }
445
446 // Get the next attribute.
447 attr = attr->Next();
448 }
449
450 m_meshGraph->SetMeshDimension(meshDimension);
451 m_meshGraph->SetSpaceDimension(spaceDimension);
452
453 ASSERTL0(meshDimension <= spaceDimension,
454 "Mesh dimension greater than space dimension");
455
457 v_ReadCurves();
458 if (meshDimension >= 2)
459 {
460 v_ReadEdges();
461 if (meshDimension == 3)
462 {
463 v_ReadFaces();
464 }
465 }
466 ReadElements();
468 ReadDomain();
469
470 if (fillGraph)
471 {
472 m_meshGraph->FillGraph();
473 }
474}
475
477{
478 auto &vertSet = m_meshGraph->GetAllPointGeoms();
479 int spaceDimension = m_meshGraph->GetSpaceDimension();
480
481 // Now read the vertices
482 TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
483 ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
484
485 NekDouble xscale, yscale, zscale;
486
487 // check to see if any scaling parameters are in
488 // attributes and determine these values
489 LibUtilities::Interpreter expEvaluator;
490 const char *xscal = element->Attribute("XSCALE");
491 if (!xscal)
492 {
493 xscale = 1.0;
494 }
495 else
496 {
497 std::string xscalstr = xscal;
498 int expr_id = expEvaluator.DefineFunction("", xscalstr);
499 xscale = expEvaluator.Evaluate(expr_id);
500 }
501
502 const char *yscal = element->Attribute("YSCALE");
503 if (!yscal)
504 {
505 yscale = 1.0;
506 }
507 else
508 {
509 std::string yscalstr = yscal;
510 int expr_id = expEvaluator.DefineFunction("", yscalstr);
511 yscale = expEvaluator.Evaluate(expr_id);
512 }
513
514 const char *zscal = element->Attribute("ZSCALE");
515 if (!zscal)
516 {
517 zscale = 1.0;
518 }
519 else
520 {
521 std::string zscalstr = zscal;
522 int expr_id = expEvaluator.DefineFunction("", zscalstr);
523 zscale = expEvaluator.Evaluate(expr_id);
524 }
525
526 NekDouble xmove, ymove, zmove;
527
528 // check to see if any moving parameters are in
529 // attributes and determine these values
530
531 const char *xmov = element->Attribute("XMOVE");
532 if (!xmov)
533 {
534 xmove = 0.0;
535 }
536 else
537 {
538 std::string xmovstr = xmov;
539 int expr_id = expEvaluator.DefineFunction("", xmovstr);
540 xmove = expEvaluator.Evaluate(expr_id);
541 }
542
543 const char *ymov = element->Attribute("YMOVE");
544 if (!ymov)
545 {
546 ymove = 0.0;
547 }
548 else
549 {
550 std::string ymovstr = ymov;
551 int expr_id = expEvaluator.DefineFunction("", ymovstr);
552 ymove = expEvaluator.Evaluate(expr_id);
553 }
554
555 const char *zmov = element->Attribute("ZMOVE");
556 if (!zmov)
557 {
558 zmove = 0.0;
559 }
560 else
561 {
562 std::string zmovstr = zmov;
563 int expr_id = expEvaluator.DefineFunction("", zmovstr);
564 zmove = expEvaluator.Evaluate(expr_id);
565 }
566
567 TiXmlElement *vertex = element->FirstChildElement("V");
568
569 int indx;
570
571 while (vertex)
572 {
573 TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
574 std::string attrName(vertexAttr->Name());
575
576 ASSERTL0(attrName == "ID",
577 (std::string("Unknown attribute name: ") + attrName).c_str());
578
579 int err = vertexAttr->QueryIntValue(&indx);
580 ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
581
582 // Now read body of vertex
583 std::string vertexBodyStr;
584
585 TiXmlNode *vertexBody = vertex->FirstChild();
586
587 while (vertexBody)
588 {
589 // Accumulate all non-comment body data.
590 if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
591 {
592 vertexBodyStr += vertexBody->ToText()->Value();
593 vertexBodyStr += " ";
594 }
595
596 vertexBody = vertexBody->NextSibling();
597 }
598
599 ASSERTL0(!vertexBodyStr.empty(),
600 "Vertex definitions must contain vertex data.");
601
602 // Get vertex data from the data string.
603 NekDouble xval, yval, zval;
604 std::istringstream vertexDataStrm(vertexBodyStr.c_str());
605
606 try
607 {
608 while (!vertexDataStrm.fail())
609 {
610 vertexDataStrm >> xval >> yval >> zval;
611
612 xval = xval * xscale + xmove;
613 yval = yval * yscale + ymove;
614 zval = zval * zscale + zmove;
615
616 // Need to check it here because we may not be
617 // good after the read indicating that there
618 // was nothing to read.
619 if (!vertexDataStrm.fail())
620 {
623 spaceDimension, indx, xval, yval, zval));
624 vertSet[indx] = vert;
625 }
626 }
627 }
628 catch (...)
629 {
630 ASSERTL0(false, "Unable to read VERTEX data.");
631 }
632
633 vertex = vertex->NextSiblingElement("V");
634 }
635}
636
638{
639 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
640 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
641 int meshDimension = m_meshGraph->GetMeshDimension();
642
643 // check to see if any scaling parameters are in
644 // attributes and determine these values
645 TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
646 ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
647
648 NekDouble xscale, yscale, zscale;
649
650 LibUtilities::Interpreter expEvaluator;
651 const char *xscal = element->Attribute("XSCALE");
652 if (!xscal)
653 {
654 xscale = 1.0;
655 }
656 else
657 {
658 std::string xscalstr = xscal;
659 int expr_id = expEvaluator.DefineFunction("", xscalstr);
660 xscale = expEvaluator.Evaluate(expr_id);
661 }
662
663 const char *yscal = element->Attribute("YSCALE");
664 if (!yscal)
665 {
666 yscale = 1.0;
667 }
668 else
669 {
670 std::string yscalstr = yscal;
671 int expr_id = expEvaluator.DefineFunction("", yscalstr);
672 yscale = expEvaluator.Evaluate(expr_id);
673 }
674
675 const char *zscal = element->Attribute("ZSCALE");
676 if (!zscal)
677 {
678 zscale = 1.0;
679 }
680 else
681 {
682 std::string zscalstr = zscal;
683 int expr_id = expEvaluator.DefineFunction("", zscalstr);
684 zscale = expEvaluator.Evaluate(expr_id);
685 }
686
687 NekDouble xmove, ymove, zmove;
688
689 // check to see if any moving parameters are in
690 // attributes and determine these values
691
692 const char *xmov = element->Attribute("XMOVE");
693 if (!xmov)
694 {
695 xmove = 0.0;
696 }
697 else
698 {
699 std::string xmovstr = xmov;
700 int expr_id = expEvaluator.DefineFunction("", xmovstr);
701 xmove = expEvaluator.Evaluate(expr_id);
702 }
703
704 const char *ymov = element->Attribute("YMOVE");
705 if (!ymov)
706 {
707 ymove = 0.0;
708 }
709 else
710 {
711 std::string ymovstr = ymov;
712 int expr_id = expEvaluator.DefineFunction("", ymovstr);
713 ymove = expEvaluator.Evaluate(expr_id);
714 }
715
716 const char *zmov = element->Attribute("ZMOVE");
717 if (!zmov)
718 {
719 zmove = 0.0;
720 }
721 else
722 {
723 std::string zmovstr = zmov;
724 int expr_id = expEvaluator.DefineFunction("", zmovstr);
725 zmove = expEvaluator.Evaluate(expr_id);
726 }
727
728 int err;
729
730 /// Look for elements in CURVE block.
731 TiXmlElement *field = m_xmlGeom->FirstChildElement("CURVED");
732
733 if (!field) // return if no curved entities
734 {
735 return;
736 }
737
738 /// All curves are of the form: "<? ID="#" TYPE="GLL OR other
739 /// points type" NUMPOINTS="#"> ... </?>", with ? being an
740 /// element type (either E or F).
741
742 TiXmlElement *edgelement = field->FirstChildElement("E");
743
744 int edgeindx, edgeid;
745
746 while (edgelement)
747 {
748 std::string edge(edgelement->ValueStr());
749 ASSERTL0(edge == "E",
750 (std::string("Unknown 3D curve type:") + edge).c_str());
751
752 /// Read id attribute.
753 err = edgelement->QueryIntAttribute("ID", &edgeindx);
754 ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute ID.");
755
756 /// Read edge id attribute.
757 err = edgelement->QueryIntAttribute("EDGEID", &edgeid);
758 ASSERTL0(err == TIXML_SUCCESS,
759 "Unable to read curve attribute EDGEID.");
760
761 /// Read text edgelement description.
762 std::string elementStr;
763 TiXmlNode *elementChild = edgelement->FirstChild();
764
765 while (elementChild)
766 {
767 // Accumulate all non-comment element data
768 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
769 {
770 elementStr += elementChild->ToText()->ValueStr();
771 elementStr += " ";
772 }
773 elementChild = elementChild->NextSibling();
774 }
775
776 ASSERTL0(!elementStr.empty(), "Unable to read curve description body.");
777
778 /// Parse out the element components corresponding to type of
779 /// element.
780 if (edge == "E")
781 {
782 int numPts = 0;
783 // Determine the points type
784 std::string typeStr = edgelement->Attribute("TYPE");
785 ASSERTL0(!typeStr.empty(), "TYPE must be specified in "
786 "points definition");
787
789 const std::string *begStr = LibUtilities::kPointsTypeStr;
790 const std::string *endStr =
792 const std::string *ptsStr = std::find(begStr, endStr, typeStr);
793
794 ASSERTL0(ptsStr != endStr, "Invalid points type.");
795 type = (LibUtilities::PointsType)(ptsStr - begStr);
796
797 // Determine the number of points
798 err = edgelement->QueryIntAttribute("NUMPOINTS", &numPts);
799 ASSERTL0(err == TIXML_SUCCESS,
800 "Unable to read curve attribute NUMPOINTS.");
801 CurveSharedPtr curve(
803
804 // Read points (x, y, z)
805 NekDouble xval, yval, zval;
806 std::istringstream elementDataStrm(elementStr.c_str());
807 try
808 {
809 while (!elementDataStrm.fail())
810 {
811 elementDataStrm >> xval >> yval >> zval;
812
813 xval = xval * xscale + xmove;
814 yval = yval * yscale + ymove;
815 zval = zval * zscale + zmove;
816
817 // Need to check it here because we may not be
818 // good after the read indicating that there
819 // was nothing to read.
820 if (!elementDataStrm.fail())
821 {
824 meshDimension, edgeindx, xval, yval, zval));
825
826 curve->m_points.push_back(vert);
827 }
828 }
829 }
830 catch (...)
831 {
833 (std::string("Unable to read curve data for EDGE: ") +
834 elementStr)
835 .c_str());
836 }
837
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));
843
844 curvedEdges[edgeid] = curve;
845
846 edgelement = edgelement->NextSiblingElement("E");
847
848 } // end if-loop
849
850 } // end while-loop
851
852 TiXmlElement *facelement = field->FirstChildElement("F");
853 int faceindx, faceid;
854
855 while (facelement)
856 {
857 std::string face(facelement->ValueStr());
858 ASSERTL0(face == "F",
859 (std::string("Unknown 3D curve type: ") + face).c_str());
860
861 /// Read id attribute.
862 err = facelement->QueryIntAttribute("ID", &faceindx);
863 ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute ID.");
864
865 /// Read face id attribute.
866 err = facelement->QueryIntAttribute("FACEID", &faceid);
867 ASSERTL0(err == TIXML_SUCCESS,
868 "Unable to read curve attribute FACEID.");
869
870 /// Read text face element description.
871 std::string elementStr;
872 TiXmlNode *elementChild = facelement->FirstChild();
873
874 while (elementChild)
875 {
876 // Accumulate all non-comment element data
877 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
878 {
879 elementStr += elementChild->ToText()->ValueStr();
880 elementStr += " ";
881 }
882 elementChild = elementChild->NextSibling();
883 }
884
885 ASSERTL0(!elementStr.empty(), "Unable to read curve description body.");
886
887 /// Parse out the element components corresponding to type of
888 /// element.
889 if (face == "F")
890 {
891 std::string typeStr = facelement->Attribute("TYPE");
892 ASSERTL0(!typeStr.empty(), "TYPE must be specified in "
893 "points definition");
895 const std::string *begStr = LibUtilities::kPointsTypeStr;
896 const std::string *endStr =
898 const std::string *ptsStr = std::find(begStr, endStr, typeStr);
899
900 ASSERTL0(ptsStr != endStr, "Invalid points type.");
901 type = (LibUtilities::PointsType)(ptsStr - begStr);
902
903 std::string numptsStr = facelement->Attribute("NUMPOINTS");
904 ASSERTL0(!numptsStr.empty(),
905 "NUMPOINTS must be specified in points definition");
906 int numPts = 0;
907 std::stringstream s;
908 s << numptsStr;
909 s >> numPts;
910
911 CurveSharedPtr curve(
913
914 ASSERTL0(numPts >= 3, "NUMPOINTS for face must be greater than 2");
915
916 if (numPts == 3)
917 {
918 ASSERTL0(ptsStr != endStr, "Invalid points type.");
919 }
920
921 // Read points (x, y, z)
922 NekDouble xval, yval, zval;
923 std::istringstream elementDataStrm(elementStr.c_str());
924 try
925 {
926 while (!elementDataStrm.fail())
927 {
928 elementDataStrm >> xval >> yval >> zval;
929
930 // Need to check it here because we
931 // may not be good after the read
932 // indicating that there was nothing
933 // to read.
934 if (!elementDataStrm.fail())
935 {
938 meshDimension, faceindx, xval, yval, zval));
939 curve->m_points.push_back(vert);
940 }
941 }
942 }
943 catch (...)
944 {
946 (std::string("Unable to read curve data for FACE: ") +
947 elementStr)
948 .c_str());
949 }
950 curvedFaces[faceid] = curve;
951
952 facelement = facelement->NextSiblingElement("F");
953 }
954 }
955}
956
958{
959 auto &domain = m_meshGraph->GetDomain();
960
961 TiXmlElement *element = nullptr;
962 /// Look for data in DOMAIN block.
963 element = m_xmlGeom->FirstChildElement("DOMAIN");
964
965 ASSERTL0(element, "Unable to find DOMAIN tag in file.");
966
967 /// Elements are of the form: "<D ID = "N"> ... </D>".
968 /// Read the ID field first.
969 TiXmlElement *multidomains = element->FirstChildElement("D");
970
971 if (multidomains)
972 {
973 while (multidomains)
974 {
975 int indx;
976 int err = multidomains->QueryIntAttribute("ID", &indx);
977 ASSERTL0(err == TIXML_SUCCESS,
978 "Unable to read attribute ID in Domain.");
979
980 TiXmlNode *elementChild = multidomains->FirstChild();
981 while (elementChild &&
982 elementChild->Type() != TiXmlNode::TINYXML_TEXT)
983 {
984 elementChild = elementChild->NextSibling();
985 }
986
987 ASSERTL0(elementChild, "Unable to read DOMAIN body.");
988 std::string elementStr = elementChild->ToText()->ValueStr();
989
990 elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
991
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);
996
997 ASSERTL0(
998 !indxStr.empty(),
999 "Unable to read domain's composite index (index missing?).");
1000
1001 // Read the domain composites.
1002 // Parse the composites into a list.
1003 std::map<int, CompositeSharedPtr> unrollDomain;
1004 m_meshGraph->GetCompositeList(indxStr, unrollDomain);
1005 domain[indx] = unrollDomain;
1006
1007 ASSERTL0(!domain[indx].empty(),
1008 (std::string(
1009 "Unable to obtain domain's referenced composite: ") +
1010 indxStr)
1011 .c_str());
1012
1013 /// Keep looking
1014 multidomains = multidomains->NextSiblingElement("D");
1015 }
1016 }
1017 else // previous definition of just one composite
1018 {
1019
1020 // find the non comment portion of the body.
1021 TiXmlNode *elementChild = element->FirstChild();
1022 while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1023 {
1024 elementChild = elementChild->NextSibling();
1025 }
1026
1027 ASSERTL0(elementChild, "Unable to read DOMAIN body.");
1028 std::string elementStr = elementChild->ToText()->ValueStr();
1029
1030 elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
1031
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);
1035
1036 ASSERTL0(!indxStr.empty(),
1037 "Unable to read domain's composite index (index missing?).");
1038
1039 // Read the domain composites.
1040 // Parse the composites into a list.
1041 std::map<int, CompositeSharedPtr> fullDomain;
1042 m_meshGraph->GetCompositeList(indxStr, fullDomain);
1043 domain[0] = fullDomain;
1044
1045 ASSERTL0(
1046 !domain[0].empty(),
1047 (std::string("Unable to obtain domain's referenced composite: ") +
1048 indxStr)
1049 .c_str());
1050 }
1051}
1052
1054{
1055 auto &segGeoms = m_meshGraph->GetAllSegGeoms();
1056 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
1057 int spaceDimension = m_meshGraph->GetSpaceDimension();
1058
1059 CurveMap::iterator it;
1060
1061 /// Look for elements in ELEMENT block.
1062 TiXmlElement *field = m_xmlGeom->FirstChildElement("EDGE");
1063
1064 ASSERTL0(field, "Unable to find EDGE tag in file.");
1065
1066 /// All elements are of the form: "<E ID="#"> ... </E>", with
1067 /// ? being the element type.
1068 /// Read the ID field first.
1069 TiXmlElement *edgeElement = field->FirstChildElement("E");
1070
1071 /// Since all edge data is one big text block, we need to accumulate
1072 /// all TINYXML_TEXT data and then parse it. This approach effectively
1073 /// skips
1074 /// all comments or other node types since we only care about the
1075 /// edge list. We cannot handle missing edge numbers as we could
1076 /// with missing element numbers due to the text block format.
1077 std::string edgeStr;
1078 int indx;
1079
1080 while (edgeElement)
1081 {
1082 int err = edgeElement->QueryIntAttribute("ID", &indx);
1083 ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
1084
1085 TiXmlNode *child = edgeElement->FirstChild();
1086 edgeStr.clear();
1087 if (child->Type() == TiXmlNode::TINYXML_TEXT)
1088 {
1089 edgeStr += child->ToText()->ValueStr();
1090 }
1091
1092 /// Now parse out the edges, three fields at a time.
1093 int vertex1, vertex2;
1094 std::istringstream edgeDataStrm(edgeStr.c_str());
1095
1096 try
1097 {
1098 while (!edgeDataStrm.fail())
1099 {
1100 edgeDataStrm >> vertex1 >> vertex2;
1101
1102 // Must check after the read because we
1103 // may be at the end and not know it. If
1104 // we are at the end we will add a
1105 // duplicate of the last entry if we don't
1106 // check here.
1107 if (!edgeDataStrm.fail())
1108 {
1109 PointGeomSharedPtr vertices[2] = {
1110 m_meshGraph->GetVertex(vertex1),
1111 m_meshGraph->GetVertex(vertex2)};
1112 it = curvedEdges.find(indx);
1113
1114 if (it == curvedEdges.end())
1115 {
1116 segGeoms[indx] =
1118 indx, spaceDimension, vertices);
1119 }
1120 else
1121 {
1122 segGeoms[indx] =
1124 indx, spaceDimension, vertices, it->second);
1125 }
1126 }
1127 }
1128 }
1129 catch (...)
1130 {
1131 NEKERROR(
1133 (std::string("Unable to read edge data: ") + edgeStr).c_str());
1134 }
1135
1136 edgeElement = edgeElement->NextSiblingElement("E");
1137 }
1138}
1139
1141{
1142 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
1143 auto &triGeoms = m_meshGraph->GetAllTriGeoms();
1144 auto &quadGeoms = m_meshGraph->GetAllQuadGeoms();
1145
1146 /// Look for elements in FACE block.
1147 TiXmlElement *field = m_xmlGeom->FirstChildElement("FACE");
1148
1149 ASSERTL0(field, "Unable to find FACE tag in file.");
1150
1151 /// All faces are of the form: "<? ID="#"> ... </?>", with
1152 /// ? being an element type (either Q or T).
1153 /// They might be in compressed format and so then need upacking.
1154
1155 TiXmlElement *element = field->FirstChildElement();
1156 CurveMap::iterator it;
1157
1158 while (element)
1159 {
1160 std::string elementType(element->ValueStr());
1161
1162 ASSERTL0(elementType == "Q" || elementType == "T",
1163 (std::string("Unknown 3D face type: ") + elementType).c_str());
1164
1165 /// Read id attribute.
1166 int indx;
1167 int err = element->QueryIntAttribute("ID", &indx);
1168 ASSERTL0(err == TIXML_SUCCESS, "Unable to read face attribute ID.");
1169
1170 /// See if this face has curves.
1171 it = curvedFaces.find(indx);
1172
1173 /// Read text element description.
1174 TiXmlNode *elementChild = element->FirstChild();
1175 std::string elementStr;
1176 while (elementChild)
1177 {
1178 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1179 {
1180 elementStr += elementChild->ToText()->ValueStr();
1181 }
1182 elementChild = elementChild->NextSibling();
1183 }
1184
1185 ASSERTL0(!elementStr.empty(), "Unable to read face description body.");
1186
1187 /// Parse out the element components corresponding to type of
1188 /// element.
1189 if (elementType == "T")
1190 {
1191 // Read three edge numbers
1192 int edge1, edge2, edge3;
1193 std::istringstream elementDataStrm(elementStr.c_str());
1194
1195 try
1196 {
1197 elementDataStrm >> edge1;
1198 elementDataStrm >> edge2;
1199 elementDataStrm >> edge3;
1200
1201 ASSERTL0(
1202 !elementDataStrm.fail(),
1203 (std::string("Unable to read face data for TRIANGLE: ") +
1204 elementStr)
1205 .c_str());
1206
1207 /// Create a TriGeom to hold the new definition.
1209 m_meshGraph->GetSegGeom(edge1),
1210 m_meshGraph->GetSegGeom(edge2),
1211 m_meshGraph->GetSegGeom(edge3)};
1212
1213 if (it == curvedFaces.end())
1214 {
1215 triGeoms[indx] =
1217 }
1218 else
1219 {
1221 indx, edges, it->second);
1222 }
1223 triGeoms[indx]->SetGlobalID(indx);
1224 }
1225 catch (...)
1226 {
1227 NEKERROR(
1229 (std::string("Unable to read face data for TRIANGLE: ") +
1230 elementStr)
1231 .c_str());
1232 }
1233 }
1234 else if (elementType == "Q")
1235 {
1236 // Read four edge numbers
1237 int edge1, edge2, edge3, edge4;
1238 std::istringstream elementDataStrm(elementStr.c_str());
1239
1240 try
1241 {
1242 elementDataStrm >> edge1;
1243 elementDataStrm >> edge2;
1244 elementDataStrm >> edge3;
1245 elementDataStrm >> edge4;
1246
1247 ASSERTL0(!elementDataStrm.fail(),
1248 (std::string("Unable to read face data for QUAD: ") +
1249 elementStr)
1250 .c_str());
1251
1252 /// Create a QuadGeom to hold the new definition.
1254 m_meshGraph->GetSegGeom(edge1),
1255 m_meshGraph->GetSegGeom(edge2),
1256 m_meshGraph->GetSegGeom(edge3),
1257 m_meshGraph->GetSegGeom(edge4)};
1258
1259 QuadGeomSharedPtr quadgeom;
1260
1261 if (it == curvedFaces.end())
1262 {
1263 quadGeoms[indx] =
1265 }
1266 else
1267 {
1268 quadGeoms[indx] =
1270 it->second);
1271 }
1272 quadGeoms[indx]->SetGlobalID(indx);
1273 }
1274 catch (...)
1275 {
1277 (std::string("Unable to read face data for QUAD: ") +
1278 elementStr)
1279 .c_str());
1280 }
1281 }
1282 // Keep looking
1283 element = element->NextSiblingElement();
1284 }
1285}
1286
1288{
1289 int meshDimension = m_meshGraph->GetMeshDimension();
1290
1291 switch (meshDimension)
1292 {
1293 case 1:
1295 break;
1296 case 2:
1298 break;
1299 case 3:
1301 break;
1302 }
1303}
1304
1306{
1307 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
1308 auto &segGeoms = m_meshGraph->GetAllSegGeoms();
1309 int spaceDimension = m_meshGraph->GetSpaceDimension();
1310
1311 TiXmlElement *field = nullptr;
1312
1313 /// Look for elements in ELEMENT block.
1314 field = m_xmlGeom->FirstChildElement("ELEMENT");
1315
1316 ASSERTL0(field, "Unable to find ELEMENT tag in file.");
1317
1318 /// All elements are of the form: "<S ID = n> ... </S>", with
1319 /// ? being the element type.
1320
1321 TiXmlElement *segment = field->FirstChildElement("S");
1322 CurveMap::iterator it;
1323
1324 while (segment)
1325 {
1326 int indx;
1327 int err = segment->QueryIntAttribute("ID", &indx);
1328 ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
1329
1330 TiXmlNode *elementChild = segment->FirstChild();
1331 while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1332 {
1333 elementChild = elementChild->NextSibling();
1334 }
1335
1336 ASSERTL0(elementChild, "Unable to read element description body.");
1337 std::string elementStr = elementChild->ToText()->ValueStr();
1338
1339 /// Parse out the element components corresponding to type of
1340 /// element.
1341 /// Read two vertex numbers
1342 int vertex1, vertex2;
1343 std::istringstream elementDataStrm(elementStr.c_str());
1344
1345 try
1346 {
1347 elementDataStrm >> vertex1;
1348 elementDataStrm >> vertex2;
1349
1350 ASSERTL0(!elementDataStrm.fail(),
1351 (std::string("Unable to read element data for SEGMENT: ") +
1352 elementStr)
1353 .c_str());
1354
1355 PointGeomSharedPtr vertices[2] = {m_meshGraph->GetVertex(vertex1),
1356 m_meshGraph->GetVertex(vertex2)};
1357 it = curvedEdges.find(indx);
1358
1359 if (it == curvedEdges.end())
1360 {
1362 indx, spaceDimension, vertices);
1363 }
1364 else
1365 {
1367 indx, spaceDimension, vertices, it->second);
1368 }
1369 segGeoms[indx]->SetGlobalID(indx);
1370 }
1371 catch (...)
1372 {
1374 (std::string("Unable to read element data for segment: ") +
1375 elementStr)
1376 .c_str());
1377 }
1378 /// Keep looking for additional segments
1379 segment = segment->NextSiblingElement("S");
1380 }
1381}
1382
1384{
1385 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
1386 auto &triGeoms = m_meshGraph->GetAllTriGeoms();
1387 auto &quadGeoms = m_meshGraph->GetAllQuadGeoms();
1388
1389 /// Look for elements in ELEMENT block.
1390 TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
1391
1392 ASSERTL0(field, "Unable to find ELEMENT tag in file.");
1393
1394 // Set up curve map for curved elements on an embedded manifold.
1395 CurveMap::iterator it;
1396
1397 /// All elements are of the form: "<? ID="#"> ... </?>", with
1398 /// ? being the element type.
1399
1400 TiXmlElement *element = field->FirstChildElement();
1401
1402 while (element)
1403 {
1404 std::string elementType(element->ValueStr());
1405
1406 ASSERTL0(
1407 elementType == "Q" || elementType == "T",
1408 (std::string("Unknown 2D element type: ") + elementType).c_str());
1409
1410 /// Read id attribute.
1411 int indx;
1412 int err = element->QueryIntAttribute("ID", &indx);
1413 ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
1414
1415 it = curvedFaces.find(indx);
1416
1417 /// Read text element description.
1418 TiXmlNode *elementChild = element->FirstChild();
1419 std::string elementStr;
1420 while (elementChild)
1421 {
1422 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1423 {
1424 elementStr += elementChild->ToText()->ValueStr();
1425 }
1426 elementChild = elementChild->NextSibling();
1427 }
1428
1429 ASSERTL0(!elementStr.empty(),
1430 "Unable to read element description body.");
1431
1432 /// Parse out the element components corresponding to type of
1433 /// element.
1434 if (elementType == "T")
1435 {
1436 // Read three edge numbers
1437 int edge1, edge2, edge3;
1438 std::istringstream elementDataStrm(elementStr.c_str());
1439
1440 try
1441 {
1442 elementDataStrm >> edge1;
1443 elementDataStrm >> edge2;
1444 elementDataStrm >> edge3;
1445
1446 ASSERTL0(
1447 !elementDataStrm.fail(),
1448 (std::string("Unable to read element data for TRIANGLE: ") +
1449 elementStr)
1450 .c_str());
1451
1452 /// Create a TriGeom to hold the new definition.
1454 m_meshGraph->GetSegGeom(edge1),
1455 m_meshGraph->GetSegGeom(edge2),
1456 m_meshGraph->GetSegGeom(edge3)};
1457
1458 if (it == curvedFaces.end())
1459 {
1460 triGeoms[indx] =
1462 }
1463 else
1464 {
1466 indx, edges, it->second);
1467 }
1468 triGeoms[indx]->SetGlobalID(indx);
1469 }
1470 catch (...)
1471 {
1472 NEKERROR(
1474 (std::string("Unable to read element data for TRIANGLE: ") +
1475 elementStr)
1476 .c_str());
1477 }
1478 }
1479 else if (elementType == "Q")
1480 {
1481 // Read four edge numbers
1482 int edge1, edge2, edge3, edge4;
1483 std::istringstream elementDataStrm(elementStr.c_str());
1484
1485 try
1486 {
1487 elementDataStrm >> edge1;
1488 elementDataStrm >> edge2;
1489 elementDataStrm >> edge3;
1490 elementDataStrm >> edge4;
1491
1492 ASSERTL0(
1493 !elementDataStrm.fail(),
1494 (std::string("Unable to read element data for QUAD: ") +
1495 elementStr)
1496 .c_str());
1497
1498 /// Create a QuadGeom to hold the new definition.
1500 m_meshGraph->GetSegGeom(edge1),
1501 m_meshGraph->GetSegGeom(edge2),
1502 m_meshGraph->GetSegGeom(edge3),
1503 m_meshGraph->GetSegGeom(edge4)};
1504
1505 QuadGeomSharedPtr quadgeom;
1506 if (it == curvedFaces.end())
1507 {
1508 quadGeoms[indx] =
1510 }
1511 else
1512 {
1513 quadGeoms[indx] =
1515 it->second);
1516 }
1517 quadGeoms[indx]->SetGlobalID(indx);
1518 }
1519 catch (...)
1520 {
1521 NEKERROR(
1523 (std::string("Unable to read element data for QUAD: ") +
1524 elementStr)
1525 .c_str());
1526 }
1527 }
1528 /// Keep looking
1529 element = element->NextSiblingElement();
1530 }
1531}
1532
1534{
1535 auto &tetGeoms = m_meshGraph->GetAllTetGeoms();
1536 auto &pyrGeoms = m_meshGraph->GetAllPyrGeoms();
1537 auto &prismGeoms = m_meshGraph->GetAllPrismGeoms();
1538 auto &hexGeoms = m_meshGraph->GetAllHexGeoms();
1539
1540 /// Look for elements in ELEMENT block.
1541 TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
1542
1543 ASSERTL0(field, "Unable to find ELEMENT tag in file.");
1544
1545 /// All elements are of the form: "<? ID="#"> ... </?>", with
1546 /// ? being the element type.
1547
1548 TiXmlElement *element = field->FirstChildElement();
1549
1550 while (element)
1551 {
1552 std::string elementType(element->ValueStr());
1553
1554 // A - tet, P - pyramid, R - prism, H - hex
1555 ASSERTL0(
1556 elementType == "A" || elementType == "P" || elementType == "R" ||
1557 elementType == "H",
1558 (std::string("Unknown 3D element type: ") + elementType).c_str());
1559
1560 /// Read id attribute.
1561 int indx;
1562 int err = element->QueryIntAttribute("ID", &indx);
1563 ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
1564
1565 /// Read text element description.
1566 TiXmlNode *elementChild = element->FirstChild();
1567 std::string elementStr;
1568 while (elementChild)
1569 {
1570 if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1571 {
1572 elementStr += elementChild->ToText()->ValueStr();
1573 }
1574 elementChild = elementChild->NextSibling();
1575 }
1576
1577 ASSERTL0(!elementStr.empty(),
1578 "Unable to read element description body.");
1579
1580 std::istringstream elementDataStrm(elementStr.c_str());
1581
1582 /// Parse out the element components corresponding to type of
1583 /// element.
1584
1585 // Tetrahedral
1586 if (elementType == "A")
1587 {
1588 try
1589 {
1590 /// Create arrays for the tri and quad faces.
1591 const int kNfaces = TetGeom::kNfaces;
1592 const int kNtfaces = TetGeom::kNtfaces;
1593 const int kNqfaces = TetGeom::kNqfaces;
1594 TriGeomSharedPtr tfaces[kNtfaces];
1595 int Ntfaces = 0;
1596 int Nqfaces = 0;
1597
1598 /// Fill the arrays and make sure there aren't too many
1599 /// faces.
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++)
1605 {
1606 int faceID;
1607 elementDataStrm >> faceID;
1608 Geometry2DSharedPtr face =
1609 m_meshGraph->GetGeometry2D(faceID);
1610 if (face == Geometry2DSharedPtr() ||
1611 (face->GetShapeType() != LibUtilities::eTriangle &&
1612 face->GetShapeType() != LibUtilities::eQuadrilateral))
1613 {
1614 std::stringstream errorstring;
1615 errorstring << "Element " << indx
1616 << " has invalid face: " << faceID;
1617 ASSERTL0(false, errorstring.str().c_str());
1618 }
1619 else if (face->GetShapeType() == LibUtilities::eTriangle)
1620 {
1621 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1622 tfaces[Ntfaces++] =
1623 std::static_pointer_cast<TriGeom>(face);
1624 }
1625 else if (face->GetShapeType() ==
1627 {
1628 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1629 }
1630 }
1631
1632 /// Make sure all of the face indicies could be read, and
1633 /// that there weren't too few.
1634 ASSERTL0(!elementDataStrm.fail(),
1635 (std::string(
1636 "Unable to read element data for TETRAHEDRON: ") +
1637 elementStr)
1638 .c_str());
1639 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1640 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1641
1642 TetGeomSharedPtr tetgeom(
1644
1645 tetGeoms[indx] = tetgeom;
1646 m_meshGraph->PopulateFaceToElMap(tetgeom, kNfaces);
1647 }
1648 catch (...)
1649 {
1651 (std::string(
1652 "Unable to read element data for TETRAHEDRON: ") +
1653 elementStr)
1654 .c_str());
1655 }
1656 }
1657 // Pyramid
1658 else if (elementType == "P")
1659 {
1660 try
1661 {
1662 /// Create arrays for the tri and quad faces.
1663 const int kNfaces = PyrGeom::kNfaces;
1664 const int kNtfaces = PyrGeom::kNtfaces;
1665 const int kNqfaces = PyrGeom::kNqfaces;
1666 Geometry2DSharedPtr faces[kNfaces];
1667 int Nfaces = 0;
1668 int Ntfaces = 0;
1669 int Nqfaces = 0;
1670
1671 /// Fill the arrays and make sure there aren't too many
1672 /// faces.
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++)
1678 {
1679 int faceID;
1680 elementDataStrm >> faceID;
1681 Geometry2DSharedPtr face =
1682 m_meshGraph->GetGeometry2D(faceID);
1683 if (face == Geometry2DSharedPtr() ||
1684 (face->GetShapeType() != LibUtilities::eTriangle &&
1685 face->GetShapeType() != LibUtilities::eQuadrilateral))
1686 {
1687 std::stringstream errorstring;
1688 errorstring << "Element " << indx
1689 << " has invalid face: " << faceID;
1690 ASSERTL0(false, errorstring.str().c_str());
1691 }
1692 else if (face->GetShapeType() == LibUtilities::eTriangle)
1693 {
1694 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1695 faces[Nfaces++] =
1696 std::static_pointer_cast<TriGeom>(face);
1697 Ntfaces++;
1698 }
1699 else if (face->GetShapeType() ==
1701 {
1702 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1703 faces[Nfaces++] =
1704 std::static_pointer_cast<QuadGeom>(face);
1705 Nqfaces++;
1706 }
1707 }
1708
1709 /// Make sure all of the face indicies could be read, and
1710 /// that there weren't too few.
1711 ASSERTL0(
1712 !elementDataStrm.fail(),
1713 (std::string("Unable to read element data for PYRAMID: ") +
1714 elementStr)
1715 .c_str());
1716 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1717 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1718
1719 PyrGeomSharedPtr pyrgeom(
1721
1722 pyrGeoms[indx] = pyrgeom;
1723 m_meshGraph->PopulateFaceToElMap(pyrgeom, kNfaces);
1724 }
1725 catch (...)
1726 {
1727 NEKERROR(
1729 (std::string("Unable to read element data for PYRAMID: ") +
1730 elementStr)
1731 .c_str());
1732 }
1733 }
1734 // Prism
1735 else if (elementType == "R")
1736 {
1737 try
1738 {
1739 /// Create arrays for the tri and quad faces.
1740 const int kNfaces = PrismGeom::kNfaces;
1741 const int kNtfaces = PrismGeom::kNtfaces;
1742 const int kNqfaces = PrismGeom::kNqfaces;
1743 Geometry2DSharedPtr faces[kNfaces];
1744 int Ntfaces = 0;
1745 int Nqfaces = 0;
1746 int Nfaces = 0;
1747
1748 /// Fill the arrays and make sure there aren't too many
1749 /// faces.
1750 std::stringstream errorstring;
1751 errorstring << "Element " << indx << " must have " << kNtfaces
1752 << " triangle face(s), and " << kNqfaces
1753 << " quadrilateral face(s).";
1754
1755 for (int i = 0; i < kNfaces; i++)
1756 {
1757 int faceID;
1758 elementDataStrm >> faceID;
1759 Geometry2DSharedPtr face =
1760 m_meshGraph->GetGeometry2D(faceID);
1761 if (face == Geometry2DSharedPtr() ||
1762 (face->GetShapeType() != LibUtilities::eTriangle &&
1763 face->GetShapeType() != LibUtilities::eQuadrilateral))
1764 {
1765 std::stringstream errorstring;
1766 errorstring << "Element " << indx
1767 << " has invalid face: " << faceID;
1768 ASSERTL0(false, errorstring.str().c_str());
1769 }
1770 else if (face->GetShapeType() == LibUtilities::eTriangle)
1771 {
1772 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1773 faces[Nfaces++] =
1774 std::static_pointer_cast<TriGeom>(face);
1775 Ntfaces++;
1776 }
1777 else if (face->GetShapeType() ==
1779 {
1780 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1781 faces[Nfaces++] =
1782 std::static_pointer_cast<QuadGeom>(face);
1783 Nqfaces++;
1784 }
1785 }
1786
1787 /// Make sure all of the face indices could be read, and
1788 /// that there weren't too few.
1789 ASSERTL0(
1790 !elementDataStrm.fail(),
1791 (std::string("Unable to read element data for PRISM: ") +
1792 elementStr)
1793 .c_str());
1794 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1795 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1796
1797 PrismGeomSharedPtr prismgeom(
1799
1800 prismGeoms[indx] = prismgeom;
1801 m_meshGraph->PopulateFaceToElMap(prismgeom, kNfaces);
1802 }
1803 catch (...)
1804 {
1805 NEKERROR(
1807 (std::string("Unable to read element data for PRISM: ") +
1808 elementStr)
1809 .c_str());
1810 }
1811 }
1812 // Hexahedral
1813 else if (elementType == "H")
1814 {
1815 try
1816 {
1817 /// Create arrays for the tri and quad faces.
1818 const int kNfaces = HexGeom::kNfaces;
1819 const int kNtfaces = HexGeom::kNtfaces;
1820 const int kNqfaces = HexGeom::kNqfaces;
1821 // TriGeomSharedPtr tfaces[kNtfaces];
1822 QuadGeomSharedPtr qfaces[kNqfaces];
1823 int Ntfaces = 0;
1824 int Nqfaces = 0;
1825
1826 /// Fill the arrays and make sure there aren't too many
1827 /// faces.
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++)
1833 {
1834 int faceID;
1835 elementDataStrm >> faceID;
1836 Geometry2DSharedPtr face =
1837 m_meshGraph->GetGeometry2D(faceID);
1838 if (face == Geometry2DSharedPtr() ||
1839 (face->GetShapeType() != LibUtilities::eTriangle &&
1840 face->GetShapeType() != LibUtilities::eQuadrilateral))
1841 {
1842 std::stringstream errorstring;
1843 errorstring << "Element " << indx
1844 << " has invalid face: " << faceID;
1845 ASSERTL0(false, errorstring.str().c_str());
1846 }
1847 else if (face->GetShapeType() == LibUtilities::eTriangle)
1848 {
1849 ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
1850 }
1851 else if (face->GetShapeType() ==
1853 {
1854 ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
1855 qfaces[Nqfaces++] =
1856 std::static_pointer_cast<QuadGeom>(face);
1857 }
1858 }
1859
1860 /// Make sure all of the face indicies could be read, and
1861 /// that there weren't too few.
1862 ASSERTL0(!elementDataStrm.fail(),
1863 (std::string(
1864 "Unable to read element data for HEXAHEDRAL: ") +
1865 elementStr)
1866 .c_str());
1867 ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
1868 ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
1869
1870 HexGeomSharedPtr hexgeom(
1872
1873 hexGeoms[indx] = hexgeom;
1874 m_meshGraph->PopulateFaceToElMap(hexgeom, kNfaces);
1875 }
1876 catch (...)
1877 {
1879 (std::string(
1880 "Unable to read element data for HEXAHEDRAL: ") +
1881 elementStr)
1882 .c_str());
1883 }
1884 }
1885 /// Keep looking
1886 element = element->NextSiblingElement();
1887 }
1888}
1889
1891{
1892 auto &meshComposites = m_meshGraph->GetComposites();
1893 auto &compositesLabels = m_meshGraph->GetCompositesLabels();
1894
1895 TiXmlElement *field = nullptr;
1896
1897 /// Look for elements in ELEMENT block.
1898 field = m_xmlGeom->FirstChildElement("COMPOSITE");
1899
1900 ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
1901
1902 TiXmlElement *node = field->FirstChildElement("C");
1903
1904 // Sequential counter for the composite numbers.
1905 int nextCompositeNumber = -1;
1906
1907 while (node)
1908 {
1909 /// All elements are of the form: "<? ID="#"> ... </?>", with
1910 /// ? being the element type.
1911
1912 nextCompositeNumber++;
1913
1914 int indx;
1915 int err = node->QueryIntAttribute("ID", &indx);
1916 ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
1917 // ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with
1918 // zero and be sequential.");
1919
1920 TiXmlNode *compositeChild = node->FirstChild();
1921 // This is primarily to skip comments that may be present.
1922 // Comments appear as nodes just like elements.
1923 // We are specifically looking for text in the body
1924 // of the definition.
1925 while (compositeChild &&
1926 compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
1927 {
1928 compositeChild = compositeChild->NextSibling();
1929 }
1930
1931 ASSERTL0(compositeChild, "Unable to read composite definition body.");
1932 std::string compositeStr = compositeChild->ToText()->ValueStr();
1933
1934 /// Parse out the element components corresponding to type of element.
1935
1936 std::istringstream compositeDataStrm(compositeStr.c_str());
1937
1938 try
1939 {
1940 bool first = true;
1941 std::string prevCompositeElementStr;
1942
1943 while (!compositeDataStrm.fail())
1944 {
1945 std::string compositeElementStr;
1946 compositeDataStrm >> compositeElementStr;
1947
1948 if (!compositeDataStrm.fail())
1949 {
1950 if (first)
1951 {
1952 first = false;
1953 meshComposites[indx] =
1955 }
1956
1957 if (compositeElementStr.length() > 0)
1958 {
1959 ResolveGeomRef(prevCompositeElementStr,
1960 compositeElementStr,
1961 meshComposites[indx]);
1962 }
1963 prevCompositeElementStr = compositeElementStr;
1964 }
1965 }
1966 }
1967 catch (...)
1968 {
1969 NEKERROR(
1971 (std::string("Unable to read COMPOSITE data for composite: ") +
1972 compositeStr)
1973 .c_str());
1974 }
1975
1976 // Read optional name as string and save to compositeLabels if exists
1977 std::string name;
1978 err = node->QueryStringAttribute("NAME", &name);
1979 if (err == TIXML_SUCCESS)
1980 {
1981 compositesLabels[indx] = name;
1982 }
1983
1984 /// Keep looking for additional composite definitions.
1985 node = node->NextSiblingElement("C");
1986 }
1987
1988 ASSERTL0(nextCompositeNumber >= 0,
1989 "At least one composite must be specified.");
1990}
1991
1992void MeshGraphIOXml::ResolveGeomRef(const std::string &prevToken,
1993 const std::string &token,
1994 CompositeSharedPtr &composite)
1995{
1996 int meshDimension = m_meshGraph->GetMeshDimension();
1997
1998 switch (meshDimension)
1999 {
2000 case 1:
2001 ResolveGeomRef1D(prevToken, token, composite);
2002 break;
2003 case 2:
2004 ResolveGeomRef2D(prevToken, token, composite);
2005 break;
2006 case 3:
2007 ResolveGeomRef3D(prevToken, token, composite);
2008 break;
2009 }
2010}
2011
2012void MeshGraphIOXml::ResolveGeomRef1D(const std::string &prevToken,
2013 const std::string &token,
2014 CompositeSharedPtr &composite)
2015{
2016 auto &vertSet = m_meshGraph->GetAllPointGeoms();
2017 auto &segGeoms = m_meshGraph->GetAllSegGeoms();
2018
2019 try
2020 {
2021 std::istringstream tokenStream(token);
2022 std::istringstream prevTokenStream(prevToken);
2023
2024 char type;
2025 char prevType;
2026
2027 tokenStream >> type;
2028
2029 std::string::size_type indxBeg = token.find_first_of('[') + 1;
2030 std::string::size_type indxEnd = token.find_last_of(']') - 1;
2031
2032 ASSERTL0(
2033 indxBeg <= indxEnd,
2034 (std::string("Error reading index definition:") + token).c_str());
2035
2036 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2037
2038 typedef std::vector<unsigned int> SeqVectorType;
2039 SeqVectorType seqVector;
2040
2041 if (!ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector))
2042 {
2044 (std::string("Ill-formed sequence definition: ") + indxStr)
2045 .c_str());
2046 }
2047
2048 prevTokenStream >> prevType;
2049
2050 // All composites must be of the same dimension.
2051 bool validSequence =
2052 (prevToken.empty() || // No previous, then current is just fine.
2053 (type == 'V' && prevType == 'V') ||
2054 (type == 'S' && prevType == 'S'));
2055
2056 ASSERTL0(validSequence,
2057 std::string("Invalid combination of composite items: ") +
2058 type + " and " + prevType + ".");
2059
2060 switch (type)
2061 {
2062 case 'V': // Vertex
2063 for (SeqVectorType::iterator iter = seqVector.begin();
2064 iter != seqVector.end(); ++iter)
2065 {
2066 if (vertSet.find(*iter) == vertSet.end())
2067 {
2069 "Unknown vertex index: " +
2070 std::to_string(*iter));
2071 }
2072 else
2073 {
2074 composite->m_geomVec.push_back(vertSet[*iter]);
2075 }
2076 }
2077 break;
2078
2079 case 'S': // Segment
2080 for (SeqVectorType::iterator iter = seqVector.begin();
2081 iter != seqVector.end(); ++iter)
2082 {
2083 if (segGeoms.find(*iter) == segGeoms.end())
2084 {
2086 "Unknown segment index: " +
2087 std::to_string(*iter));
2088 }
2089 else
2090 {
2091 composite->m_geomVec.push_back(segGeoms[*iter]);
2092 }
2093 }
2094 break;
2095
2096 default:
2098 "Unrecognized composite token: " + token);
2099 }
2100 }
2101 catch (...)
2102 {
2104 "Problem processing composite token: " + token);
2105 }
2106
2107 return;
2108}
2109
2110void MeshGraphIOXml::ResolveGeomRef2D(const std::string &prevToken,
2111 const std::string &token,
2112 CompositeSharedPtr &composite)
2113{
2114 PointGeomMap &vertSet = m_meshGraph->GetAllPointGeoms();
2115 SegGeomMap &segGeoms = m_meshGraph->GetAllSegGeoms();
2116 TriGeomMap &triGeoms = m_meshGraph->GetAllTriGeoms();
2117 QuadGeomMap &quadGeoms = m_meshGraph->GetAllQuadGeoms();
2118
2119 try
2120 {
2121 std::istringstream tokenStream(token);
2122 std::istringstream prevTokenStream(prevToken);
2123
2124 char type;
2125 char prevType;
2126
2127 tokenStream >> type;
2128
2129 std::string::size_type indxBeg = token.find_first_of('[') + 1;
2130 std::string::size_type indxEnd = token.find_last_of(']') - 1;
2131
2132 ASSERTL0(
2133 indxBeg <= indxEnd,
2134 (std::string("Error reading index definition:") + token).c_str());
2135
2136 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2137 std::vector<unsigned int> seqVector;
2138 std::vector<unsigned int>::iterator seqIter;
2139
2140 bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
2141
2142 ASSERTL0(err,
2143 (std::string("Error reading composite elements: ") + indxStr)
2144 .c_str());
2145
2146 prevTokenStream >> prevType;
2147
2148 // All composites must be of the same dimension.
2149 bool validSequence =
2150 (prevToken.empty() || // No previous, then current is just fine.
2151 (type == 'V' && prevType == 'V') ||
2152 (type == 'E' && prevType == 'E') ||
2153 ((type == 'T' || type == 'Q') &&
2154 (prevType == 'T' || prevType == 'Q')));
2155
2156 ASSERTL0(validSequence,
2157 std::string("Invalid combination of composite items: ") +
2158 type + " and " + prevType + ".");
2159
2160 switch (type)
2161 {
2162 case 'E': // Edge
2163 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2164 ++seqIter)
2165 {
2166 if (segGeoms.find(*seqIter) == segGeoms.end())
2167 {
2169 "Unknown edge index: " +
2170 std::to_string(*seqIter));
2171 }
2172 else
2173 {
2174 composite->m_geomVec.push_back(segGeoms[*seqIter]);
2175 }
2176 }
2177 break;
2178
2179 case 'T': // Triangle
2180 {
2181 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2182 ++seqIter)
2183 {
2184 if (triGeoms.count(*seqIter) == 0)
2185 {
2187 "Unknown triangle index: " +
2188 std::to_string(*seqIter));
2189 }
2190 else
2191 {
2192 if (m_meshGraph->CheckRange(*triGeoms[*seqIter]))
2193 {
2194 composite->m_geomVec.push_back(triGeoms[*seqIter]);
2195 }
2196 }
2197 }
2198 }
2199 break;
2200
2201 case 'Q': // Quad
2202 {
2203 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2204 ++seqIter)
2205 {
2206 if (quadGeoms.count(*seqIter) == 0)
2207 {
2208 NEKERROR(
2210 "Unknown quad index: " + std::to_string(*seqIter) +
2211 " in Composite section");
2212 }
2213 else
2214 {
2215 if (m_meshGraph->CheckRange(*quadGeoms[*seqIter]))
2216 {
2217 composite->m_geomVec.push_back(quadGeoms[*seqIter]);
2218 }
2219 }
2220 }
2221 }
2222 break;
2223
2224 case 'V': // Vertex
2225 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2226 ++seqIter)
2227 {
2228 if (*seqIter >= vertSet.size())
2229 {
2231 "Unknown vertex index: " +
2232 std::to_string(*seqIter));
2233 }
2234 else
2235 {
2236 composite->m_geomVec.push_back(vertSet[*seqIter]);
2237 }
2238 }
2239 break;
2240
2241 default:
2243 "Unrecognized composite token: " + token);
2244 }
2245 }
2246 catch (...)
2247 {
2249 "Problem processing composite token: " + token);
2250 }
2251
2252 return;
2253}
2254
2255void MeshGraphIOXml::ResolveGeomRef3D(const std::string &prevToken,
2256 const std::string &token,
2257 CompositeSharedPtr &composite)
2258{
2259 PointGeomMap &vertSet = m_meshGraph->GetAllPointGeoms();
2260 SegGeomMap &segGeoms = m_meshGraph->GetAllSegGeoms();
2261 TriGeomMap &triGeoms = m_meshGraph->GetAllTriGeoms();
2262 QuadGeomMap &quadGeoms = m_meshGraph->GetAllQuadGeoms();
2263 TetGeomMap &tetGeoms = m_meshGraph->GetAllTetGeoms();
2264 PyrGeomMap &pyrGeoms = m_meshGraph->GetAllPyrGeoms();
2265 PrismGeomMap &prismGeoms = m_meshGraph->GetAllPrismGeoms();
2266 HexGeomMap &hexGeoms = m_meshGraph->GetAllHexGeoms();
2267
2268 try
2269 {
2270 std::istringstream tokenStream(token);
2271 std::istringstream prevTokenStream(prevToken);
2272
2273 char type;
2274 char prevType;
2275
2276 tokenStream >> type;
2277
2278 std::string::size_type indxBeg = token.find_first_of('[') + 1;
2279 std::string::size_type indxEnd = token.find_last_of(']') - 1;
2280
2281 ASSERTL0(indxBeg <= indxEnd,
2282 "Error reading index definition: " + token);
2283
2284 std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
2285
2286 std::vector<unsigned int> seqVector;
2287 std::vector<unsigned int>::iterator seqIter;
2288
2289 bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
2290
2291 ASSERTL0(err, "Error reading composite elements: " + indxStr);
2292
2293 prevTokenStream >> prevType;
2294
2295 // All composites must be of the same dimension. This map makes things
2296 // clean to compare.
2297 std::map<char, int> typeMap;
2298 typeMap['V'] = 1; // Vertex
2299 typeMap['E'] = 1; // Edge
2300 typeMap['T'] = 2; // Triangle
2301 typeMap['Q'] = 2; // Quad
2302 typeMap['A'] = 3; // Tet
2303 typeMap['P'] = 3; // Pyramid
2304 typeMap['R'] = 3; // Prism
2305 typeMap['H'] = 3; // Hex
2306
2307 // Make sure only geoms of the same dimension are combined.
2308 bool validSequence =
2309 (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
2310
2311 ASSERTL0(validSequence,
2312 std::string("Invalid combination of composite items: ") +
2313 type + " and " + prevType + ".");
2314
2315 switch (type)
2316 {
2317 case 'V': // Vertex
2318 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2319 ++seqIter)
2320 {
2321 if (vertSet.find(*seqIter) == vertSet.end())
2322 {
2324 "Unknown vertex index: " +
2325 std::to_string(*seqIter));
2326 }
2327 else
2328 {
2329 composite->m_geomVec.push_back(vertSet[*seqIter]);
2330 }
2331 }
2332 break;
2333
2334 case 'E': // Edge
2335 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2336 ++seqIter)
2337 {
2338 if (segGeoms.find(*seqIter) == segGeoms.end())
2339 {
2341 "Unknown edge index: " +
2342 std::to_string(*seqIter));
2343 }
2344 else
2345 {
2346 composite->m_geomVec.push_back(segGeoms[*seqIter]);
2347 }
2348 }
2349 break;
2350
2351 case 'F': // Face
2352 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2353 ++seqIter)
2354 {
2355 Geometry2DSharedPtr face =
2356 m_meshGraph->GetGeometry2D(*seqIter);
2357 if (face == Geometry2DSharedPtr())
2358 {
2360 "Unknown face index: " +
2361 std::to_string(*seqIter));
2362 }
2363 else
2364 {
2365 if (m_meshGraph->CheckRange(*face))
2366 {
2367 composite->m_geomVec.push_back(face);
2368 }
2369 }
2370 }
2371 break;
2372
2373 case 'T': // Triangle
2374 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2375 ++seqIter)
2376 {
2377 if (triGeoms.find(*seqIter) == triGeoms.end())
2378 {
2380 "Unknown triangle index: " +
2381 std::to_string(*seqIter));
2382 }
2383 else
2384 {
2385 if (m_meshGraph->CheckRange(*triGeoms[*seqIter]))
2386 {
2387 composite->m_geomVec.push_back(triGeoms[*seqIter]);
2388 }
2389 }
2390 }
2391 break;
2392
2393 case 'Q': // Quad
2394 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2395 ++seqIter)
2396 {
2397 if (quadGeoms.find(*seqIter) == quadGeoms.end())
2398 {
2400 "Unknown quad index: " +
2401 std::to_string(*seqIter));
2402 }
2403 else
2404 {
2405 if (m_meshGraph->CheckRange(*quadGeoms[*seqIter]))
2406 {
2407 composite->m_geomVec.push_back(quadGeoms[*seqIter]);
2408 }
2409 }
2410 }
2411 break;
2412
2413 // Tetrahedron
2414 case 'A':
2415 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2416 ++seqIter)
2417 {
2418 if (tetGeoms.find(*seqIter) == tetGeoms.end())
2419 {
2421 "Unknown tet index: " +
2422 std::to_string(*seqIter));
2423 }
2424 else
2425 {
2426 if (m_meshGraph->CheckRange(*tetGeoms[*seqIter]))
2427 {
2428 composite->m_geomVec.push_back(tetGeoms[*seqIter]);
2429 }
2430 }
2431 }
2432 break;
2433
2434 // Pyramid
2435 case 'P':
2436 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2437 ++seqIter)
2438 {
2439 if (pyrGeoms.find(*seqIter) == pyrGeoms.end())
2440 {
2442 "Unknown pyramid index: " +
2443 std::to_string(*seqIter));
2444 }
2445 else
2446 {
2447 if (m_meshGraph->CheckRange(*pyrGeoms[*seqIter]))
2448 {
2449 composite->m_geomVec.push_back(pyrGeoms[*seqIter]);
2450 }
2451 }
2452 }
2453 break;
2454
2455 // Prism
2456 case 'R':
2457 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2458 ++seqIter)
2459 {
2460 if (prismGeoms.find(*seqIter) == prismGeoms.end())
2461 {
2463 "Unknown prism index: " +
2464 std::to_string(*seqIter));
2465 }
2466 else
2467 {
2468 if (m_meshGraph->CheckRange(*prismGeoms[*seqIter]))
2469 {
2470 composite->m_geomVec.push_back(
2471 prismGeoms[*seqIter]);
2472 }
2473 }
2474 }
2475 break;
2476
2477 // Hex
2478 case 'H':
2479 for (seqIter = seqVector.begin(); seqIter != seqVector.end();
2480 ++seqIter)
2481 {
2482 if (hexGeoms.find(*seqIter) == hexGeoms.end())
2483 {
2485 "Unknown hex index: " +
2486 std::to_string(*seqIter));
2487 }
2488 else
2489 {
2490 if (m_meshGraph->CheckRange(*hexGeoms[*seqIter]))
2491 {
2492 composite->m_geomVec.push_back(hexGeoms[*seqIter]);
2493 }
2494 }
2495 }
2496 break;
2497
2498 default:
2500 "Unrecognized composite token: " + token);
2501 }
2502 }
2503 catch (...)
2504 {
2506 "Problem processing composite token: " + token);
2507 }
2508
2509 return;
2510}
2511
2512void MeshGraphIOXml::v_WriteVertices(TiXmlElement *geomTag, PointGeomMap &verts)
2513{
2514 TiXmlElement *vertTag = new TiXmlElement("VERTEX");
2515
2516 for (auto &i : verts)
2517 {
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);
2525 }
2526
2527 geomTag->LinkEndChild(vertTag);
2528}
2529
2530void MeshGraphIOXml::v_WriteEdges(TiXmlElement *geomTag, SegGeomMap &edges)
2531{
2532 int meshDimension = m_meshGraph->GetMeshDimension();
2533
2534 TiXmlElement *edgeTag =
2535 new TiXmlElement(meshDimension == 1 ? "ELEMENT" : "EDGE");
2536 std::string tag = meshDimension == 1 ? "S" : "E";
2537
2538 for (auto &i : edges)
2539 {
2540 std::stringstream s;
2541 SegGeomSharedPtr seg = i.second;
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);
2547 }
2548
2549 geomTag->LinkEndChild(edgeTag);
2550}
2551
2552void MeshGraphIOXml::v_WriteTris(TiXmlElement *faceTag, TriGeomMap &tris)
2553{
2554 std::string tag = "T";
2555
2556 for (auto &i : tris)
2557 {
2558 std::stringstream s;
2559 TriGeomSharedPtr tri = i.second;
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);
2565 }
2566}
2567
2568void MeshGraphIOXml::v_WriteQuads(TiXmlElement *faceTag, QuadGeomMap &quads)
2569{
2570 std::string tag = "Q";
2571
2572 for (auto &i : quads)
2573 {
2574 std::stringstream s;
2575 QuadGeomSharedPtr quad = i.second;
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);
2582 }
2583}
2584
2585void MeshGraphIOXml::v_WriteHexs(TiXmlElement *elmtTag, HexGeomMap &hexs)
2586{
2587 std::string tag = "H";
2588
2589 for (auto &i : hexs)
2590 {
2591 std::stringstream s;
2592 HexGeomSharedPtr hex = i.second;
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);
2600 }
2601}
2602
2603void MeshGraphIOXml::v_WritePrisms(TiXmlElement *elmtTag, PrismGeomMap &pris)
2604{
2605 std::string tag = "R";
2606
2607 for (auto &i : pris)
2608 {
2609 std::stringstream s;
2610 PrismGeomSharedPtr prism = i.second;
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);
2618 }
2619}
2620
2621void MeshGraphIOXml::v_WritePyrs(TiXmlElement *elmtTag, PyrGeomMap &pyrs)
2622{
2623 std::string tag = "P";
2624
2625 for (auto &i : pyrs)
2626 {
2627 std::stringstream s;
2628 PyrGeomSharedPtr pyr = i.second;
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);
2635 }
2636}
2637
2638void MeshGraphIOXml::v_WriteTets(TiXmlElement *elmtTag, TetGeomMap &tets)
2639{
2640 std::string tag = "A";
2641
2642 for (auto &i : tets)
2643 {
2644 std::stringstream s;
2645 TetGeomSharedPtr tet = i.second;
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);
2652 }
2653}
2654
2655void MeshGraphIOXml::v_WriteCurves(TiXmlElement *geomTag, CurveMap &edges,
2656 CurveMap &faces)
2657{
2658 TiXmlElement *curveTag = new TiXmlElement("CURVED");
2659 CurveMap::iterator curveIt;
2660 int curveId = 0;
2661
2662 for (curveIt = edges.begin(); curveIt != edges.end(); ++curveIt)
2663 {
2664 CurveSharedPtr curve = curveIt->second;
2665 TiXmlElement *c = new TiXmlElement("E");
2666 std::stringstream s;
2667 s.precision(8);
2668
2669 for (int j = 0; j < curve->m_points.size(); ++j)
2670 {
2671 SpatialDomains::PointGeomSharedPtr p = curve->m_points[j];
2672 s << std::scientific << (*p)(0) << " " << (*p)(1) << " " << (*p)(2)
2673 << " ";
2674 }
2675
2676 c->SetAttribute("ID", curveId++);
2677 c->SetAttribute("EDGEID", curve->m_curveID);
2678 c->SetAttribute("NUMPOINTS", curve->m_points.size());
2679 c->SetAttribute("TYPE", LibUtilities::kPointsTypeStr[curve->m_ptype]);
2680 c->LinkEndChild(new TiXmlText(s.str()));
2681 curveTag->LinkEndChild(c);
2682 }
2683
2684 for (curveIt = faces.begin(); curveIt != faces.end(); ++curveIt)
2685 {
2686 CurveSharedPtr curve = curveIt->second;
2687 TiXmlElement *c = new TiXmlElement("F");
2688 std::stringstream s;
2689 s.precision(8);
2690
2691 for (int j = 0; j < curve->m_points.size(); ++j)
2692 {
2693 SpatialDomains::PointGeomSharedPtr p = curve->m_points[j];
2694 s << std::scientific << (*p)(0) << " " << (*p)(1) << " " << (*p)(2)
2695 << " ";
2696 }
2697
2698 c->SetAttribute("ID", curveId++);
2699 c->SetAttribute("FACEID", curve->m_curveID);
2700 c->SetAttribute("NUMPOINTS", curve->m_points.size());
2701 c->SetAttribute("TYPE", LibUtilities::kPointsTypeStr[curve->m_ptype]);
2702 c->LinkEndChild(new TiXmlText(s.str()));
2703 curveTag->LinkEndChild(c);
2704 }
2705
2706 geomTag->LinkEndChild(curveTag);
2707}
2708
2709void MeshGraphIOXml::WriteComposites(TiXmlElement *geomTag, CompositeMap &comps,
2710 std::map<int, std::string> &compLabels)
2711{
2712 auto &compositesLabels = m_meshGraph->GetCompositesLabels();
2713 TiXmlElement *compTag = new TiXmlElement("COMPOSITE");
2714
2715 for (auto &cIt : comps)
2716 {
2717 if (cIt.second->m_geomVec.size() == 0)
2718 {
2719 continue;
2720 }
2721
2722 TiXmlElement *c = new TiXmlElement("C");
2723 c->SetAttribute("ID", cIt.first);
2724 if (!compositesLabels[cIt.first].empty())
2725 {
2726 c->SetAttribute("NAME", compLabels[cIt.first]);
2727 }
2728 c->LinkEndChild(new TiXmlText(GetCompositeString(cIt.second)));
2729 compTag->LinkEndChild(c);
2730 }
2731
2732 geomTag->LinkEndChild(compTag);
2733}
2734
2735void MeshGraphIOXml::WriteDomain(TiXmlElement *geomTag,
2736 std::map<int, CompositeMap> &domainMap)
2737{
2738 TiXmlElement *domTag = new TiXmlElement("DOMAIN");
2739
2740 std::vector<unsigned int> idxList;
2741 for (auto &domain : domainMap)
2742 {
2743 TiXmlElement *c = new TiXmlElement("D");
2744 idxList.clear();
2745 std::stringstream s;
2746 s << " "
2747 << "C"
2748 << "[";
2749
2750 for (const auto &elem : domain.second)
2751 {
2752 idxList.push_back(elem.first);
2753 }
2754
2755 s << ParseUtils::GenerateSeqString(idxList) << "] ";
2756 c->SetAttribute("ID", domain.first);
2757 c->LinkEndChild(new TiXmlText(s.str()));
2758 domTag->LinkEndChild(c);
2759 }
2760
2761 geomTag->LinkEndChild(domTag);
2762}
2763
2765{
2766 int meshDimension = m_meshGraph->GetMeshDimension();
2767 auto &meshComposites = m_meshGraph->GetComposites();
2768
2769 TiXmlElement *expTag = new TiXmlElement("EXPANSIONS");
2770
2771 for (auto it = meshComposites.begin(); it != meshComposites.end(); it++)
2772 {
2773 if (it->second->m_geomVec[0]->GetShapeDim() == meshDimension)
2774 {
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");
2781
2782 expTag->LinkEndChild(exp);
2783 }
2784 }
2785 root->LinkEndChild(expTag);
2786}
2787
2788/**
2789 * @brief Write out an XML file containing the GEOMETRY block
2790 * representing this MeshGraph instance inside a NEKTAR tag.
2791 */
2793 const std::string &outfilename, bool defaultExp,
2794 const LibUtilities::FieldMetaDataMap &metadata)
2795{
2796 int meshDimension = m_meshGraph->GetMeshDimension();
2797 int spaceDimension = m_meshGraph->GetSpaceDimension();
2798
2799 // Create empty TinyXML document.
2800 TiXmlDocument doc;
2801 TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", "");
2802 doc.LinkEndChild(decl);
2803
2804 TiXmlElement *root = new TiXmlElement("NEKTAR");
2805 doc.LinkEndChild(root);
2806 TiXmlElement *geomTag = new TiXmlElement("GEOMETRY");
2807 root->LinkEndChild(geomTag);
2808
2809 // Add provenance information using FieldIO library.
2811 new LibUtilities::XmlTagWriter(root)),
2812 metadata);
2813
2814 // Update attributes with dimensions.
2815 geomTag->SetAttribute("DIM", meshDimension);
2816 geomTag->SetAttribute("SPACE", spaceDimension);
2817
2818 if (m_session != nullptr && !m_session->GetComm()->IsSerial())
2819 {
2820 geomTag->SetAttribute("PARTITION", m_session->GetComm()->GetRank());
2821 }
2822
2823 // Clear existing elements.
2824 geomTag->Clear();
2825
2826 // Write out informatio
2827 v_WriteVertices(geomTag, m_meshGraph->GetAllPointGeoms());
2828 v_WriteEdges(geomTag, m_meshGraph->GetAllSegGeoms());
2829 if (meshDimension > 1)
2830 {
2831 TiXmlElement *faceTag =
2832 new TiXmlElement(meshDimension == 2 ? "ELEMENT" : "FACE");
2833
2834 v_WriteTris(faceTag, m_meshGraph->GetAllTriGeoms());
2835 v_WriteQuads(faceTag, m_meshGraph->GetAllQuadGeoms());
2836 geomTag->LinkEndChild(faceTag);
2837 }
2838 if (meshDimension > 2)
2839 {
2840 TiXmlElement *elmtTag = new TiXmlElement("ELEMENT");
2841
2842 v_WriteHexs(elmtTag, m_meshGraph->GetAllHexGeoms());
2843 v_WritePyrs(elmtTag, m_meshGraph->GetAllPyrGeoms());
2844 v_WritePrisms(elmtTag, m_meshGraph->GetAllPrismGeoms());
2845 v_WriteTets(elmtTag, m_meshGraph->GetAllTetGeoms());
2846
2847 geomTag->LinkEndChild(elmtTag);
2848 }
2849 v_WriteCurves(geomTag, m_meshGraph->GetCurvedEdges(),
2850 m_meshGraph->GetCurvedFaces());
2851 WriteComposites(geomTag, m_meshGraph->GetComposites(),
2852 m_meshGraph->GetCompositesLabels());
2853 WriteDomain(geomTag, m_meshGraph->GetDomain());
2854
2855 if (defaultExp)
2856 {
2858 }
2859
2860 auto &movement = m_meshGraph->GetMovement();
2861 if (movement)
2862 {
2863 movement->WriteMovement(root);
2864 }
2865
2866 // Save file.
2867 doc.SaveFile(outfilename);
2868}
2869
2871 std::string outname, std::vector<std::set<unsigned int>> elements,
2872 std::vector<unsigned int> partitions)
2873{
2874 // so in theory this function is used by the mesh partitioner
2875 // giving instructions on how to write out a paritioned mesh.
2876 // the theory goes that the elements stored in the meshgraph are the
2877 // "whole" mesh so based on the information from the elmements list
2878 // we can filter the mesh entities and write some individual files
2879 // hopefully
2880
2881 // this is xml so we are going to write a directory with lots of
2882 // xml files
2883 std::string dirname = outname + "_xml";
2884 fs::path pdirname(dirname);
2885
2886 if (!fs::is_directory(dirname))
2887 {
2888 fs::create_directory(dirname);
2889 }
2890
2891 ASSERTL0(elements.size() == partitions.size(),
2892 "error in partitioned information");
2893
2894 for (int i = 0; i < partitions.size(); i++)
2895 {
2896 TiXmlDocument doc;
2897 TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", "");
2898 doc.LinkEndChild(decl);
2899
2900 TiXmlElement *root = doc.FirstChildElement("NEKTAR");
2901 TiXmlElement *geomTag;
2902
2903 // Try to find existing NEKTAR tag.
2904 if (!root)
2905 {
2906 root = new TiXmlElement("NEKTAR");
2907 doc.LinkEndChild(root);
2908
2909 geomTag = new TiXmlElement("GEOMETRY");
2910 root->LinkEndChild(geomTag);
2911 }
2912 else
2913 {
2914 // Try to find existing GEOMETRY tag.
2915 geomTag = root->FirstChildElement("GEOMETRY");
2916
2917 if (!geomTag)
2918 {
2919 geomTag = new TiXmlElement("GEOMETRY");
2920 root->LinkEndChild(geomTag);
2921 }
2922 }
2923
2924 int meshDimension = m_meshGraph->GetMeshDimension();
2925 int spaceDimension = m_meshGraph->GetSpaceDimension();
2926
2927 geomTag->SetAttribute("DIM", meshDimension);
2928 geomTag->SetAttribute("SPACE", spaceDimension);
2929 geomTag->SetAttribute("PARTITION", partitions[i]);
2930
2931 // Add Mesh //
2932 // Get the elements
2933 auto &vertSet = m_meshGraph->GetAllPointGeoms();
2934 auto &segGeoms = m_meshGraph->GetAllSegGeoms();
2935 auto &triGeoms = m_meshGraph->GetAllTriGeoms();
2936 auto &quadGeoms = m_meshGraph->GetAllQuadGeoms();
2937 auto &tetGeoms = m_meshGraph->GetAllTetGeoms();
2938 auto &pyrGeoms = m_meshGraph->GetAllPyrGeoms();
2939 auto &prismGeoms = m_meshGraph->GetAllPrismGeoms();
2940 auto &hexGeoms = m_meshGraph->GetAllHexGeoms();
2941 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
2942 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
2943 auto &meshComposites = m_meshGraph->GetComposites();
2944 auto &globalDomain = m_meshGraph->GetDomain();
2945
2946 HexGeomMap localHex;
2947 PyrGeomMap localPyr;
2948 PrismGeomMap localPrism;
2949 TetGeomMap localTet;
2950 TriGeomMap localTri;
2951 QuadGeomMap localQuad;
2952 SegGeomMap localEdge;
2953 PointGeomMap localVert;
2954 CurveMap localCurveEdge;
2955 CurveMap localCurveFace;
2956
2957 std::vector<std::set<unsigned int>> entityIds(4);
2958 entityIds[meshDimension] = elements[i];
2959
2960 switch (meshDimension)
2961 {
2962 case 3:
2963 {
2964 for (auto &j : entityIds[3])
2965 {
2967 if (hexGeoms.count(j))
2968 {
2969 g = hexGeoms[j];
2970 localHex[j] = hexGeoms[j];
2971 }
2972 else if (pyrGeoms.count(j))
2973 {
2974 g = pyrGeoms[j];
2975 localPyr[j] = pyrGeoms[j];
2976 }
2977 else if (prismGeoms.count(j))
2978 {
2979 g = prismGeoms[j];
2980 localPrism[j] = prismGeoms[j];
2981 }
2982 else if (tetGeoms.count(j))
2983 {
2984 g = tetGeoms[j];
2985 localTet[j] = tetGeoms[j];
2986 }
2987 else
2988 {
2989 ASSERTL0(false, "element in partition not found");
2990 }
2991
2992 for (int k = 0; k < g->GetNumFaces(); k++)
2993 {
2994 entityIds[2].insert(g->GetFid(k));
2995 }
2996 for (int k = 0; k < g->GetNumEdges(); k++)
2997 {
2998 entityIds[1].insert(g->GetEid(k));
2999 }
3000 for (int k = 0; k < g->GetNumVerts(); k++)
3001 {
3002 entityIds[0].insert(g->GetVid(k));
3003 }
3004 }
3005 }
3006 break;
3007 case 2:
3008 {
3009 for (auto &j : entityIds[2])
3010 {
3012 if (triGeoms.count(j))
3013 {
3014 g = triGeoms[j];
3015 localTri[j] = triGeoms[j];
3016 }
3017 else if (quadGeoms.count(j))
3018 {
3019 g = quadGeoms[j];
3020 localQuad[j] = quadGeoms[j];
3021 }
3022 else
3023 {
3024 ASSERTL0(false, "element in partition not found");
3025 }
3026
3027 for (int k = 0; k < g->GetNumEdges(); k++)
3028 {
3029 entityIds[1].insert(g->GetEid(k));
3030 }
3031 for (int k = 0; k < g->GetNumVerts(); k++)
3032 {
3033 entityIds[0].insert(g->GetVid(k));
3034 }
3035 }
3036 }
3037 break;
3038 case 1:
3039 {
3040 for (auto &j : entityIds[1])
3041 {
3043 if (segGeoms.count(j))
3044 {
3045 g = segGeoms[j];
3046 localEdge[j] = segGeoms[j];
3047 }
3048 else
3049 {
3050 ASSERTL0(false, "element in partition not found");
3051 }
3052
3053 for (int k = 0; k < g->GetNumVerts(); k++)
3054 {
3055 entityIds[0].insert(g->GetVid(k));
3056 }
3057 }
3058 }
3059 }
3060
3061 if (meshDimension > 2)
3062 {
3063 for (auto &j : entityIds[2])
3064 {
3065 if (triGeoms.count(j))
3066 {
3067 localTri[j] = triGeoms[j];
3068 }
3069 else if (quadGeoms.count(j))
3070 {
3071 localQuad[j] = quadGeoms[j];
3072 }
3073 else
3074 {
3075 ASSERTL0(false, "face not found");
3076 }
3077 }
3078 }
3079
3080 if (meshDimension > 1)
3081 {
3082 for (auto &j : entityIds[1])
3083 {
3084 if (segGeoms.count(j))
3085 {
3086 localEdge[j] = segGeoms[j];
3087 }
3088 else
3089 {
3090 ASSERTL0(false, "edge not found");
3091 }
3092 }
3093 }
3094
3095 for (auto &j : entityIds[0])
3096 {
3097 if (vertSet.count(j))
3098 {
3099 localVert[j] = vertSet[j];
3100 }
3101 else
3102 {
3103 ASSERTL0(false, "vert not found");
3104 }
3105 }
3106
3107 v_WriteVertices(geomTag, localVert);
3108 v_WriteEdges(geomTag, localEdge);
3109 if (meshDimension > 1)
3110 {
3111 TiXmlElement *faceTag =
3112 new TiXmlElement(meshDimension == 2 ? "ELEMENT" : "FACE");
3113
3114 v_WriteTris(faceTag, localTri);
3115 v_WriteQuads(faceTag, localQuad);
3116 geomTag->LinkEndChild(faceTag);
3117 }
3118 if (meshDimension > 2)
3119 {
3120 TiXmlElement *elmtTag = new TiXmlElement("ELEMENT");
3121
3122 v_WriteHexs(elmtTag, localHex);
3123 v_WritePyrs(elmtTag, localPyr);
3124 v_WritePrisms(elmtTag, localPrism);
3125 v_WriteTets(elmtTag, localTet);
3126
3127 geomTag->LinkEndChild(elmtTag);
3128 }
3129
3130 for (auto &j : localTri)
3131 {
3132 if (curvedFaces.count(j.first))
3133 {
3134 localCurveFace[j.first] = curvedFaces[j.first];
3135 }
3136 }
3137 for (auto &j : localQuad)
3138 {
3139 if (curvedFaces.count(j.first))
3140 {
3141 localCurveFace[j.first] = curvedFaces[j.first];
3142 }
3143 }
3144 for (auto &j : localEdge)
3145 {
3146 if (curvedEdges.count(j.first))
3147 {
3148 localCurveEdge[j.first] = curvedEdges[j.first];
3149 }
3150 }
3151
3152 v_WriteCurves(geomTag, localCurveEdge, localCurveFace);
3153
3154 CompositeMap localComp;
3155 std::map<int, std::string> localCompLabels;
3156
3157 for (auto &j : meshComposites)
3158 {
3160 int dim = j.second->m_geomVec[0]->GetShapeDim();
3161
3162 for (int k = 0; k < j.second->m_geomVec.size(); k++)
3163 {
3164 if (entityIds[dim].count(j.second->m_geomVec[k]->GetGlobalID()))
3165 {
3166 comp->m_geomVec.push_back(j.second->m_geomVec[k]);
3167 }
3168 }
3169
3170 if (comp->m_geomVec.size())
3171 {
3172 localComp[j.first] = comp;
3173 auto compositesLabels = m_meshGraph->GetCompositesLabels();
3174 if (!compositesLabels[j.first].empty())
3175 {
3176 localCompLabels[j.first] = compositesLabels[j.first];
3177 }
3178 }
3179 }
3180
3181 WriteComposites(geomTag, localComp, localCompLabels);
3182
3183 std::map<int, CompositeMap> domain;
3184 for (auto &j : localComp)
3185 {
3186 for (auto &dom : globalDomain)
3187 {
3188 for (auto &dIt : dom.second)
3189 {
3190 if (j.first == dIt.first)
3191 {
3192 domain[dom.first][j.first] = j.second;
3193 break;
3194 }
3195 }
3196 }
3197 }
3198
3199 WriteDomain(geomTag, domain);
3200
3201 if (m_session->DefinesElement("NEKTAR/CONDITIONS"))
3202 {
3203 std::set<int> vBndRegionIdList;
3204 TiXmlElement *vConditions =
3205 new TiXmlElement(*m_session->GetElement("Nektar/Conditions"));
3206 TiXmlElement *vBndRegions =
3207 vConditions->FirstChildElement("BOUNDARYREGIONS");
3208 // Use fine-level for mesh partition (Parallel-in-Time)
3210 TiXmlElement *vBndConditions =
3211 vConditions->FirstChildElement("BOUNDARYCONDITIONS");
3212 // Use fine-level for mesh partition (Parallel-in-Time)
3214 0);
3215 TiXmlElement *vItem;
3216
3217 if (vBndRegions)
3218 {
3219 // Check for parallel-in-time
3220 bool multiLevel =
3221 vConditions->FirstChildElement("BOUNDARYREGIONS")
3222 ->FirstChildElement("TIMELEVEL") != nullptr;
3223
3224 TiXmlElement *vNewBndRegions =
3225 multiLevel ? new TiXmlElement("TIMELEVEL")
3226 : new TiXmlElement("BOUNDARYREGIONS");
3227 vItem = vBndRegions->FirstChildElement();
3228 auto graph_bndRegOrder = m_meshGraph->GetBndRegionOrdering();
3229 while (vItem)
3230 {
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;
3239 ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
3240
3241 std::vector<unsigned int> idxList;
3242
3243 for (unsigned int i = 0; i < vSeq.size(); ++i)
3244 {
3245 if (localComp.find(vSeq[i]) != localComp.end())
3246 {
3247 idxList.push_back(vSeq[i]);
3248 }
3249 }
3250 int p = atoi(vItem->Attribute("ID"));
3251
3252 std::string vListStr =
3254
3255 if (vListStr.length() == 0)
3256 {
3257 TiXmlElement *tmp = vItem;
3258 vItem = vItem->NextSiblingElement();
3259 vBndRegions->RemoveChild(tmp);
3260 }
3261 else
3262 {
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();
3271 }
3272
3273 // store original bnd region order
3274 m_bndRegOrder[p] = vSeq;
3275 graph_bndRegOrder[p] = vSeq;
3276 }
3277 if (multiLevel)
3278 {
3279 // Use fine-level for mesh partition (Parallel-in-Time)
3280 size_t timeLevel = 0;
3281 while (vBndRegions)
3282 {
3283 vNewBndRegions->SetAttribute("VALUE", timeLevel);
3284 vConditions->FirstChildElement("BOUNDARYREGIONS")
3285 ->ReplaceChild(vBndRegions, *vNewBndRegions);
3286 vBndRegions = vBndRegions->NextSiblingElement();
3287 timeLevel++;
3288 }
3289 }
3290 else
3291 {
3292 vConditions->ReplaceChild(vBndRegions, *vNewBndRegions);
3293 }
3294 }
3295
3296 if (vBndConditions)
3297 {
3298 vItem = vBndConditions->FirstChildElement();
3299 while (vItem)
3300 {
3301 std::set<int>::iterator x;
3302 if ((x = vBndRegionIdList.find(atoi(vItem->Attribute(
3303 "REF")))) != vBndRegionIdList.end())
3304 {
3305 vItem->SetAttribute("REF", *x);
3306 vItem = vItem->NextSiblingElement();
3307 }
3308 else
3309 {
3310 TiXmlElement *tmp = vItem;
3311 vItem = vItem->NextSiblingElement();
3312 vBndConditions->RemoveChild(tmp);
3313 }
3314 }
3315 }
3316 root->LinkEndChild(vConditions);
3317 }
3318
3319 // Distribute other sections of the XML to each process as is.
3320 TiXmlElement *vSrc =
3321 m_session->GetElement("Nektar")->FirstChildElement();
3322 while (vSrc)
3323 {
3324 std::string vName = boost::to_upper_copy(vSrc->ValueStr());
3325 if (vName != "GEOMETRY" && vName != "CONDITIONS")
3326 {
3327 root->LinkEndChild(new TiXmlElement(*vSrc));
3328 }
3329 vSrc = vSrc->NextSiblingElement();
3330 }
3331
3332 // Save Mesh
3333
3334 boost::format pad("P%1$07d.xml");
3335 pad % partitions[i];
3336 fs::path pFilename(pad.str());
3337
3338 fs::path fullpath = pdirname / pFilename;
3339 doc.SaveFile(LibUtilities::PortablePath(fullpath));
3340 }
3341}
3342
3344{
3345 auto &meshComposites = m_meshGraph->GetComposites();
3346 auto &domain = m_meshGraph->GetDomain();
3347
3349
3350 for (auto &c : meshComposites)
3351 {
3352 bool fillComp = true;
3353 for (auto &d : domain[0])
3354 {
3355 if (c.second == d.second)
3356 {
3357 fillComp = false;
3358 }
3359 }
3360 if (fillComp)
3361 {
3362 std::vector<unsigned int> ids;
3363 for (auto &elmt : c.second->m_geomVec)
3364 {
3365 ids.push_back(elmt->GetGlobalID());
3366 }
3367 ret[c.first] = ids;
3368 }
3369 }
3370
3371 return ret;
3372}
3373
3374} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
static void AddInfoTag(TagWriterSharedPtr root, const FieldMetaDataMap &fieldmetadatamap)
Add provenance information to the field metadata map.
Definition: FieldIO.cpp:341
Interpreter class for the evaluation of mathematical expressions.
Definition: Interpreter.h:76
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.
Definition: ParseUtils.h:72
static bool GenerateSeqVector(const std::string &str, std::vector< unsigned int > &out)
Takes a comma-separated compressed string and converts it to entries in a vector.
Definition: ParseUtils.cpp:104
static const int kNfaces
Definition: HexGeom.h:58
static const int kNqfaces
Definition: HexGeom.h:56
static const int kNtfaces
Definition: HexGeom.h:57
LibUtilities::SessionReaderSharedPtr m_session
Definition: MeshGraphIO.h:89
void ReadGeometry(LibUtilities::DomainRangeShPtr rng, bool fillGraph)
Definition: MeshGraphIO.h:73
std::string GetCompositeString(CompositeSharedPtr comp)
Returns a string representation of a composite.
CompositeDescriptor CreateCompositeDescriptor()
std::map< int, MeshEntity > CreateMeshEntities()
Create mesh entities for this graph.
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_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)
void v_PartitionMesh(LibUtilities::SessionReaderSharedPtr session) override
void WriteDefaultExpansion(TiXmlElement *root)
void v_ReadGeometry(LibUtilities::DomainRangeShPtr rng, bool fillGraph) override
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_WritePrisms(TiXmlElement *elmtTag, PrismGeomMap &pris)
void WriteComposites(TiXmlElement *geomTag, CompositeMap &comps, std::map< int, std::string > &compLabels)
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)
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 kNfaces
Definition: PyrGeom.h:55
static const int kNtfaces
Definition: PyrGeom.h:54
static const int kNqfaces
Definition: PyrGeom.h:53
static const int kNedges
Definition: QuadGeom.h:75
static const int kNfaces
Definition: TetGeom.h:56
static const int kNqfaces
Definition: TetGeom.h:54
static const int kNtfaces
Definition: TetGeom.h:55
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:70
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:50
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:52
static std::string PortablePath(const fs::path &path)
create portable path on different platforms for std::filesystem path.
Definition: Filesystem.hpp:66
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< DomainRange > DomainRangeShPtr
Definition: DomainRange.h:64
std::shared_ptr< XmlTagWriter > XmlTagWriterSharedPtr
Definition: FieldIOXml.h:166
static DomainRangeShPtr NullDomainRangeShPtr
Definition: DomainRange.h:65
@ SIZE_PointsType
Length of enum list.
Definition: PointsType.h:95
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:57
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshGraph.h:107
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:45
std::map< int, PyrGeomSharedPtr > PyrGeomMap
Definition: PyrGeom.h:76
std::map< int, QuadGeomSharedPtr > QuadGeomMap
Definition: QuadGeom.h:53
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:82
MeshPartitionFactory & GetMeshPartitionFactory()
std::shared_ptr< Composite > CompositeSharedPtr
Definition: MeshGraph.h:135
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:50
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:59
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:58
std::map< int, TetGeomSharedPtr > TetGeomMap
Definition: TetGeom.h:84
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:75
MeshGraphIOFactory & GetMeshGraphIOFactory()
Definition: MeshGraphIO.cpp:45
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::map< int, PrismGeomSharedPtr > PrismGeomMap
Definition: PrismGeom.h:83
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::shared_ptr< MeshPartition > MeshPartitionSharedPtr
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:51
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:56
std::map< int, HexGeomSharedPtr > HexGeomMap
Definition: HexGeom.h:85
std::map< int, PointGeomSharedPtr > PointGeomMap
Definition: PointGeom.h:51
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:136
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:475
std::vector< double > d(NPUPPER *NPUPPER)
std::vector< double > q(NPUPPER *NPUPPER)
double NekDouble