Nektar++
MeshGraph3D.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MeshGraph3D.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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description:
33 //
34 //
35 ////////////////////////////////////////////////////////////////////////////////
36 
38 #include <SpatialDomains/TriGeom.h>
40 #include <tinyxml.h>
41 
42 namespace Nektar
43 {
44  namespace SpatialDomains
45  {
47  {
48  }
49 
51  const DomainRangeShPtr &rng)
52  : MeshGraph(pSession,rng)
53  {
54  ReadGeometry(pSession->GetDocument());
55  ReadExpansions(pSession->GetDocument());
56  }
57 
59  {
60  }
61 
62  void MeshGraph3D::ReadGeometry(const std::string &infilename)
63  {
64  TiXmlDocument doc(infilename);
65  bool loadOkay = doc.LoadFile();
66 
67  std::stringstream errstr;
68  errstr << "Unable to load file: " << infilename << "\n";
69  errstr << doc.ErrorDesc() << " (Line " << doc.ErrorRow()
70  << ", Column " << doc.ErrorCol() << ")";
71  ASSERTL0(loadOkay, errstr.str());
72 
73  ReadGeometry(doc);
74  }
75 
76  // \brief Read segments (and general MeshGraph) given TiXmlDocument.
77  void MeshGraph3D::ReadGeometry(TiXmlDocument &doc)
78  {
79  // Read mesh first
81  TiXmlHandle docHandle(&doc);
82 
83  TiXmlElement* mesh = NULL;
84 
85  /// Look for all geometry related data in GEOMETRY block.
86  mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
87 
88  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
89 
90  ReadCurves(doc);
91  ReadEdges(doc);
92  ReadFaces(doc);
93  ReadElements(doc);
94  ReadComposites(doc);
95  ReadDomain(doc);
96  }
97 
98  void MeshGraph3D::ReadEdges(TiXmlDocument &doc)
99  {
100  /// We know we have it since we made it this far.
101  TiXmlHandle docHandle(&doc);
102  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
103  TiXmlElement* field = NULL;
104 
105  /// Look for elements in ELEMENT block.
106  field = mesh->FirstChildElement("EDGE");
107 
108  ASSERTL0(field, "Unable to find EDGE tag in file.");
109 
110  /// All elements are of the form: "<E ID="#"> ... </E>", with
111  /// ? being the element type.
112  /// Read the ID field first.
113  TiXmlElement *edge = field->FirstChildElement("E");
114 
115  /// Since all edge data is one big text block, we need to accumulate
116  /// all TINYXML_TEXT data and then parse it. This approach effectively skips
117  /// all comments or other node types since we only care about the
118  /// edge list. We cannot handle missing edge numbers as we could
119  /// with missing element numbers due to the text block format.
120  std::string edgeStr;
121  int indx;
123 
124  while(edge)
125  {
126  int err = edge->QueryIntAttribute("ID",&indx);
127  ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
128 
129  TiXmlNode *child = edge->FirstChild();
130  edgeStr.clear();
131  if (child->Type() == TiXmlNode::TINYXML_TEXT)
132  {
133  edgeStr += child->ToText()->ValueStr();
134  }
135 
136  /// Now parse out the edges, three fields at a time.
137  int vertex1, vertex2;
138  std::istringstream edgeDataStrm(edgeStr.c_str());
139 
140  try
141  {
142  while (!edgeDataStrm.fail())
143  {
144  edgeDataStrm >> vertex1 >> vertex2;
145 
146  // Must check after the read because we may be at the end and not know it.
147  // If we are at the end we will add a duplicate of the last entry if we
148  // don't check here.
149  if (!edgeDataStrm.fail())
150  {
151  PointGeomSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
152  SegGeomSharedPtr edge;
153  it = m_curvedEdges.find(indx);
154 
155  if (it == m_curvedEdges.end())
156  {
158  }
159  else
160  {
161  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices, it->second);
162  }
163 
164  m_segGeoms[indx] = edge;
165  }
166  }
167  }
168  catch(...)
169  {
170  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
171  }
172 
173  edge = edge->NextSiblingElement("E");
174  }
175  }
176 
177  void MeshGraph3D::ReadFaces(TiXmlDocument &doc)
178  {
179  /// We know we have it since we made it this far.
180  TiXmlHandle docHandle(&doc);
181  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
182  TiXmlElement* field = NULL;
183 
184  /// Look for elements in FACE block.
185  field = mesh->FirstChildElement("FACE");
186 
187  ASSERTL0(field, "Unable to find FACE tag in file.");
188 
189  /// All faces are of the form: "<? ID="#"> ... </?>", with
190  /// ? being an element type (either Q or T).
191 
192  TiXmlElement *element = field->FirstChildElement();
194 
195  while (element)
196  {
197  std::string elementType(element->ValueStr());
198 
199  ASSERTL0(elementType == "Q" || elementType == "T",
200  (std::string("Unknown 3D face type: ") + elementType).c_str());
201 
202  /// Read id attribute.
203  int indx;
204  int err = element->QueryIntAttribute("ID", &indx);
205  ASSERTL0(err == TIXML_SUCCESS, "Unable to read face attribute ID.");
206 
207  /// See if this face has curves.
208  it = m_curvedFaces.find(indx);
209 
210  /// Read text element description.
211  TiXmlNode* elementChild = element->FirstChild();
212  std::string elementStr;
213  while(elementChild)
214  {
215  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
216  {
217  elementStr += elementChild->ToText()->ValueStr();
218  }
219  elementChild = elementChild->NextSibling();
220  }
221 
222  ASSERTL0(!elementStr.empty(), "Unable to read face description body.");
223 
224  /// Parse out the element components corresponding to type of element.
225  if (elementType == "T")
226  {
227  // Read three edge numbers
228  int edge1, edge2, edge3;
229  std::istringstream elementDataStrm(elementStr.c_str());
230 
231  try
232  {
233  elementDataStrm >> edge1;
234  elementDataStrm >> edge2;
235  elementDataStrm >> edge3;
236 
237  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read face data for TRIANGLE: ") + elementStr).c_str());
238 
239  /// Create a TriGeom to hold the new definition.
241  {
242  GetSegGeom(edge1),
243  GetSegGeom(edge2),
244  GetSegGeom(edge3)
245  };
246 
248  {
249  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
250  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
251  SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
252  };
253 
254  TriGeomSharedPtr trigeom;
255 
256  if (it == m_curvedFaces.end())
257  {
258  trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, edgeorient);
259  }
260  else
261  {
262  trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, edgeorient, it->second);
263  }
264 
265  trigeom->SetGlobalID(indx);
266 
267  m_triGeoms[indx] = trigeom;
268  }
269  catch(...)
270  {
272  (std::string("Unable to read face data for TRIANGLE: ") + elementStr).c_str());
273  }
274  }
275  else if (elementType == "Q")
276  {
277  // Read four edge numbers
278  int edge1, edge2, edge3, edge4;
279  std::istringstream elementDataStrm(elementStr.c_str());
280 
281  try
282  {
283  elementDataStrm >> edge1;
284  elementDataStrm >> edge2;
285  elementDataStrm >> edge3;
286  elementDataStrm >> edge4;
287 
288  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read face data for QUAD: ") + elementStr).c_str());
289 
290  /// Create a QuadGeom to hold the new definition.
292  {GetSegGeom(edge1),GetSegGeom(edge2),
293  GetSegGeom(edge3),GetSegGeom(edge4)};
294 
296  {
297  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
298  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
299  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
300  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
301  };
302 
303  QuadGeomSharedPtr quadgeom;
304 
305  if (it == m_curvedEdges.end())
306  {
307  quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient);
308  }
309  else
310  {
311  quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient, it->second);
312  }
313  quadgeom->SetGlobalID(indx);
314 
315  m_quadGeoms[indx] = quadgeom;
316 
317  }
318  catch(...)
319  {
320  NEKERROR(ErrorUtil::efatal,(std::string("Unable to read face data for QUAD: ") + elementStr).c_str());
321  }
322  }
323 
324  /// Keep looking
325  element = element->NextSiblingElement();
326  }
327  }
328 
329  void MeshGraph3D::ReadElements(TiXmlDocument &doc)
330  {
331  /// We know we have it since we made it this far.
332  TiXmlHandle docHandle(&doc);
333  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
334  TiXmlElement* field = NULL;
335 
336  /// Look for elements in ELEMENT block.
337  field = mesh->FirstChildElement("ELEMENT");
338 
339  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
340 
341  int nextElementNumber = -1;
342 
343  /// All elements are of the form: "<? ID="#"> ... </?>", with
344  /// ? being the element type.
345 
346  TiXmlElement *element = field->FirstChildElement();
347 
348  while (element)
349  {
350  std::string elementType(element->ValueStr());
351 
352  //A - tet, P - pyramid, R - prism, H - hex
353  ASSERTL0(elementType == "A" || elementType == "P" || elementType == "R" || elementType == "H",
354  (std::string("Unknown 3D element type: ") + elementType).c_str());
355 
356  /// These should be ordered.
357  nextElementNumber++;
358 
359  /// Read id attribute.
360  int indx;
361  int err = element->QueryIntAttribute("ID", &indx);
362  ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
363 // ASSERTL0(indx == nextElementNumber, "Element IDs must begin with zero and be sequential.");
364 
365  /// Read text element description.
366  TiXmlNode* elementChild = element->FirstChild();
367  std::string elementStr;
368  while(elementChild)
369  {
370  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
371  {
372  elementStr += elementChild->ToText()->ValueStr();
373  }
374  elementChild = elementChild->NextSibling();
375  }
376 
377  ASSERTL0(!elementStr.empty(), "Unable to read element description body.");
378 
379  std::istringstream elementDataStrm(elementStr.c_str());
380 
381  /// Parse out the element components corresponding to type of element.
382 
383  // Tetrahedral
384  if (elementType == "A")
385  {
386  try
387  {
388  /// Create arrays for the tri and quad faces.
389  const int kNfaces = TetGeom::kNfaces;
390  const int kNtfaces = TetGeom::kNtfaces;
391  const int kNqfaces = TetGeom::kNqfaces;
392  TriGeomSharedPtr tfaces[kNtfaces];
393  //QuadGeomSharedPtr qfaces[kNqfaces];
394  int Ntfaces = 0;
395  int Nqfaces = 0;
396 
397  /// Fill the arrays and make sure there aren't too many faces.
398  std::stringstream errorstring;
399  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
400  for (int i = 0; i < kNfaces; i++)
401  {
402  int faceID;
403  elementDataStrm >> faceID;
404  Geometry2DSharedPtr face = GetGeometry2D(faceID);
405  if (face == Geometry2DSharedPtr() ||
406  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
407  {
408  std::stringstream errorstring;
409  errorstring << "Element " << indx << " has invalid face: " << faceID;
410  ASSERTL0(false, errorstring.str().c_str());
411  }
412  else if (face->GetShapeType() == LibUtilities::eTriangle)
413  {
414  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
415  tfaces[Ntfaces++] = boost::static_pointer_cast<TriGeom>(face);
416  }
417  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
418  {
419  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
420  }
421  }
422 
423  /// Make sure all of the face indicies could be read, and that there weren't too few.
424  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TETRAHEDRON: ") + elementStr).c_str());
425  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
426  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
427 
429  tetgeom->SetGlobalID(indx);
430 
431  m_tetGeoms[indx] = tetgeom;
432  PopulateFaceToElMap(tetgeom, kNfaces);
433  }
434  catch(...)
435  {
437  (std::string("Unable to read element data for TETRAHEDRON: ") + elementStr).c_str());
438  }
439  }
440  // Pyramid
441  else if (elementType == "P")
442  {
443  try
444  {
445  /// Create arrays for the tri and quad faces.
446  const int kNfaces = PyrGeom::kNfaces;
447  const int kNtfaces = PyrGeom::kNtfaces;
448  const int kNqfaces = PyrGeom::kNqfaces;
449  Geometry2DSharedPtr faces[kNfaces];
450  int Nfaces = 0;
451  int Ntfaces = 0;
452  int Nqfaces = 0;
453 
454  /// Fill the arrays and make sure there aren't too many faces.
455  std::stringstream errorstring;
456  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
457  for (int i = 0; i < kNfaces; i++)
458  {
459  int faceID;
460  elementDataStrm >> faceID;
461  Geometry2DSharedPtr face = GetGeometry2D(faceID);
462  if (face == Geometry2DSharedPtr() ||
463  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
464  {
465  std::stringstream errorstring;
466  errorstring << "Element " << indx << " has invalid face: " << faceID;
467  ASSERTL0(false, errorstring.str().c_str());
468  }
469  else if (face->GetShapeType() == LibUtilities::eTriangle)
470  {
471  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
472  faces[Nfaces++] = boost::static_pointer_cast<TriGeom>(face);
473  Ntfaces++;
474  }
475  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
476  {
477  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
478  faces[Nfaces++] = boost::static_pointer_cast<QuadGeom>(face);
479  Nqfaces++;
480  }
481  }
482 
483  /// Make sure all of the face indicies could be read, and that there weren't too few.
484  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for PYRAMID: ") + elementStr).c_str());
485  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
486  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
487 
489  pyrgeom->SetGlobalID(indx);
490 
491  m_pyrGeoms[indx] = pyrgeom;
492  PopulateFaceToElMap(pyrgeom, kNfaces);
493  }
494  catch(...)
495  {
497  (std::string("Unable to read element data for PYRAMID: ") + elementStr).c_str());
498  }
499  }
500  // Prism
501  else if (elementType == "R")
502  {
503  try
504  {
505  /// Create arrays for the tri and quad faces.
506  const int kNfaces = PrismGeom::kNfaces;
507  const int kNtfaces = PrismGeom::kNtfaces;
508  const int kNqfaces = PrismGeom::kNqfaces;
509  Geometry2DSharedPtr faces[kNfaces];
510  int Ntfaces = 0;
511  int Nqfaces = 0;
512  int Nfaces = 0;
513 
514  /// Fill the arrays and make sure there aren't too many faces.
515  std::stringstream errorstring;
516  errorstring << "Element " << indx << " must have "
517  << kNtfaces << " triangle face(s), and "
518  << kNqfaces << " quadrilateral face(s).";
519 
520  for (int i = 0; i < kNfaces; i++)
521  {
522  int faceID;
523  elementDataStrm >> faceID;
524  Geometry2DSharedPtr face = GetGeometry2D(faceID);
525  if (face == Geometry2DSharedPtr() ||
526  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
527  {
528  std::stringstream errorstring;
529  errorstring << "Element " << indx << " has invalid face: " << faceID;
530  ASSERTL0(false, errorstring.str().c_str());
531  }
532  else if (face->GetShapeType() == LibUtilities::eTriangle)
533  {
534  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
535  faces[Nfaces++] = boost::static_pointer_cast<TriGeom>(face);
536  Ntfaces++;
537  }
538  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
539  {
540  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
541  faces[Nfaces++] = boost::static_pointer_cast<QuadGeom>(face);
542  Nqfaces++;
543  }
544  }
545 
546  /// Make sure all of the face indicies could be read, and that there weren't too few.
547  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for PRISM: ") + elementStr).c_str());
548  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
549  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
550 
552  prismgeom->SetGlobalID(indx);
553 
554  m_prismGeoms[indx] = prismgeom;
555  PopulateFaceToElMap(prismgeom, kNfaces);
556  }
557  catch(...)
558  {
560  (std::string("Unable to read element data for PRISM: ") + elementStr).c_str());
561  }
562  }
563  // Hexahedral
564  else if (elementType == "H")
565  {
566  try
567  {
568  /// Create arrays for the tri and quad faces.
569  const int kNfaces = HexGeom::kNfaces;
570  const int kNtfaces = HexGeom::kNtfaces;
571  const int kNqfaces = HexGeom::kNqfaces;
572  //TriGeomSharedPtr tfaces[kNtfaces];
573  QuadGeomSharedPtr qfaces[kNqfaces];
574  int Ntfaces = 0;
575  int Nqfaces = 0;
576 
577  /// Fill the arrays and make sure there aren't too many faces.
578  std::stringstream errorstring;
579  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
580  for (int i = 0; i < kNfaces; i++)
581  {
582  int faceID;
583  elementDataStrm >> faceID;
584  Geometry2DSharedPtr face = GetGeometry2D(faceID);
585  if (face == Geometry2DSharedPtr() ||
586  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
587  {
588  std::stringstream errorstring;
589  errorstring << "Element " << indx << " has invalid face: " << faceID;
590  ASSERTL0(false, errorstring.str().c_str());
591  }
592  else if (face->GetShapeType() == LibUtilities::eTriangle)
593  {
594  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
595  //tfaces[Ntfaces++] = boost::static_pointer_cast<TriGeom>(face);
596  }
597  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
598  {
599  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
600  qfaces[Nqfaces++] = boost::static_pointer_cast<QuadGeom>(face);
601  }
602  }
603 
604  /// Make sure all of the face indicies could be read, and that there weren't too few.
605  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for HEXAHEDRAL: ") + elementStr).c_str());
606  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
607  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
608 
610  hexgeom->SetGlobalID(indx);
611 
612  m_hexGeoms[indx] = hexgeom;
613  PopulateFaceToElMap(hexgeom, kNfaces);
614  }
615  catch(...)
616  {
618  (std::string("Unable to read element data for HEXAHEDRAL: ") + elementStr).c_str());
619  }
620  }
621 
622  /// Keep looking
623  element = element->NextSiblingElement();
624  }
625  }
626 
627  void MeshGraph3D::ReadComposites(TiXmlDocument &doc)
628  {
629  TiXmlHandle docHandle(&doc);
630 
631  /// We know we have it since we made it this far.
632  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
633  TiXmlElement* field = NULL;
634 
635  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
636 
637  /// Look for elements in ELEMENT block.
638  field = mesh->FirstChildElement("COMPOSITE");
639 
640  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
641 
642  int nextCompositeNumber = -1;
643 
644  /// All elements are of the form: "<C ID = "N"> ... </C>".
645 
646  /// Read the ID field first.
647  TiXmlElement *composite = field->FirstChildElement("C");
648 
649  while (composite)
650  {
651  nextCompositeNumber++;
652 
653  int indx;
654  int err = composite->QueryIntAttribute("ID", &indx);
655  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
656 // ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with zero and be sequential.");
657 
658  TiXmlNode* compositeChild = composite->FirstChild();
659  // This is primarily to skip comments that may be present.
660  // Comments appear as nodes just like elements.
661  // We are specifically looking for text in the body
662  // of the definition.
663  while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
664  {
665  compositeChild = compositeChild->NextSibling();
666  }
667 
668  ASSERTL0(compositeChild, "Unable to read composite definition body.");
669  std::string compositeStr = compositeChild->ToText()->ValueStr();
670 
671  /// Parse out the element components corresponding to type of element.
672 
673  std::istringstream compositeDataStrm(compositeStr.c_str());
674 
675  try
676  {
677  bool first = true;
678  std::string prevCompositeElementStr;
679 
680  while (!compositeDataStrm.fail())
681  {
682  std::string compositeElementStr;
683  compositeDataStrm >> compositeElementStr;
684 
685  if (!compositeDataStrm.fail())
686  {
687  if (first)
688  {
689  first = false;
690 
692  m_meshComposites[indx] = curVector;
693  }
694 
695  if (compositeElementStr.length() > 0)
696  {
697  ResolveGeomRef(prevCompositeElementStr, compositeElementStr, m_meshComposites[indx]);
698  }
699  prevCompositeElementStr = compositeElementStr;
700  }
701  }
702  }
703  catch(...)
704  {
706  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
707  }
708 
709  /// Keep looking
710  composite = composite->NextSiblingElement("C");
711  }
712  }
713 
714 
716  {
717  SegGeomSharedPtr returnval;
718  SegGeomMap::iterator x = m_segGeoms.find(eID);
719  ASSERTL0(x != m_segGeoms.end(), "Segment "
720  + boost::lexical_cast<string>(eID) + " not found.");
721  return x->second;
722  };
723 
725  {
726  TriGeomMapIter it1;
727  QuadGeomMapIter it2;
728 
729  it1 = m_triGeoms.find(gID);
730  if (it1 != m_triGeoms.end())
731  return it1->second;
732 
733  it2 = m_quadGeoms.find(gID);
734  if (it2 != m_quadGeoms.end())
735  return it2->second;
736 
737  return Geometry2DSharedPtr();
738  };
739 
740  // Take the string that is the composite reference and find the
741  // pointer to the Geometry object corresponding to it.
742 
743  // The only allowable combinations of previous and current items
744  // are V (0D); E (1D); and T and Q (2D); A (Tet, 3D), P (Pyramid, 3D), R (Prism, 3D), H (Hex, 3D).
745  // Only elements of the same dimension are allowed to be grouped.
746  void MeshGraph3D::ResolveGeomRef(const std::string &prevToken, const std::string &token,
747  Composite& composite)
748  {
749  try
750  {
751  std::istringstream tokenStream(token);
752  std::istringstream prevTokenStream(prevToken);
753 
754  char type;
755  char prevType;
756 
757  tokenStream >> type;
758 
759  std::string::size_type indxBeg = token.find_first_of('[') + 1;
760  std::string::size_type indxEnd = token.find_last_of(']') - 1;
761 
762  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());
763 
764  std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
765 
766  std::vector<unsigned int> seqVector;
768 
769  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
770 
771  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
772 
773  prevTokenStream >> prevType;
774 
775  // All composites must be of the same dimension. This map makes things clean to compare.
776  map<char, int> typeMap;
777  typeMap['V'] = 1; // Vertex
778  typeMap['E'] = 1; // Edge
779  typeMap['T'] = 2; // Triangle
780  typeMap['Q'] = 2; // Quad
781  typeMap['A'] = 3; // Tet
782  typeMap['P'] = 3; // Pyramid
783  typeMap['R'] = 3; // Prism
784  typeMap['H'] = 3; // Hex
785 
786  // Make sure only geoms of the same dimension are combined.
787  bool validSequence = (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
788 
789  ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
790  + type + " and " + prevType + ".").c_str());
791 
792  switch(type)
793  {
794  case 'V': // Vertex
795  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
796  {
797  if (m_vertSet.find(*seqIter) == m_vertSet.end())
798  {
799  char errStr[16] = "";
800  ::sprintf(errStr, "%d", *seqIter);
801  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
802  }
803  else
804  {
805  composite->push_back(m_vertSet[*seqIter]);
806  }
807  }
808  break;
809 
810  case 'E': // Edge
811  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
812  {
813  if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
814  {
815  char errStr[16] = "";
816  ::sprintf(errStr, "%d", *seqIter);
817  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
818  }
819  else
820  {
821  composite->push_back(m_segGeoms[*seqIter]);
822  }
823  }
824  break;
825 
826  case 'F': // Face
827  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
828  {
829  Geometry2DSharedPtr face = GetGeometry2D(*seqIter);
830  if (face == Geometry2DSharedPtr())
831  {
832  char errStr[16] = "";
833  ::sprintf(errStr, "%d", *seqIter);
834  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown face index: ") + errStr).c_str());
835  }
836  else
837  {
838  if(CheckRange(*face))
839  {
840  composite->push_back(face);
841  }
842  }
843  }
844  break;
845 
846  case 'T': // Triangle
847  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
848  {
849  if (m_triGeoms.find(*seqIter) == m_triGeoms.end())
850  {
851  char errStr[16] = "";
852  ::sprintf(errStr, "%d", *seqIter);
853  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
854  }
855  else
856  {
857  if(CheckRange(*m_triGeoms[*seqIter]))
858  {
859  composite->push_back(m_triGeoms[*seqIter]);
860  }
861  }
862  }
863  break;
864 
865  case 'Q': // Quad
866  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
867  {
868  if (m_quadGeoms.find(*seqIter) == m_quadGeoms.end())
869  {
870  char errStr[16] = "";
871  ::sprintf(errStr, "%d", *seqIter);
872  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr).c_str());
873  }
874  else
875  {
876  if(CheckRange(*m_quadGeoms[*seqIter]))
877  {
878  composite->push_back(m_quadGeoms[*seqIter]);
879  }
880  }
881  }
882  break;
883 
884  // Tetrahedron
885  case 'A':
886  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
887  {
888  if (m_tetGeoms.find(*seqIter) == m_tetGeoms.end())
889  {
890  char errStr[16] = "";
891  ::sprintf(errStr, "%d", *seqIter);
892  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown tet index: ") + errStr).c_str());
893  }
894  else
895  {
896  if(CheckRange(*m_tetGeoms[*seqIter]))
897  {
898  composite->push_back(m_tetGeoms[*seqIter]);
899  }
900  }
901  }
902  break;
903 
904  // Pyramid
905  case 'P':
906  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
907  {
908  if (m_pyrGeoms.find(*seqIter) == m_pyrGeoms.end())
909  {
910  char errStr[16] = "";
911  ::sprintf(errStr, "%d", *seqIter);
912  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown pyramid index: ") + errStr).c_str());
913  }
914  else
915  {
916  if(CheckRange(*m_pyrGeoms[*seqIter]))
917  {
918  composite->push_back(m_pyrGeoms[*seqIter]);
919  }
920  }
921  }
922  break;
923 
924  // Prism
925  case 'R':
926  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
927  {
928  if (m_prismGeoms.find(*seqIter) == m_prismGeoms.end())
929  {
930  char errStr[16] = "";
931  ::sprintf(errStr, "%d", *seqIter);
932  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown prism index: ") + errStr).c_str());
933  }
934  else
935  {
936  if(CheckRange(*m_prismGeoms[*seqIter]))
937  {
938  composite->push_back(m_prismGeoms[*seqIter]);
939  }
940  }
941  }
942  break;
943 
944  // Hex
945  case 'H':
946  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
947  {
948  if (m_hexGeoms.find(*seqIter) == m_hexGeoms.end())
949  {
950  char errStr[16] = "";
951  ::sprintf(errStr, "%d", *seqIter);
952  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown hex index: ") + errStr).c_str());
953  }
954  else
955  {
956  if(CheckRange(*m_hexGeoms[*seqIter]))
957  {
958  composite->push_back(m_hexGeoms[*seqIter]);
959  }
960  }
961  }
962  break;
963 
964  default:
965  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
966  }
967  }
968  catch(...)
969  {
970  NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
971  }
972 
973  return;
974  }
975 
976 
978  {
980  m_faceToElMap.find(face->GetGlobalID());
981 
982  ASSERTL0(it != m_faceToElMap.end(), "Unable to find corresponding face!");
983 
984  return it->second;
985  }
986 
987 
988  /**
989  * Retrieve the basis key for a given face direction.
990  */
992  Geometry2DSharedPtr face,
993  const int facedir,
994  const std::string variable)
995  {
996  // Retrieve the list of elements and the associated face index
997  // to which the face geometry belongs.
999 
1000  ASSERTL0(elements->size() > 0, "No elements for the given face."
1001  " Check all elements belong to the domain composite.");
1002 
1003  // Perhaps, a check should be done here to ensure that in case
1004  // elements->size!=1, all elements to which the edge belongs have
1005  // the same type and order of expansion such that no confusion can
1006  // arise.
1007 
1008  // Get the Expansion structure detailing the basis keys used for
1009  // this element.
1010  ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element,
1011  variable);
1012 
1013  ASSERTL0(expansion, "Could not find expansion connected to face "+
1014  boost::lexical_cast<string>(face->GetGlobalID()));
1015 
1016  // Retrieve the geometry object of the element as a Geometry3D.
1017  Geometry3DSharedPtr geom3d =
1018  boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
1019  expansion->m_geomShPtr);
1020 
1021  // Use the geometry of the element to calculate the coordinate
1022  // direction of the element which corresponds to the requested
1023  // coordinate direction of the given face.
1024  int dir = geom3d->GetDir((*elements)[0]->m_FaceIndx, facedir);
1025 
1026  if(face->GetNumVerts() == 3)
1027  {
1028  return StdRegions::EvaluateTriFaceBasisKey(facedir,
1029  expansion->m_basisKeyVector[dir].GetBasisType(),
1030  expansion->m_basisKeyVector[dir].GetNumPoints(),
1031  expansion->m_basisKeyVector[dir].GetNumModes());
1032  }
1033  else
1034  {
1035  return StdRegions::EvaluateQuadFaceBasisKey(facedir,
1036  expansion->m_basisKeyVector[dir].GetBasisType(),
1037  expansion->m_basisKeyVector[dir].GetNumPoints(),
1038  expansion->m_basisKeyVector[dir].GetNumModes());
1039  }
1040 
1041  // Keep things happy by returning a value.
1043  }
1044 
1045 
1046  /**
1047  * @brief Given a 3D geometry object #element, populate the face to
1048  * element map #m_faceToElMap which maps faces to their corresponding
1049  * element(s).
1050  *
1051  * @param element Element to process.
1052  * @param kNfaces Number of faces of #element. Should be removed and
1053  * put into Geometry3D as a virtual member function.
1054  */
1056  {
1057  // Set up face -> element map
1058  for (int i = 0; i < kNfaces; ++i)
1059  {
1060  int faceId = element->GetFace(i)->GetGlobalID();
1061  ElementFaceSharedPtr elementFace =
1063 
1064  elementFace->m_Element = element;
1065  elementFace->m_FaceIndx = i;
1066 
1067  // Search map to see if face already exists.
1069  m_faceToElMap.find(faceId);
1070 
1071  if (it == m_faceToElMap.end())
1072  {
1075  tmp->push_back(elementFace);
1076  m_faceToElMap[faceId] = tmp;
1077  }
1078  else
1079  {
1080  ElementFaceVectorSharedPtr tmp = it->second;
1081  tmp->push_back(elementFace);
1082  }
1083  }
1084  }
1085  }; //end of namespace
1086 }; //end of namespace
boost::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:83
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
boost::shared_ptr< ElementFaceVector > ElementFaceVectorSharedPtr
Definition: MeshGraph.h:135
static const int kNfaces
Definition: PyrGeom.h:60
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void ReadCurves(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1103
int GetDir(const int faceidx, const int facedir) const
Returns the element coordinate direction corresponding to a given face coordinate direction...
Definition: Geometry3D.cpp:88
void ReadFaces(TiXmlDocument &doc)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static const int kNtfaces
Definition: HexGeom.h:63
void PopulateFaceToElMap(Geometry3DSharedPtr element, int kNfaces)
Given a 3D geometry object #element, populate the face to element map m_faceToElMap which maps faces ...
void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph3D.cpp:62
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
PointGeomSharedPtr GetVertex(int id)
Definition: MeshGraph.h:567
boost::unordered_map< int, ElementFaceVectorSharedPtr > m_faceToElMap
Definition: MeshGraph3D.h:220
static const int kNfaces
Definition: HexGeom.h:64
std::map< int, QuadGeomSharedPtr >::iterator QuadGeomMapIter
Definition: QuadGeom.h:58
static const int kNtfaces
Definition: TetGeom.h:59
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Geometry2DSharedPtr GetGeometry2D(int gID)
void ResolveGeomRef(const std::string &prevToken, const std::string &token, Composite &composite)
static const int kNfaces
Definition: TetGeom.h:60
SegGeomSharedPtr GetSegGeom(int eID)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:78
virtual void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph.cpp:230
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:293
std::map< int, TriGeomSharedPtr >::iterator TriGeomMapIter
Definition: TriGeom.h:63
ElementFaceVectorSharedPtr GetElementsFromFace(Geometry2DSharedPtr face)
Return the elements (shared ptrs) that have this face.
void ReadComposites(TiXmlDocument &doc)
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
boost::shared_ptr< DomainRange > DomainRangeShPtr
Definition: MeshGraph.h:154
void ReadElements(TiXmlDocument &doc)
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:111
boost::shared_ptr< ElementFace > ElementFaceSharedPtr
Definition: MeshGraph.h:133
static const int kNqfaces
Definition: TetGeom.h:58
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:98
static const int kNtfaces
Definition: PyrGeom.h:59
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
LibUtilities::BasisKey GetFaceBasisKey(Geometry2DSharedPtr face, const int flag, const std::string variable="DefaultVar")
Return the BasisKey corresponding to a face of an element.
void ReadEdges(TiXmlDocument &doc)
Definition: MeshGraph3D.cpp:98
static const int kNqfaces
Definition: PyrGeom.h:58
3D geometry information
Definition: Geometry3D.h:70
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:170
Base class for a spectral/hp element mesh.
Definition: MeshGraph.h:183
void ReadDomain(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1009
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
bool CheckRange(Geometry2D &geom)
Check if goemetry is in range definition if activated.
Definition: MeshGraph.cpp:1813
ExpansionShPtr GetExpansion(GeometrySharedPtr geom, const std::string variable="DefaultVar")
Definition: MeshGraph.cpp:2136
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
void ReadExpansions(const std::string &infilename)
Read the expansions given the XML file path.
Definition: MeshGraph.cpp:525
static const int kNqfaces
Definition: HexGeom.h:62
Describes the specification for a Basis.
Definition: Basis.h:50
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52