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