Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 i,indx;
122  int nextEdgeNumber = -1;
123 
124  // Curved Edges
125  map<int, int> edge_curved;
126  for(i = 0; i < m_curvedEdges.size(); ++i)
127  {
128  edge_curved[m_curvedEdges[i]->m_curveID] = i;
129  }
130 
131  while(edge)
132  {
133  nextEdgeNumber++;
134 
135  int err = edge->QueryIntAttribute("ID",&indx);
136  ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
137 
138  TiXmlNode *child = edge->FirstChild();
139  edgeStr.clear();
140  if (child->Type() == TiXmlNode::TINYXML_TEXT)
141  {
142  edgeStr += child->ToText()->ValueStr();
143  }
144 
145  /// Now parse out the edges, three fields at a time.
146  int vertex1, vertex2;
147  std::istringstream edgeDataStrm(edgeStr.c_str());
148 
149  try
150  {
151  while (!edgeDataStrm.fail())
152  {
153  edgeDataStrm >> vertex1 >> vertex2;
154 
155  // Must check after the read because we may be at the end and not know it.
156  // If we are at the end we will add a duplicate of the last entry if we
157  // don't check here.
158  if (!edgeDataStrm.fail())
159  {
160  PointGeomSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
161  SegGeomSharedPtr edge;
162 
163  if (edge_curved.count(indx) == 0)
164  {
166  }
167  else
168  {
169  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices, m_curvedEdges[edge_curved.find(indx)->second]);
170  }
171 
172  m_segGeoms[indx] = edge;
173  }
174  }
175  }
176  catch(...)
177  {
178  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
179  }
180 
181  edge = edge->NextSiblingElement("E");
182  }
183  }
184 
185  void MeshGraph3D::ReadFaces(TiXmlDocument &doc)
186  {
187  /// We know we have it since we made it this far.
188  TiXmlHandle docHandle(&doc);
189  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
190  TiXmlElement* field = NULL;
191 
192  /// Look for elements in FACE block.
193  field = mesh->FirstChildElement("FACE");
194 
195  ASSERTL0(field, "Unable to find FACE tag in file.");
196 
197  // Curved faces
198  map<int, int> face_curved;
199  for (int i = 0; i < m_curvedFaces.size(); ++i)
200  {
201  face_curved[m_curvedFaces[i]->m_curveID] = i;
202  }
203 
204  /// All faces are of the form: "<? ID="#"> ... </?>", with
205  /// ? being an element type (either Q or T).
206 
207  TiXmlElement *element = field->FirstChildElement();
208 
209  while (element)
210  {
211  std::string elementType(element->ValueStr());
212 
213  ASSERTL0(elementType == "Q" || elementType == "T",
214  (std::string("Unknown 3D face type: ") + elementType).c_str());
215 
216  /// Read id attribute.
217  int indx;
218  int err = element->QueryIntAttribute("ID", &indx);
219  ASSERTL0(err == TIXML_SUCCESS, "Unable to read face attribute ID.");
220 
221  /// Read text element description.
222  TiXmlNode* elementChild = element->FirstChild();
223  std::string elementStr;
224  while(elementChild)
225  {
226  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
227  {
228  elementStr += elementChild->ToText()->ValueStr();
229  }
230  elementChild = elementChild->NextSibling();
231  }
232 
233  ASSERTL0(!elementStr.empty(), "Unable to read face description body.");
234 
235  /// Parse out the element components corresponding to type of element.
236  if (elementType == "T")
237  {
238  // Read three edge numbers
239  int edge1, edge2, edge3;
240  std::istringstream elementDataStrm(elementStr.c_str());
241 
242  try
243  {
244  elementDataStrm >> edge1;
245  elementDataStrm >> edge2;
246  elementDataStrm >> edge3;
247 
248  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read face data for TRIANGLE: ") + elementStr).c_str());
249 
250  /// Create a TriGeom to hold the new definition.
252  {
253  GetSegGeom(edge1),
254  GetSegGeom(edge2),
255  GetSegGeom(edge3)
256  };
257 
259  {
260  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
261  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
262  SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
263  };
264 
265  TriGeomSharedPtr trigeom;
266 
267  if (face_curved.count(indx) == 0)
268  {
269  trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, edgeorient);
270  }
271  else
272  {
273  trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, edgeorient, m_curvedFaces[face_curved.find(indx)->second]);
274  }
275 
276  trigeom->SetGlobalID(indx);
277 
278  m_triGeoms[indx] = trigeom;
279  }
280  catch(...)
281  {
283  (std::string("Unable to read face data for TRIANGLE: ") + elementStr).c_str());
284  }
285  }
286  else if (elementType == "Q")
287  {
288  // Read four edge numbers
289  int edge1, edge2, edge3, edge4;
290  std::istringstream elementDataStrm(elementStr.c_str());
291 
292  try
293  {
294  elementDataStrm >> edge1;
295  elementDataStrm >> edge2;
296  elementDataStrm >> edge3;
297  elementDataStrm >> edge4;
298 
299  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read face data for QUAD: ") + elementStr).c_str());
300 
301  /// Create a QuadGeom to hold the new definition.
303  {GetSegGeom(edge1),GetSegGeom(edge2),
304  GetSegGeom(edge3),GetSegGeom(edge4)};
305 
307  {
308  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
309  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
310  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
311  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
312  };
313 
314  QuadGeomSharedPtr quadgeom;
315 
316  if (face_curved.count(indx) == 0)
317  {
318  quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient);
319  } else {
320  quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient, m_curvedFaces[face_curved.find(indx)->second]);
321  }
322  quadgeom->SetGlobalID(indx);
323 
324  m_quadGeoms[indx] = quadgeom;
325 
326  }
327  catch(...)
328  {
329  NEKERROR(ErrorUtil::efatal,(std::string("Unable to read face data for QUAD: ") + elementStr).c_str());
330  }
331  }
332 
333  /// Keep looking
334  element = element->NextSiblingElement();
335  }
336  }
337 
338  void MeshGraph3D::ReadElements(TiXmlDocument &doc)
339  {
340  /// We know we have it since we made it this far.
341  TiXmlHandle docHandle(&doc);
342  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
343  TiXmlElement* field = NULL;
344 
345  /// Look for elements in ELEMENT block.
346  field = mesh->FirstChildElement("ELEMENT");
347 
348  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
349 
350  int nextElementNumber = -1;
351 
352  /// All elements are of the form: "<? ID="#"> ... </?>", with
353  /// ? being the element type.
354 
355  TiXmlElement *element = field->FirstChildElement();
356 
357  while (element)
358  {
359  std::string elementType(element->ValueStr());
360 
361  //A - tet, P - pyramid, R - prism, H - hex
362  ASSERTL0(elementType == "A" || elementType == "P" || elementType == "R" || elementType == "H",
363  (std::string("Unknown 3D element type: ") + elementType).c_str());
364 
365  /// These should be ordered.
366  nextElementNumber++;
367 
368  /// Read id attribute.
369  int indx;
370  int err = element->QueryIntAttribute("ID", &indx);
371  ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
372 // ASSERTL0(indx == nextElementNumber, "Element IDs must begin with zero and be sequential.");
373 
374  /// Read text element description.
375  TiXmlNode* elementChild = element->FirstChild();
376  std::string elementStr;
377  while(elementChild)
378  {
379  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
380  {
381  elementStr += elementChild->ToText()->ValueStr();
382  }
383  elementChild = elementChild->NextSibling();
384  }
385 
386  ASSERTL0(!elementStr.empty(), "Unable to read element description body.");
387 
388  std::istringstream elementDataStrm(elementStr.c_str());
389 
390  /// Parse out the element components corresponding to type of element.
391 
392  // Tetrahedral
393  if (elementType == "A")
394  {
395  try
396  {
397  /// Create arrays for the tri and quad faces.
398  const int kNfaces = TetGeom::kNfaces;
399  const int kNtfaces = TetGeom::kNtfaces;
400  const int kNqfaces = TetGeom::kNqfaces;
401  TriGeomSharedPtr tfaces[kNtfaces];
402  //QuadGeomSharedPtr qfaces[kNqfaces];
403  int Ntfaces = 0;
404  int Nqfaces = 0;
405 
406  /// Fill the arrays and make sure there aren't too many faces.
407  std::stringstream errorstring;
408  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
409  for (int i = 0; i < kNfaces; i++)
410  {
411  int faceID;
412  elementDataStrm >> faceID;
413  Geometry2DSharedPtr face = GetGeometry2D(faceID);
414  if (face == Geometry2DSharedPtr() ||
415  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
416  {
417  std::stringstream errorstring;
418  errorstring << "Element " << indx << " has invalid face: " << faceID;
419  ASSERTL0(false, errorstring.str().c_str());
420  }
421  else if (face->GetShapeType() == LibUtilities::eTriangle)
422  {
423  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
424  tfaces[Ntfaces++] = boost::static_pointer_cast<TriGeom>(face);
425  }
426  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
427  {
428  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
429  }
430  }
431 
432  /// Make sure all of the face indicies could be read, and that there weren't too few.
433  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TETRAHEDRON: ") + elementStr).c_str());
434  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
435  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
436 
438  tetgeom->SetGlobalID(indx);
439 
440  m_tetGeoms[indx] = tetgeom;
441  PopulateFaceToElMap(tetgeom, kNfaces);
442  }
443  catch(...)
444  {
446  (std::string("Unable to read element data for TETRAHEDRON: ") + elementStr).c_str());
447  }
448  }
449  // Pyramid
450  else if (elementType == "P")
451  {
452  try
453  {
454  /// Create arrays for the tri and quad faces.
455  const int kNfaces = PyrGeom::kNfaces;
456  const int kNtfaces = PyrGeom::kNtfaces;
457  const int kNqfaces = PyrGeom::kNqfaces;
458  Geometry2DSharedPtr faces[kNfaces];
459  int Nfaces = 0;
460  int Ntfaces = 0;
461  int Nqfaces = 0;
462 
463  /// Fill the arrays and make sure there aren't too many faces.
464  std::stringstream errorstring;
465  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
466  for (int i = 0; i < kNfaces; i++)
467  {
468  int faceID;
469  elementDataStrm >> faceID;
470  Geometry2DSharedPtr face = GetGeometry2D(faceID);
471  if (face == Geometry2DSharedPtr() ||
472  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
473  {
474  std::stringstream errorstring;
475  errorstring << "Element " << indx << " has invalid face: " << faceID;
476  ASSERTL0(false, errorstring.str().c_str());
477  }
478  else if (face->GetShapeType() == LibUtilities::eTriangle)
479  {
480  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
481  faces[Nfaces++] = boost::static_pointer_cast<TriGeom>(face);
482  Ntfaces++;
483  }
484  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
485  {
486  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
487  faces[Nfaces++] = boost::static_pointer_cast<QuadGeom>(face);
488  Nqfaces++;
489  }
490  }
491 
492  /// Make sure all of the face indicies could be read, and that there weren't too few.
493  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for PYRAMID: ") + elementStr).c_str());
494  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
495  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
496 
498  pyrgeom->SetGlobalID(indx);
499 
500  m_pyrGeoms[indx] = pyrgeom;
501  PopulateFaceToElMap(pyrgeom, kNfaces);
502  }
503  catch(...)
504  {
506  (std::string("Unable to read element data for PYRAMID: ") + elementStr).c_str());
507  }
508  }
509  // Prism
510  else if (elementType == "R")
511  {
512  try
513  {
514  /// Create arrays for the tri and quad faces.
515  const int kNfaces = PrismGeom::kNfaces;
516  const int kNtfaces = PrismGeom::kNtfaces;
517  const int kNqfaces = PrismGeom::kNqfaces;
518  Geometry2DSharedPtr faces[kNfaces];
519  int Ntfaces = 0;
520  int Nqfaces = 0;
521  int Nfaces = 0;
522 
523  /// Fill the arrays and make sure there aren't too many faces.
524  std::stringstream errorstring;
525  errorstring << "Element " << indx << " must have "
526  << kNtfaces << " triangle face(s), and "
527  << kNqfaces << " quadrilateral face(s).";
528 
529  for (int i = 0; i < kNfaces; i++)
530  {
531  int faceID;
532  elementDataStrm >> faceID;
533  Geometry2DSharedPtr face = GetGeometry2D(faceID);
534  if (face == Geometry2DSharedPtr() ||
535  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
536  {
537  std::stringstream errorstring;
538  errorstring << "Element " << indx << " has invalid face: " << faceID;
539  ASSERTL0(false, errorstring.str().c_str());
540  }
541  else if (face->GetShapeType() == LibUtilities::eTriangle)
542  {
543  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
544  faces[Nfaces++] = boost::static_pointer_cast<TriGeom>(face);
545  Ntfaces++;
546  }
547  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
548  {
549  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
550  faces[Nfaces++] = boost::static_pointer_cast<QuadGeom>(face);
551  Nqfaces++;
552  }
553  }
554 
555  /// Make sure all of the face indicies could be read, and that there weren't too few.
556  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for PRISM: ") + elementStr).c_str());
557  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
558  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
559 
561  prismgeom->SetGlobalID(indx);
562 
563  m_prismGeoms[indx] = prismgeom;
564  PopulateFaceToElMap(prismgeom, kNfaces);
565  }
566  catch(...)
567  {
569  (std::string("Unable to read element data for PRISM: ") + elementStr).c_str());
570  }
571  }
572  // Hexahedral
573  else if (elementType == "H")
574  {
575  try
576  {
577  /// Create arrays for the tri and quad faces.
578  const int kNfaces = HexGeom::kNfaces;
579  const int kNtfaces = HexGeom::kNtfaces;
580  const int kNqfaces = HexGeom::kNqfaces;
581  //TriGeomSharedPtr tfaces[kNtfaces];
582  QuadGeomSharedPtr qfaces[kNqfaces];
583  int Ntfaces = 0;
584  int Nqfaces = 0;
585 
586  /// Fill the arrays and make sure there aren't too many faces.
587  std::stringstream errorstring;
588  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
589  for (int i = 0; i < kNfaces; i++)
590  {
591  int faceID;
592  elementDataStrm >> faceID;
593  Geometry2DSharedPtr face = GetGeometry2D(faceID);
594  if (face == Geometry2DSharedPtr() ||
595  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
596  {
597  std::stringstream errorstring;
598  errorstring << "Element " << indx << " has invalid face: " << faceID;
599  ASSERTL0(false, errorstring.str().c_str());
600  }
601  else if (face->GetShapeType() == LibUtilities::eTriangle)
602  {
603  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
604  //tfaces[Ntfaces++] = boost::static_pointer_cast<TriGeom>(face);
605  }
606  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
607  {
608  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
609  qfaces[Nqfaces++] = boost::static_pointer_cast<QuadGeom>(face);
610  }
611  }
612 
613  /// Make sure all of the face indicies could be read, and that there weren't too few.
614  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for HEXAHEDRAL: ") + elementStr).c_str());
615  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
616  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
617 
619  hexgeom->SetGlobalID(indx);
620 
621  m_hexGeoms[indx] = hexgeom;
622  PopulateFaceToElMap(hexgeom, kNfaces);
623  }
624  catch(...)
625  {
627  (std::string("Unable to read element data for HEXAHEDRAL: ") + elementStr).c_str());
628  }
629  }
630 
631  /// Keep looking
632  element = element->NextSiblingElement();
633  }
634  }
635 
636  void MeshGraph3D::ReadComposites(TiXmlDocument &doc)
637  {
638  TiXmlHandle docHandle(&doc);
639 
640  /// We know we have it since we made it this far.
641  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
642  TiXmlElement* field = NULL;
643 
644  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
645 
646  /// Look for elements in ELEMENT block.
647  field = mesh->FirstChildElement("COMPOSITE");
648 
649  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
650 
651  int nextCompositeNumber = -1;
652 
653  /// All elements are of the form: "<C ID = "N"> ... </C>".
654 
655  /// Read the ID field first.
656  TiXmlElement *composite = field->FirstChildElement("C");
657 
658  while (composite)
659  {
660  nextCompositeNumber++;
661 
662  int indx;
663  int err = composite->QueryIntAttribute("ID", &indx);
664  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
665 // ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with zero and be sequential.");
666 
667  TiXmlNode* compositeChild = composite->FirstChild();
668  // This is primarily to skip comments that may be present.
669  // Comments appear as nodes just like elements.
670  // We are specifically looking for text in the body
671  // of the definition.
672  while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
673  {
674  compositeChild = compositeChild->NextSibling();
675  }
676 
677  ASSERTL0(compositeChild, "Unable to read composite definition body.");
678  std::string compositeStr = compositeChild->ToText()->ValueStr();
679 
680  /// Parse out the element components corresponding to type of element.
681 
682  std::istringstream compositeDataStrm(compositeStr.c_str());
683 
684  try
685  {
686  bool first = true;
687  std::string prevCompositeElementStr;
688 
689  while (!compositeDataStrm.fail())
690  {
691  std::string compositeElementStr;
692  compositeDataStrm >> compositeElementStr;
693 
694  if (!compositeDataStrm.fail())
695  {
696  if (first)
697  {
698  first = false;
699 
701  m_meshComposites[indx] = curVector;
702  }
703 
704  if (compositeElementStr.length() > 0)
705  {
706  ResolveGeomRef(prevCompositeElementStr, compositeElementStr, m_meshComposites[indx]);
707  }
708  prevCompositeElementStr = compositeElementStr;
709  }
710  }
711  }
712  catch(...)
713  {
715  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
716  }
717 
718  /// Keep looking
719  composite = composite->NextSiblingElement("C");
720  }
721  }
722 
723 
725  {
726  SegGeomSharedPtr returnval;
727  SegGeomMap::iterator x = m_segGeoms.find(eID);
728  ASSERTL0(x != m_segGeoms.end(), "Segment "
729  + boost::lexical_cast<string>(eID) + " not found.");
730  return x->second;
731  };
732 
734  {
735  TriGeomMapIter it1;
736  QuadGeomMapIter it2;
737 
738  it1 = m_triGeoms.find(gID);
739  if (it1 != m_triGeoms.end())
740  return it1->second;
741 
742  it2 = m_quadGeoms.find(gID);
743  if (it2 != m_quadGeoms.end())
744  return it2->second;
745 
746  return Geometry2DSharedPtr();
747  };
748 
749  // Take the string that is the composite reference and find the
750  // pointer to the Geometry object corresponding to it.
751 
752  // The only allowable combinations of previous and current items
753  // are V (0D); E (1D); and T and Q (2D); A (Tet, 3D), P (Pyramid, 3D), R (Prism, 3D), H (Hex, 3D).
754  // Only elements of the same dimension are allowed to be grouped.
755  void MeshGraph3D::ResolveGeomRef(const std::string &prevToken, const std::string &token,
756  Composite& composite)
757  {
758  try
759  {
760  std::istringstream tokenStream(token);
761  std::istringstream prevTokenStream(prevToken);
762 
763  char type;
764  char prevType;
765 
766  tokenStream >> type;
767 
768  std::string::size_type indxBeg = token.find_first_of('[') + 1;
769  std::string::size_type indxEnd = token.find_last_of(']') - 1;
770 
771  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());
772 
773  std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
774 
775  std::vector<unsigned int> seqVector;
777 
778  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
779 
780  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
781 
782  prevTokenStream >> prevType;
783 
784  // All composites must be of the same dimension. This map makes things clean to compare.
785  map<char, int> typeMap;
786  typeMap['V'] = 1; // Vertex
787  typeMap['E'] = 1; // Edge
788  typeMap['T'] = 2; // Triangle
789  typeMap['Q'] = 2; // Quad
790  typeMap['A'] = 3; // Tet
791  typeMap['P'] = 3; // Pyramid
792  typeMap['R'] = 3; // Prism
793  typeMap['H'] = 3; // Hex
794 
795  // Make sure only geoms of the same dimension are combined.
796  bool validSequence = (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
797 
798  ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
799  + type + " and " + prevType + ".").c_str());
800 
801  switch(type)
802  {
803  case 'V': // Vertex
804  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
805  {
806  if (m_vertSet.find(*seqIter) == m_vertSet.end())
807  {
808  char errStr[16] = "";
809  ::sprintf(errStr, "%d", *seqIter);
810  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
811  }
812  else
813  {
814  composite->push_back(m_vertSet[*seqIter]);
815  }
816  }
817  break;
818 
819  case 'E': // Edge
820  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
821  {
822  if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
823  {
824  char errStr[16] = "";
825  ::sprintf(errStr, "%d", *seqIter);
826  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
827  }
828  else
829  {
830  composite->push_back(m_segGeoms[*seqIter]);
831  }
832  }
833  break;
834 
835  case 'F': // Face
836  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
837  {
838  Geometry2DSharedPtr face = GetGeometry2D(*seqIter);
839  if (face == Geometry2DSharedPtr())
840  {
841  char errStr[16] = "";
842  ::sprintf(errStr, "%d", *seqIter);
843  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown face index: ") + errStr).c_str());
844  }
845  else
846  {
847  if(CheckRange(*face))
848  {
849  composite->push_back(face);
850  }
851  }
852  }
853  break;
854 
855  case 'T': // Triangle
856  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
857  {
858  if (m_triGeoms.find(*seqIter) == m_triGeoms.end())
859  {
860  char errStr[16] = "";
861  ::sprintf(errStr, "%d", *seqIter);
862  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
863  }
864  else
865  {
866  if(CheckRange(*m_triGeoms[*seqIter]))
867  {
868  composite->push_back(m_triGeoms[*seqIter]);
869  }
870  }
871  }
872  break;
873 
874  case 'Q': // Quad
875  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
876  {
877  if (m_quadGeoms.find(*seqIter) == m_quadGeoms.end())
878  {
879  char errStr[16] = "";
880  ::sprintf(errStr, "%d", *seqIter);
881  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr).c_str());
882  }
883  else
884  {
885  if(CheckRange(*m_quadGeoms[*seqIter]))
886  {
887  composite->push_back(m_quadGeoms[*seqIter]);
888  }
889  }
890  }
891  break;
892 
893  // Tetrahedron
894  case 'A':
895  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
896  {
897  if (m_tetGeoms.find(*seqIter) == m_tetGeoms.end())
898  {
899  char errStr[16] = "";
900  ::sprintf(errStr, "%d", *seqIter);
901  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown tet index: ") + errStr).c_str());
902  }
903  else
904  {
905  if(CheckRange(*m_tetGeoms[*seqIter]))
906  {
907  composite->push_back(m_tetGeoms[*seqIter]);
908  }
909  }
910  }
911  break;
912 
913  // Pyramid
914  case 'P':
915  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
916  {
917  if (m_pyrGeoms.find(*seqIter) == m_pyrGeoms.end())
918  {
919  char errStr[16] = "";
920  ::sprintf(errStr, "%d", *seqIter);
921  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown pyramid index: ") + errStr).c_str());
922  }
923  else
924  {
925  if(CheckRange(*m_pyrGeoms[*seqIter]))
926  {
927  composite->push_back(m_pyrGeoms[*seqIter]);
928  }
929  }
930  }
931  break;
932 
933  // Prism
934  case 'R':
935  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
936  {
937  if (m_prismGeoms.find(*seqIter) == m_prismGeoms.end())
938  {
939  char errStr[16] = "";
940  ::sprintf(errStr, "%d", *seqIter);
941  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown prism index: ") + errStr).c_str());
942  }
943  else
944  {
945  if(CheckRange(*m_prismGeoms[*seqIter]))
946  {
947  composite->push_back(m_prismGeoms[*seqIter]);
948  }
949  }
950  }
951  break;
952 
953  // Hex
954  case 'H':
955  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
956  {
957  if (m_hexGeoms.find(*seqIter) == m_hexGeoms.end())
958  {
959  char errStr[16] = "";
960  ::sprintf(errStr, "%d", *seqIter);
961  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown hex index: ") + errStr).c_str());
962  }
963  else
964  {
965  if(CheckRange(*m_hexGeoms[*seqIter]))
966  {
967  composite->push_back(m_hexGeoms[*seqIter]);
968  }
969  }
970  }
971  break;
972 
973  default:
974  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
975  }
976  }
977  catch(...)
978  {
979  NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
980  }
981 
982  return;
983  }
984 
985 
987  {
989  m_faceToElMap.find(face->GetGlobalID());
990 
991  ASSERTL0(it != m_faceToElMap.end(), "Unable to find corresponding face!");
992 
993  return it->second;
994  }
995 
996 
997  /**
998  * Retrieve the basis key for a given face direction.
999  */
1001  Geometry2DSharedPtr face,
1002  const int facedir,
1003  const std::string variable)
1004  {
1005  // Retrieve the list of elements and the associated face index
1006  // to which the face geometry belongs.
1008 
1009  ASSERTL0(elements->size() > 0, "No elements for the given face."
1010  " Check all elements belong to the domain composite.");
1011 
1012  // Perhaps, a check should be done here to ensure that in case
1013  // elements->size!=1, all elements to which the edge belongs have
1014  // the same type and order of expansion such that no confusion can
1015  // arise.
1016 
1017  // Get the Expansion structure detailing the basis keys used for
1018  // this element.
1019  ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element,
1020  variable);
1021 
1022  // Retrieve the geometry object of the element as a Geometry3D.
1023  Geometry3DSharedPtr geom3d =
1024  boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
1025  expansion->m_geomShPtr);
1026 
1027  // Use the geometry of the element to calculate the coordinate
1028  // direction of the element which corresponds to the requested
1029  // coordinate direction of the given face.
1030  int dir = geom3d->GetDir((*elements)[0]->m_FaceIndx, facedir);
1031 
1032  if(face->GetNumVerts() == 3)
1033  {
1034  return StdRegions::EvaluateTriFaceBasisKey(facedir,
1035  expansion->m_basisKeyVector[dir].GetBasisType(),
1036  expansion->m_basisKeyVector[dir].GetNumPoints(),
1037  expansion->m_basisKeyVector[dir].GetNumModes());
1038  }
1039  else
1040  {
1041  return StdRegions::EvaluateQuadFaceBasisKey(facedir,
1042  expansion->m_basisKeyVector[dir].GetBasisType(),
1043  expansion->m_basisKeyVector[dir].GetNumPoints(),
1044  expansion->m_basisKeyVector[dir].GetNumModes());
1045  }
1046 
1047  // Keep things happy by returning a value.
1049  }
1050 
1051 
1052  /**
1053  * @brief Given a 3D geometry object #element, populate the face to
1054  * element map #m_faceToElMap which maps faces to their corresponding
1055  * element(s).
1056  *
1057  * @param element Element to process.
1058  * @param kNfaces Number of faces of #element. Should be removed and
1059  * put into Geometry3D as a virtual member function.
1060  */
1062  {
1063  // Set up face -> element map
1064  for (int i = 0; i < kNfaces; ++i)
1065  {
1066  int faceId = element->GetFace(i)->GetGlobalID();
1067  ElementFaceSharedPtr elementFace =
1069 
1070  elementFace->m_Element = element;
1071  elementFace->m_FaceIndx = i;
1072 
1073  // Search map to see if face already exists.
1075  m_faceToElMap.find(faceId);
1076 
1077  if (it == m_faceToElMap.end())
1078  {
1081  tmp->push_back(elementFace);
1082  m_faceToElMap[faceId] = tmp;
1083  }
1084  else
1085  {
1086  ElementFaceVectorSharedPtr tmp = it->second;
1087  tmp->push_back(elementFace);
1088  }
1089  }
1090  }
1091  }; //end of namespace
1092 }; //end of namespace