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 
37 #include <SpatialDomains/MeshGraph3D.h>
38 #include <SpatialDomains/TriGeom.h>
39 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
40 #include <tinyxml.h>
41 
42 namespace Nektar
43 {
44  namespace SpatialDomains
45  {
46  MeshGraph3D::MeshGraph3D() : MeshGraph(3,3)
47  {
48  }
49 
50  MeshGraph3D::MeshGraph3D(const LibUtilities::SessionReaderSharedPtr &pSession,
51  const DomainRangeShPtr &rng)
52  : MeshGraph(pSession,rng)
53  {
54  ReadGeometry(pSession->GetDocument());
55  ReadExpansions(pSession->GetDocument());
56  }
57 
58  MeshGraph3D::~MeshGraph3D()
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
80  MeshGraph::ReadGeometry(doc);
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  CurveMap::iterator it;
106 
107  /// Look for elements in ELEMENT block.
108  field = mesh->FirstChildElement("EDGE");
109 
110  ASSERTL0(field, "Unable to find EDGE tag in file.");
111 
112  string IsCompressed;
113  field->QueryStringAttribute("COMPRESSED",&IsCompressed);
114 
115  if(IsCompressed.size())
116  {
117  ASSERTL0(boost::iequals(IsCompressed,
118  LibUtilities::CompressData::GetCompressString()),
119  "Compressed formats do not match. Expected :"
120  + LibUtilities::CompressData::GetCompressString()
121  + " but got " + std::string(IsCompressed));
122  // Extract the edge body
123  TiXmlNode* edgeChild = field->FirstChild();
124  ASSERTL0(edgeChild, "Unable to extract the data from "
125  "the compressed edge tag.");
126 
127  std::string edgeStr;
128  if (edgeChild->Type() == TiXmlNode::TINYXML_TEXT)
129  {
130  edgeStr += edgeChild->ToText()->ValueStr();
131  }
132 
133  std::vector<LibUtilities::MeshEdge> edgeData;
134  LibUtilities::CompressData::ZlibDecodeFromBase64Str(edgeStr,
135  edgeData);
136 
137  int indx;
138  for(int i = 0; i < edgeData.size(); ++i)
139  {
140  indx = edgeData[i].id;
141  PointGeomSharedPtr vertices[2] = {
142  GetVertex(edgeData[i].v0),
143  GetVertex(edgeData[i].v1)
144  };
145  SegGeomSharedPtr edge;
146 
147  it = m_curvedEdges.find(indx);
148  if (it == m_curvedEdges.end())
149  {
150  edge = MemoryManager<SegGeom>::AllocateSharedPtr(
151  indx, m_spaceDimension, vertices);
152  }
153  else
154  {
155  edge = MemoryManager<SegGeom>::AllocateSharedPtr(
156  indx, m_spaceDimension, vertices,
157  it->second);
158  }
159  m_segGeoms[indx] = edge;
160  }
161  }
162  else
163  {
164 
165  /// All elements are of the form: "<E ID="#"> ... </E>", with
166  /// ? being the element type.
167  /// Read the ID field first.
168  TiXmlElement *edge = field->FirstChildElement("E");
169 
170  /// Since all edge data is one big text block, we need to accumulate
171  /// all TINYXML_TEXT data and then parse it. This approach effectively skips
172  /// all comments or other node types since we only care about the
173  /// edge list. We cannot handle missing edge numbers as we could
174  /// with missing element numbers due to the text block format.
175  std::string edgeStr;
176  int indx;
177 
178  while(edge)
179  {
180  int err = edge->QueryIntAttribute("ID",&indx);
181  ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
182 
183  TiXmlNode *child = edge->FirstChild();
184  edgeStr.clear();
185  if (child->Type() == TiXmlNode::TINYXML_TEXT)
186  {
187  edgeStr += child->ToText()->ValueStr();
188  }
189 
190  /// Now parse out the edges, three fields at a time.
191  int vertex1, vertex2;
192  std::istringstream edgeDataStrm(edgeStr.c_str());
193 
194  try
195  {
196  while (!edgeDataStrm.fail())
197  {
198  edgeDataStrm >> vertex1 >> vertex2;
199 
200  // Must check after the read because we
201  // may be at the end and not know it. If
202  // we are at the end we will add a
203  // duplicate of the last entry if we don't
204  // check here.
205  if (!edgeDataStrm.fail())
206  {
207  PointGeomSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
208  SegGeomSharedPtr edge;
209  it = m_curvedEdges.find(indx);
210 
211  if (it == m_curvedEdges.end())
212  {
213  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices);
214  }
215  else
216  {
217  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices, it->second);
218  }
219 
220  m_segGeoms[indx] = edge;
221  }
222  }
223  }
224  catch(...)
225  {
226  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
227  }
228 
229  edge = edge->NextSiblingElement("E");
230  }
231  }
232  }
233 
234  void MeshGraph3D::ReadFaces(TiXmlDocument &doc)
235  {
236  /// We know we have it since we made it this far.
237  TiXmlHandle docHandle(&doc);
238  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
239  TiXmlElement* field = NULL;
240 
241  /// Look for elements in FACE block.
242  field = mesh->FirstChildElement("FACE");
243 
244  ASSERTL0(field, "Unable to find FACE tag in file.");
245 
246  /// All faces are of the form: "<? ID="#"> ... </?>", with
247  /// ? being an element type (either Q or T).
248  /// They might be in compressed format and so then need upacking.
249 
250  TiXmlElement *element = field->FirstChildElement();
251  CurveMap::iterator it;
252 
253  while (element)
254  {
255  std::string elementType(element->ValueStr());
256 
257  ASSERTL0(elementType == "Q" || elementType == "T",
258  (std::string("Unknown 3D face type: ") + elementType).c_str());
259 
260 
261  string IsCompressed;
262  element->QueryStringAttribute("COMPRESSED",&IsCompressed);
263 
264  if(IsCompressed.size())
265  {
266  ASSERTL0(boost::iequals(IsCompressed,
267  LibUtilities::CompressData:: GetCompressString()),
268  "Compressed formats do not match. Expected :"
269  + LibUtilities::CompressData::GetCompressString()
270  + " but got "+ std::string(IsCompressed));
271 
272  // Extract the face body
273  TiXmlNode* faceChild = element->FirstChild();
274  ASSERTL0(faceChild, "Unable to extract the data from "
275  "the compressed face tag.");
276 
277  std::string faceStr;
278  if (faceChild->Type() == TiXmlNode::TINYXML_TEXT)
279  {
280  faceStr += faceChild->ToText()->ValueStr();
281  }
282 
283  int indx;
284  if(elementType == "T")
285  {
286  std::vector<LibUtilities::MeshTri> faceData;
287  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
288  faceStr,faceData);
289 
290  for(int i = 0; i < faceData.size(); ++i)
291  {
292  indx = faceData[i].id;
293 
294  /// See if this face has curves.
295  it = m_curvedFaces.find(indx);
296 
297  /// Create a TriGeom to hold the new definition.
298  SegGeomSharedPtr edges[TriGeom::kNedges] =
299  {
300  GetSegGeom(faceData[i].e[0]),
301  GetSegGeom(faceData[i].e[1]),
302  GetSegGeom(faceData[i].e[2])
303  };
304 
305  StdRegions::Orientation edgeorient[TriGeom::kNedges] =
306  {
307  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
308  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
309  SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
310  };
311 
312  TriGeomSharedPtr trigeom;
313  if (it == m_curvedFaces.end())
314  {
315  trigeom =
316  MemoryManager<TriGeom>::AllocateSharedPtr(
317  indx, edges, edgeorient);
318  }
319  else
320  {
321  trigeom =
322  MemoryManager<TriGeom>::AllocateSharedPtr(
323  indx, edges, edgeorient,
324  it->second);
325  }
326  trigeom->SetGlobalID(indx);
327  m_triGeoms[indx] = trigeom;
328  }
329  }
330  else if (elementType == "Q")
331  {
332  std::vector<LibUtilities::MeshQuad> faceData;
333  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
334  faceStr,faceData);
335 
336  for(int i = 0; i < faceData.size(); ++i)
337  {
338  indx = faceData[i].id;
339 
340  /// See if this face has curves.
341  it = m_curvedFaces.find(indx);
342 
343  /// Create a QuadGeom to hold the new definition.
344  SegGeomSharedPtr edges[QuadGeom::kNedges] = {
345  GetSegGeom(faceData[i].e[0]),
346  GetSegGeom(faceData[i].e[1]),
347  GetSegGeom(faceData[i].e[2]),
348  GetSegGeom(faceData[i].e[3])
349  };
350 
351  StdRegions::Orientation edgeorient[QuadGeom::kNedges] =
352  {
353  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
354  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
355  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
356  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
357  };
358 
359  QuadGeomSharedPtr quadgeom;
360  if (it == m_curvedEdges.end())
361  {
362  quadgeom =
363  MemoryManager<QuadGeom>::AllocateSharedPtr(
364  indx, edges, edgeorient);
365  }
366  else
367  {
368  quadgeom =
369  MemoryManager<QuadGeom>::AllocateSharedPtr(
370  indx, edges, edgeorient,
371  it->second);
372  }
373  quadgeom->SetGlobalID(indx);
374  m_quadGeoms[indx] = quadgeom;
375  }
376  }
377  }
378  else
379  {
380  /// Read id attribute.
381  int indx;
382  int err = element->QueryIntAttribute("ID", &indx);
383  ASSERTL0(err == TIXML_SUCCESS, "Unable to read face attribute ID.");
384 
385  /// See if this face has curves.
386  it = m_curvedFaces.find(indx);
387 
388  /// Read text element description.
389  TiXmlNode* elementChild = element->FirstChild();
390  std::string elementStr;
391  while(elementChild)
392  {
393  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
394  {
395  elementStr += elementChild->ToText()->ValueStr();
396  }
397  elementChild = elementChild->NextSibling();
398  }
399 
400  ASSERTL0(!elementStr.empty(), "Unable to read face description body.");
401 
402  /// Parse out the element components corresponding to type of element.
403  if (elementType == "T")
404  {
405  // Read three edge numbers
406  int edge1, edge2, edge3;
407  std::istringstream elementDataStrm(elementStr.c_str());
408 
409  try
410  {
411  elementDataStrm >> edge1;
412  elementDataStrm >> edge2;
413  elementDataStrm >> edge3;
414 
415  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read face data for TRIANGLE: ") + elementStr).c_str());
416 
417  /// Create a TriGeom to hold the new definition.
418  SegGeomSharedPtr edges[TriGeom::kNedges] =
419  {
420  GetSegGeom(edge1),
421  GetSegGeom(edge2),
422  GetSegGeom(edge3)
423  };
424 
425  StdRegions::Orientation edgeorient[TriGeom::kNedges] =
426  {
427  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
428  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
429  SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
430  };
431 
432  TriGeomSharedPtr trigeom;
433 
434  if (it == m_curvedFaces.end())
435  {
436  trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, edgeorient);
437  }
438  else
439  {
440  trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, edgeorient, it->second);
441  }
442 
443  trigeom->SetGlobalID(indx);
444 
445  m_triGeoms[indx] = trigeom;
446  }
447  catch(...)
448  {
449  NEKERROR(ErrorUtil::efatal,
450  (std::string("Unable to read face data for TRIANGLE: ") + elementStr).c_str());
451  }
452  }
453  else if (elementType == "Q")
454  {
455  // Read four edge numbers
456  int edge1, edge2, edge3, edge4;
457  std::istringstream elementDataStrm(elementStr.c_str());
458 
459  try
460  {
461  elementDataStrm >> edge1;
462  elementDataStrm >> edge2;
463  elementDataStrm >> edge3;
464  elementDataStrm >> edge4;
465 
466  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read face data for QUAD: ") + elementStr).c_str());
467 
468  /// Create a QuadGeom to hold the new definition.
469  SegGeomSharedPtr edges[QuadGeom::kNedges] =
470  {GetSegGeom(edge1),GetSegGeom(edge2),
471  GetSegGeom(edge3),GetSegGeom(edge4)};
472 
473  StdRegions::Orientation edgeorient[QuadGeom::kNedges] =
474  {
475  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
476  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
477  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
478  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
479  };
480 
481  QuadGeomSharedPtr quadgeom;
482 
483  if (it == m_curvedEdges.end())
484  {
485  quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient);
486  }
487  else
488  {
489  quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient, it->second);
490  }
491  quadgeom->SetGlobalID(indx);
492 
493  m_quadGeoms[indx] = quadgeom;
494 
495  }
496  catch(...)
497  {
498  NEKERROR(ErrorUtil::efatal,(std::string("Unable to read face data for QUAD: ") + elementStr).c_str());
499  }
500  }
501  }
502  /// Keep looking
503  element = element->NextSiblingElement();
504  }
505  }
506 
507  void MeshGraph3D::ReadElements(TiXmlDocument &doc)
508  {
509  /// We know we have it since we made it this far.
510  TiXmlHandle docHandle(&doc);
511  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
512  TiXmlElement* field = NULL;
513 
514  /// Look for elements in ELEMENT block.
515  field = mesh->FirstChildElement("ELEMENT");
516 
517  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
518 
519  /// All elements are of the form: "<? ID="#"> ... </?>", with
520  /// ? being the element type.
521 
522  TiXmlElement *element = field->FirstChildElement();
523 
524  while (element)
525  {
526  std::string elementType(element->ValueStr());
527 
528  //A - tet, P - pyramid, R - prism, H - hex
529  ASSERTL0(elementType == "A" || elementType == "P" || elementType == "R" || elementType == "H",
530  (std::string("Unknown 3D element type: ") + elementType).c_str());
531 
532 
533  string IsCompressed;
534  element->QueryStringAttribute("COMPRESSED",&IsCompressed);
535 
536  if(IsCompressed.size())
537  {
538  ASSERTL0(boost::iequals(IsCompressed,
539  LibUtilities::CompressData:: GetCompressString()),
540  "Compressed formats do not match. Expected :"
541  + LibUtilities::CompressData::GetCompressString()
542  + " but got " + std::string(IsCompressed));
543 
544  // Extract the face body
545  TiXmlNode* child = element->FirstChild();
546  ASSERTL0(child, "Unable to extract the data from "
547  "the compressed face tag.");
548 
549  std::string str;
550  if (child->Type() == TiXmlNode::TINYXML_TEXT)
551  {
552  str += child->ToText()->ValueStr();
553  }
554 
555  int indx;
556  if(elementType == "A")
557  {
558  std::vector<LibUtilities::MeshTet> data;
559  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
560  str,data);
561  TriGeomSharedPtr tfaces[4];
562  for(int i = 0; i < data.size(); ++i)
563  {
564  indx = data[i].id;
565  for(int j = 0; j < 4; ++j)
566  {
567  Geometry2DSharedPtr face =
568  GetGeometry2D(data[i].f[j]);
569  tfaces[j] =
570  boost::static_pointer_cast<TriGeom>(face);
571  }
572 
573  TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>
574  ::AllocateSharedPtr(tfaces));
575  tetgeom->SetGlobalID(indx);
576  m_tetGeoms[indx] = tetgeom;
577  PopulateFaceToElMap(tetgeom, 4);
578  }
579  }
580  else if (elementType == "P")
581  {
582  std::vector<LibUtilities::MeshPyr> data;
583  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
584  str,data);
585  Geometry2DSharedPtr faces[5];
586  for(int i = 0; i < data.size(); ++i)
587  {
588  indx = data[i].id;
589  int Ntfaces = 0;
590  int Nqfaces = 0;
591  for(int j = 0; j < 5; ++j)
592  {
593  Geometry2DSharedPtr face =
594  GetGeometry2D(data[i].f[j]);
595 
596  if (face == Geometry2DSharedPtr() ||
597  (face->GetShapeType() !=
598  LibUtilities::eTriangle &&
599  face->GetShapeType() !=
600  LibUtilities::eQuadrilateral))
601  {
602  std::stringstream errorstring;
603  errorstring << "Element " << indx
604  << " has invalid face: " << j;
605  ASSERTL0(false, errorstring.str().c_str());
606  }
607  else if (face->GetShapeType() ==
608  LibUtilities::eTriangle)
609  {
610  faces[j] = boost
611  ::static_pointer_cast<TriGeom>(face);
612  Ntfaces++;
613  }
614  else if (face->GetShapeType() ==
615  LibUtilities::eQuadrilateral)
616  {
617  faces[j] = boost
618  ::static_pointer_cast<QuadGeom>(face);
619  Nqfaces++;
620  }
621  }
622  ASSERTL0((Ntfaces == 4) && (Nqfaces = 1),
623  "Did not identify the correct number of "
624  "triangular and quadrilateral faces for a "
625  "pyramid");
626 
627  PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>
628  ::AllocateSharedPtr(faces));
629  pyrgeom->SetGlobalID(indx);
630  m_pyrGeoms[indx] = pyrgeom;
631  PopulateFaceToElMap(pyrgeom, 5);
632  }
633  }
634  else if (elementType == "R")
635  {
636  std::vector<LibUtilities::MeshPrism> data;
637  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
638  str,data);
639  Geometry2DSharedPtr faces[5];
640  for(int i = 0; i < data.size(); ++i)
641  {
642  indx = data[i].id;
643  int Ntfaces = 0;
644  int Nqfaces = 0;
645  for(int j = 0; j < 5; ++j)
646  {
647  Geometry2DSharedPtr face =
648  GetGeometry2D(data[i].f[j]);
649  if (face == Geometry2DSharedPtr() ||
650  (face->GetShapeType() !=
651  LibUtilities::eTriangle &&
652  face->GetShapeType() !=
653  LibUtilities::eQuadrilateral))
654  {
655  std::stringstream errorstring;
656  errorstring << "Element " << indx
657  << " has invalid face: " << j;
658  ASSERTL0(false, errorstring.str().c_str());
659  }
660  else if (face->GetShapeType() ==
661  LibUtilities::eTriangle)
662  {
663  faces[j] = boost
664  ::static_pointer_cast<TriGeom>(face);
665  Ntfaces++;
666  }
667  else if (face->GetShapeType() ==
668  LibUtilities::eQuadrilateral)
669  {
670  faces[j] = boost
671  ::static_pointer_cast<QuadGeom>(face);
672  Nqfaces++;
673  }
674  }
675  ASSERTL0((Ntfaces == 2) && (Nqfaces = 3),
676  "Did not identify the correct number of "
677  "triangular and quadrilateral faces for a "
678  "prism");
679 
680  PrismGeomSharedPtr prismgeom(
681  MemoryManager<PrismGeom>
682  ::AllocateSharedPtr(faces));
683  prismgeom->SetGlobalID(indx);
684  m_prismGeoms[indx] = prismgeom;
685  PopulateFaceToElMap(prismgeom, 5);
686  }
687  }
688  else if (elementType == "H")
689  {
690  std::vector<LibUtilities::MeshHex> data;
691  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
692  str,data);
693 
694  QuadGeomSharedPtr faces[6];
695  for(int i = 0; i < data.size(); ++i)
696  {
697  indx = data[i].id;
698  for(int j = 0; j < 6; ++j)
699  {
700  Geometry2DSharedPtr face =
701  GetGeometry2D(data[i].f[j]);
702  faces[j] = boost
703  ::static_pointer_cast<QuadGeom>(face);
704  }
705 
706  HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>
707  ::AllocateSharedPtr(faces));
708  hexgeom->SetGlobalID(indx);
709  m_hexGeoms[indx] = hexgeom;
710  PopulateFaceToElMap(hexgeom, 6);
711  }
712  }
713  }
714  else
715  {
716  /// Read id attribute.
717  int indx;
718  int err = element->QueryIntAttribute("ID", &indx);
719  ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
720 
721  /// Read text element description.
722  TiXmlNode* elementChild = element->FirstChild();
723  std::string elementStr;
724  while(elementChild)
725  {
726  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
727  {
728  elementStr += elementChild->ToText()->ValueStr();
729  }
730  elementChild = elementChild->NextSibling();
731  }
732 
733  ASSERTL0(!elementStr.empty(), "Unable to read element description body.");
734 
735  std::istringstream elementDataStrm(elementStr.c_str());
736 
737  /// Parse out the element components corresponding to type of element.
738 
739  // Tetrahedral
740  if (elementType == "A")
741  {
742  try
743  {
744  /// Create arrays for the tri and quad faces.
745  const int kNfaces = TetGeom::kNfaces;
746  const int kNtfaces = TetGeom::kNtfaces;
747  const int kNqfaces = TetGeom::kNqfaces;
748  TriGeomSharedPtr tfaces[kNtfaces];
749  int Ntfaces = 0;
750  int Nqfaces = 0;
751 
752  /// Fill the arrays and make sure there aren't too many faces.
753  std::stringstream errorstring;
754  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
755  for (int i = 0; i < kNfaces; i++)
756  {
757  int faceID;
758  elementDataStrm >> faceID;
759  Geometry2DSharedPtr face = GetGeometry2D(faceID);
760  if (face == Geometry2DSharedPtr() ||
761  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
762  {
763  std::stringstream errorstring;
764  errorstring << "Element " << indx << " has invalid face: " << faceID;
765  ASSERTL0(false, errorstring.str().c_str());
766  }
767  else if (face->GetShapeType() == LibUtilities::eTriangle)
768  {
769  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
770  tfaces[Ntfaces++] = boost::static_pointer_cast<TriGeom>(face);
771  }
772  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
773  {
774  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
775  }
776  }
777 
778  /// Make sure all of the face indicies could be read, and that there weren't too few.
779  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TETRAHEDRON: ") + elementStr).c_str());
780  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
781  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
782 
783  TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
784  tetgeom->SetGlobalID(indx);
785 
786  m_tetGeoms[indx] = tetgeom;
787  PopulateFaceToElMap(tetgeom, kNfaces);
788  }
789  catch(...)
790  {
791  NEKERROR(ErrorUtil::efatal,
792  (std::string("Unable to read element data for TETRAHEDRON: ") + elementStr).c_str());
793  }
794  }
795  // Pyramid
796  else if (elementType == "P")
797  {
798  try
799  {
800  /// Create arrays for the tri and quad faces.
801  const int kNfaces = PyrGeom::kNfaces;
802  const int kNtfaces = PyrGeom::kNtfaces;
803  const int kNqfaces = PyrGeom::kNqfaces;
804  Geometry2DSharedPtr faces[kNfaces];
805  int Nfaces = 0;
806  int Ntfaces = 0;
807  int Nqfaces = 0;
808 
809  /// Fill the arrays and make sure there aren't too many faces.
810  std::stringstream errorstring;
811  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
812  for (int i = 0; i < kNfaces; i++)
813  {
814  int faceID;
815  elementDataStrm >> faceID;
816  Geometry2DSharedPtr face = GetGeometry2D(faceID);
817  if (face == Geometry2DSharedPtr() ||
818  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
819  {
820  std::stringstream errorstring;
821  errorstring << "Element " << indx << " has invalid face: " << faceID;
822  ASSERTL0(false, errorstring.str().c_str());
823  }
824  else if (face->GetShapeType() == LibUtilities::eTriangle)
825  {
826  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
827  faces[Nfaces++] = boost::static_pointer_cast<TriGeom>(face);
828  Ntfaces++;
829  }
830  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
831  {
832  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
833  faces[Nfaces++] = boost::static_pointer_cast<QuadGeom>(face);
834  Nqfaces++;
835  }
836  }
837 
838  /// Make sure all of the face indicies could be read, and that there weren't too few.
839  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for PYRAMID: ") + elementStr).c_str());
840  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
841  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
842 
843  PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>::AllocateSharedPtr(faces));
844  pyrgeom->SetGlobalID(indx);
845 
846  m_pyrGeoms[indx] = pyrgeom;
847  PopulateFaceToElMap(pyrgeom, kNfaces);
848  }
849  catch(...)
850  {
851  NEKERROR(ErrorUtil::efatal,
852  (std::string("Unable to read element data for PYRAMID: ") + elementStr).c_str());
853  }
854  }
855  // Prism
856  else if (elementType == "R")
857  {
858  try
859  {
860  /// Create arrays for the tri and quad faces.
861  const int kNfaces = PrismGeom::kNfaces;
862  const int kNtfaces = PrismGeom::kNtfaces;
863  const int kNqfaces = PrismGeom::kNqfaces;
864  Geometry2DSharedPtr faces[kNfaces];
865  int Ntfaces = 0;
866  int Nqfaces = 0;
867  int Nfaces = 0;
868 
869  /// Fill the arrays and make sure there aren't too many faces.
870  std::stringstream errorstring;
871  errorstring << "Element " << indx << " must have "
872  << kNtfaces << " triangle face(s), and "
873  << kNqfaces << " quadrilateral face(s).";
874 
875  for (int i = 0; i < kNfaces; i++)
876  {
877  int faceID;
878  elementDataStrm >> faceID;
879  Geometry2DSharedPtr face = GetGeometry2D(faceID);
880  if (face == Geometry2DSharedPtr() ||
881  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
882  {
883  std::stringstream errorstring;
884  errorstring << "Element " << indx << " has invalid face: " << faceID;
885  ASSERTL0(false, errorstring.str().c_str());
886  }
887  else if (face->GetShapeType() == LibUtilities::eTriangle)
888  {
889  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
890  faces[Nfaces++] = boost::static_pointer_cast<TriGeom>(face);
891  Ntfaces++;
892  }
893  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
894  {
895  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
896  faces[Nfaces++] = boost::static_pointer_cast<QuadGeom>(face);
897  Nqfaces++;
898  }
899  }
900 
901  /// Make sure all of the face indicies could be read, and that there weren't too few.
902  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for PRISM: ") + elementStr).c_str());
903  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
904  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
905 
906  PrismGeomSharedPtr prismgeom(MemoryManager<PrismGeom>::AllocateSharedPtr(faces));
907  prismgeom->SetGlobalID(indx);
908 
909  m_prismGeoms[indx] = prismgeom;
910  PopulateFaceToElMap(prismgeom, kNfaces);
911  }
912  catch(...)
913  {
914  NEKERROR(ErrorUtil::efatal,
915  (std::string("Unable to read element data for PRISM: ") + elementStr).c_str());
916  }
917  }
918  // Hexahedral
919  else if (elementType == "H")
920  {
921  try
922  {
923  /// Create arrays for the tri and quad faces.
924  const int kNfaces = HexGeom::kNfaces;
925  const int kNtfaces = HexGeom::kNtfaces;
926  const int kNqfaces = HexGeom::kNqfaces;
927  //TriGeomSharedPtr tfaces[kNtfaces];
928  QuadGeomSharedPtr qfaces[kNqfaces];
929  int Ntfaces = 0;
930  int Nqfaces = 0;
931 
932  /// Fill the arrays and make sure there aren't too many faces.
933  std::stringstream errorstring;
934  errorstring << "Element " << indx << " must have " << kNtfaces << " triangle face(s), and " << kNqfaces << " quadrilateral face(s).";
935  for (int i = 0; i < kNfaces; i++)
936  {
937  int faceID;
938  elementDataStrm >> faceID;
939  Geometry2DSharedPtr face = GetGeometry2D(faceID);
940  if (face == Geometry2DSharedPtr() ||
941  (face->GetShapeType() != LibUtilities::eTriangle && face->GetShapeType() != LibUtilities::eQuadrilateral))
942  {
943  std::stringstream errorstring;
944  errorstring << "Element " << indx << " has invalid face: " << faceID;
945  ASSERTL0(false, errorstring.str().c_str());
946  }
947  else if (face->GetShapeType() == LibUtilities::eTriangle)
948  {
949  ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
950  //tfaces[Ntfaces++] = boost::static_pointer_cast<TriGeom>(face);
951  }
952  else if (face->GetShapeType() == LibUtilities::eQuadrilateral)
953  {
954  ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
955  qfaces[Nqfaces++] = boost::static_pointer_cast<QuadGeom>(face);
956  }
957  }
958 
959  /// Make sure all of the face indicies could be read, and that there weren't too few.
960  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for HEXAHEDRAL: ") + elementStr).c_str());
961  ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
962  ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
963 
964  HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>::AllocateSharedPtr(qfaces));
965  hexgeom->SetGlobalID(indx);
966 
967  m_hexGeoms[indx] = hexgeom;
968  PopulateFaceToElMap(hexgeom, kNfaces);
969  }
970  catch(...)
971  {
972  NEKERROR(ErrorUtil::efatal,
973  (std::string("Unable to read element data for HEXAHEDRAL: ") + elementStr).c_str());
974  }
975  }
976  }
977  /// Keep looking
978  element = element->NextSiblingElement();
979  }
980  }
981 
982  void MeshGraph3D::ReadComposites(TiXmlDocument &doc)
983  {
984  TiXmlHandle docHandle(&doc);
985 
986  /// We know we have it since we made it this far.
987  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
988  TiXmlElement* field = NULL;
989 
990  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
991 
992  /// Look for elements in ELEMENT block.
993  field = mesh->FirstChildElement("COMPOSITE");
994 
995  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
996 
997  int nextCompositeNumber = -1;
998 
999  /// All elements are of the form: "<C ID = "N"> ... </C>".
1000 
1001  /// Read the ID field first.
1002  TiXmlElement *composite = field->FirstChildElement("C");
1003 
1004  while (composite)
1005  {
1006  nextCompositeNumber++;
1007 
1008  int indx;
1009  int err = composite->QueryIntAttribute("ID", &indx);
1010  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
1011 
1012  // read and store label if they exist
1013  string labelstr;
1014  err = composite->QueryStringAttribute("LABEL", &labelstr);
1015  if(err == TIXML_SUCCESS)
1016  {
1017  m_compositesLabels[indx] = labelstr;
1018  }
1019 
1020  TiXmlNode* compositeChild = composite->FirstChild();
1021  // This is primarily to skip comments that may be present.
1022  // Comments appear as nodes just like elements.
1023  // We are specifically looking for text in the body
1024  // of the definition.
1025  while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
1026  {
1027  compositeChild = compositeChild->NextSibling();
1028  }
1029 
1030  ASSERTL0(compositeChild, "Unable to read composite definition body.");
1031  std::string compositeStr = compositeChild->ToText()->ValueStr();
1032 
1033  /// Parse out the element components corresponding to type of element.
1034 
1035  std::istringstream compositeDataStrm(compositeStr.c_str());
1036 
1037  try
1038  {
1039  bool first = true;
1040  std::string prevCompositeElementStr;
1041 
1042  while (!compositeDataStrm.fail())
1043  {
1044  std::string compositeElementStr;
1045  compositeDataStrm >> compositeElementStr;
1046 
1047  if (!compositeDataStrm.fail())
1048  {
1049  if (first)
1050  {
1051  first = false;
1052 
1053  Composite curVector(MemoryManager<GeometryVector>::AllocateSharedPtr());
1054  m_meshComposites[indx] = curVector;
1055  }
1056 
1057  if (compositeElementStr.length() > 0)
1058  {
1059  ResolveGeomRef(prevCompositeElementStr, compositeElementStr, m_meshComposites[indx]);
1060  }
1061  prevCompositeElementStr = compositeElementStr;
1062  }
1063  }
1064  }
1065  catch(...)
1066  {
1067  NEKERROR(ErrorUtil::efatal,
1068  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
1069  }
1070 
1071  /// Keep looking
1072  composite = composite->NextSiblingElement("C");
1073  }
1074  }
1075 
1076 
1077  SegGeomSharedPtr MeshGraph3D::GetSegGeom(int eID)
1078  {
1079  SegGeomSharedPtr returnval;
1080  SegGeomMap::iterator x = m_segGeoms.find(eID);
1081  ASSERTL0(x != m_segGeoms.end(), "Segment "
1082  + boost::lexical_cast<string>(eID) + " not found.");
1083  return x->second;
1084  };
1085 
1086  Geometry2DSharedPtr MeshGraph3D::GetGeometry2D(int gID)
1087  {
1088  TriGeomMapIter it1;
1089  QuadGeomMapIter it2;
1090 
1091  it1 = m_triGeoms.find(gID);
1092  if (it1 != m_triGeoms.end())
1093  return it1->second;
1094 
1095  it2 = m_quadGeoms.find(gID);
1096  if (it2 != m_quadGeoms.end())
1097  return it2->second;
1098 
1099  return Geometry2DSharedPtr();
1100  };
1101 
1102  // Take the string that is the composite reference and find the
1103  // pointer to the Geometry object corresponding to it.
1104 
1105  // The only allowable combinations of previous and current items
1106  // are V (0D); E (1D); and T and Q (2D); A (Tet, 3D), P (Pyramid, 3D), R (Prism, 3D), H (Hex, 3D).
1107  // Only elements of the same dimension are allowed to be grouped.
1108  void MeshGraph3D::ResolveGeomRef(const std::string &prevToken, const std::string &token,
1109  Composite& composite)
1110  {
1111  try
1112  {
1113  std::istringstream tokenStream(token);
1114  std::istringstream prevTokenStream(prevToken);
1115 
1116  char type;
1117  char prevType;
1118 
1119  tokenStream >> type;
1120 
1121  std::string::size_type indxBeg = token.find_first_of('[') + 1;
1122  std::string::size_type indxEnd = token.find_last_of(']') - 1;
1123 
1124  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());
1125 
1126  std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
1127 
1128  std::vector<unsigned int> seqVector;
1129  std::vector<unsigned int>::iterator seqIter;
1130 
1131  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
1132 
1133  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
1134 
1135  prevTokenStream >> prevType;
1136 
1137  // All composites must be of the same dimension. This map makes things clean to compare.
1138  map<char, int> typeMap;
1139  typeMap['V'] = 1; // Vertex
1140  typeMap['E'] = 1; // Edge
1141  typeMap['T'] = 2; // Triangle
1142  typeMap['Q'] = 2; // Quad
1143  typeMap['A'] = 3; // Tet
1144  typeMap['P'] = 3; // Pyramid
1145  typeMap['R'] = 3; // Prism
1146  typeMap['H'] = 3; // Hex
1147 
1148  // Make sure only geoms of the same dimension are combined.
1149  bool validSequence = (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
1150 
1151  ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
1152  + type + " and " + prevType + ".").c_str());
1153 
1154  switch(type)
1155  {
1156  case 'V': // Vertex
1157  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1158  {
1159  if (m_vertSet.find(*seqIter) == m_vertSet.end())
1160  {
1161  char errStr[16] = "";
1162  ::sprintf(errStr, "%d", *seqIter);
1163  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
1164  }
1165  else
1166  {
1167  composite->push_back(m_vertSet[*seqIter]);
1168  }
1169  }
1170  break;
1171 
1172  case 'E': // Edge
1173  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1174  {
1175  if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
1176  {
1177  char errStr[16] = "";
1178  ::sprintf(errStr, "%d", *seqIter);
1179  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
1180  }
1181  else
1182  {
1183  composite->push_back(m_segGeoms[*seqIter]);
1184  }
1185  }
1186  break;
1187 
1188  case 'F': // Face
1189  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1190  {
1191  Geometry2DSharedPtr face = GetGeometry2D(*seqIter);
1192  if (face == Geometry2DSharedPtr())
1193  {
1194  char errStr[16] = "";
1195  ::sprintf(errStr, "%d", *seqIter);
1196  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown face index: ") + errStr).c_str());
1197  }
1198  else
1199  {
1200  if(CheckRange(*face))
1201  {
1202  composite->push_back(face);
1203  }
1204  }
1205  }
1206  break;
1207 
1208  case 'T': // Triangle
1209  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1210  {
1211  if (m_triGeoms.find(*seqIter) == m_triGeoms.end())
1212  {
1213  char errStr[16] = "";
1214  ::sprintf(errStr, "%d", *seqIter);
1215  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
1216  }
1217  else
1218  {
1219  if(CheckRange(*m_triGeoms[*seqIter]))
1220  {
1221  composite->push_back(m_triGeoms[*seqIter]);
1222  }
1223  }
1224  }
1225  break;
1226 
1227  case 'Q': // Quad
1228  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1229  {
1230  if (m_quadGeoms.find(*seqIter) == m_quadGeoms.end())
1231  {
1232  char errStr[16] = "";
1233  ::sprintf(errStr, "%d", *seqIter);
1234  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr).c_str());
1235  }
1236  else
1237  {
1238  if(CheckRange(*m_quadGeoms[*seqIter]))
1239  {
1240  composite->push_back(m_quadGeoms[*seqIter]);
1241  }
1242  }
1243  }
1244  break;
1245 
1246  // Tetrahedron
1247  case 'A':
1248  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1249  {
1250  if (m_tetGeoms.find(*seqIter) == m_tetGeoms.end())
1251  {
1252  char errStr[16] = "";
1253  ::sprintf(errStr, "%d", *seqIter);
1254  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown tet index: ") + errStr).c_str());
1255  }
1256  else
1257  {
1258  if(CheckRange(*m_tetGeoms[*seqIter]))
1259  {
1260  composite->push_back(m_tetGeoms[*seqIter]);
1261  }
1262  }
1263  }
1264  break;
1265 
1266  // Pyramid
1267  case 'P':
1268  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1269  {
1270  if (m_pyrGeoms.find(*seqIter) == m_pyrGeoms.end())
1271  {
1272  char errStr[16] = "";
1273  ::sprintf(errStr, "%d", *seqIter);
1274  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown pyramid index: ") + errStr).c_str());
1275  }
1276  else
1277  {
1278  if(CheckRange(*m_pyrGeoms[*seqIter]))
1279  {
1280  composite->push_back(m_pyrGeoms[*seqIter]);
1281  }
1282  }
1283  }
1284  break;
1285 
1286  // Prism
1287  case 'R':
1288  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1289  {
1290  if (m_prismGeoms.find(*seqIter) == m_prismGeoms.end())
1291  {
1292  char errStr[16] = "";
1293  ::sprintf(errStr, "%d", *seqIter);
1294  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown prism index: ") + errStr).c_str());
1295  }
1296  else
1297  {
1298  if(CheckRange(*m_prismGeoms[*seqIter]))
1299  {
1300  composite->push_back(m_prismGeoms[*seqIter]);
1301  }
1302  }
1303  }
1304  break;
1305 
1306  // Hex
1307  case 'H':
1308  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
1309  {
1310  if (m_hexGeoms.find(*seqIter) == m_hexGeoms.end())
1311  {
1312  char errStr[16] = "";
1313  ::sprintf(errStr, "%d", *seqIter);
1314  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown hex index: ") + errStr).c_str());
1315  }
1316  else
1317  {
1318  if(CheckRange(*m_hexGeoms[*seqIter]))
1319  {
1320  composite->push_back(m_hexGeoms[*seqIter]);
1321  }
1322  }
1323  }
1324  break;
1325 
1326  default:
1327  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
1328  }
1329  }
1330  catch(...)
1331  {
1332  NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
1333  }
1334 
1335  return;
1336  }
1337 
1338 
1339  ElementFaceVectorSharedPtr MeshGraph3D::GetElementsFromFace(Geometry2DSharedPtr face)
1340  {
1341  boost::unordered_map<int, ElementFaceVectorSharedPtr>::iterator it =
1342  m_faceToElMap.find(face->GetGlobalID());
1343 
1344  ASSERTL0(it != m_faceToElMap.end(), "Unable to find corresponding face!");
1345 
1346  return it->second;
1347  }
1348 
1349 
1350  /**
1351  * Retrieve the basis key for a given face direction.
1352  */
1353  LibUtilities::BasisKey MeshGraph3D:: GetFaceBasisKey(
1354  Geometry2DSharedPtr face,
1355  const int facedir,
1356  const std::string variable)
1357  {
1358  // Retrieve the list of elements and the associated face index
1359  // to which the face geometry belongs.
1360  ElementFaceVectorSharedPtr elements = GetElementsFromFace(face);
1361 
1362  ASSERTL0(elements->size() > 0, "No elements for the given face."
1363  " Check all elements belong to the domain composite.");
1364 
1365  // Perhaps, a check should be done here to ensure that in case
1366  // elements->size!=1, all elements to which the edge belongs have
1367  // the same type and order of expansion such that no confusion can
1368  // arise.
1369 
1370  // Get the Expansion structure detailing the basis keys used for
1371  // this element.
1372  ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element,
1373  variable);
1374 
1375  ASSERTL0(expansion, "Could not find expansion connected to face "+
1376  boost::lexical_cast<string>(face->GetGlobalID()));
1377 
1378  // Retrieve the geometry object of the element as a Geometry3D.
1379  Geometry3DSharedPtr geom3d =
1380  boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
1381  expansion->m_geomShPtr);
1382 
1383  // Use the geometry of the element to calculate the coordinate
1384  // direction of the element which corresponds to the requested
1385  // coordinate direction of the given face.
1386  int dir = geom3d->GetDir((*elements)[0]->m_FaceIndx, facedir);
1387 
1388  if(face->GetNumVerts() == 3)
1389  {
1390  return StdRegions::EvaluateTriFaceBasisKey(facedir,
1391  expansion->m_basisKeyVector[dir].GetBasisType(),
1392  expansion->m_basisKeyVector[dir].GetNumPoints(),
1393  expansion->m_basisKeyVector[dir].GetNumModes());
1394  }
1395  else
1396  {
1397  return StdRegions::EvaluateQuadFaceBasisKey(facedir,
1398  expansion->m_basisKeyVector[dir].GetBasisType(),
1399  expansion->m_basisKeyVector[dir].GetNumPoints(),
1400  expansion->m_basisKeyVector[dir].GetNumModes());
1401  }
1402 
1403  // Keep things happy by returning a value.
1404  return LibUtilities::NullBasisKey;
1405  }
1406 
1407 
1408  /**
1409  * @brief Given a 3D geometry object #element, populate the face to
1410  * element map #m_faceToElMap which maps faces to their corresponding
1411  * element(s).
1412  *
1413  * @param element Element to process.
1414  * @param kNfaces Number of faces of #element. Should be removed and
1415  * put into Geometry3D as a virtual member function.
1416  */
1417  void MeshGraph3D::PopulateFaceToElMap(Geometry3DSharedPtr element, int kNfaces)
1418  {
1419  // Set up face -> element map
1420  for (int i = 0; i < kNfaces; ++i)
1421  {
1422  int faceId = element->GetFace(i)->GetGlobalID();
1423  ElementFaceSharedPtr elementFace =
1424  MemoryManager<ElementFace>::AllocateSharedPtr();
1425 
1426  elementFace->m_Element = element;
1427  elementFace->m_FaceIndx = i;
1428 
1429  // Search map to see if face already exists.
1430  boost::unordered_map<int, ElementFaceVectorSharedPtr>::iterator it =
1431  m_faceToElMap.find(faceId);
1432 
1433  if (it == m_faceToElMap.end())
1434  {
1435  ElementFaceVectorSharedPtr tmp =
1436  MemoryManager<ElementFaceVector>::AllocateSharedPtr();
1437  tmp->push_back(elementFace);
1438  m_faceToElMap[faceId] = tmp;
1439  }
1440  else
1441  {
1442  ElementFaceVectorSharedPtr tmp = it->second;
1443  tmp->push_back(elementFace);
1444  }
1445  }
1446  }
1447  }; //end of namespace
1448 }; //end of namespace
Nektar::SpatialDomains::PyrGeomSharedPtr
boost::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:84
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
Nektar::SpatialDomains::ElementFaceVectorSharedPtr
boost::shared_ptr< ElementFaceVector > ElementFaceVectorSharedPtr
Definition: MeshGraph.h:138
Nektar::SpatialDomains::PyrGeom::kNfaces
static const int kNfaces
Definition: PyrGeom.h:60
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::SpatialDomains::MeshGraph::ReadCurves
void ReadCurves(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1158
Nektar::SpatialDomains::Geometry3D::GetDir
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
Nektar::SpatialDomains::MeshGraph3D::ReadFaces
void ReadFaces(TiXmlDocument &doc)
Definition: MeshGraph3D.cpp:234
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::SpatialDomains::HexGeom::kNtfaces
static const int kNtfaces
Definition: HexGeom.h:63
Nektar::SpatialDomains::MeshGraph3D::PopulateFaceToElMap
void PopulateFaceToElMap(Geometry3DSharedPtr element, int kNfaces)
Given a 3D geometry object #element, populate the face to element map m_faceToElMap which maps faces ...
Definition: MeshGraph3D.cpp:1417
Nektar::SpatialDomains::MeshGraph3D::ReadGeometry
void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph3D.cpp:62
Nektar::SpatialDomains::QuadGeomSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
Nektar::StdRegions::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:406
Nektar::SpatialDomains::MeshGraph::GetVertex
PointGeomSharedPtr GetVertex(int id)
Definition: MeshGraph.h:587
Nektar::SpatialDomains::MeshGraph3D::m_faceToElMap
boost::unordered_map< int, ElementFaceVectorSharedPtr > m_faceToElMap
Definition: MeshGraph3D.h:220
Nektar::SpatialDomains::HexGeom::kNfaces
static const int kNfaces
Definition: HexGeom.h:64
Nektar::SpatialDomains::QuadGeomMapIter
std::map< int, QuadGeomSharedPtr >::iterator QuadGeomMapIter
Definition: QuadGeom.h:58
Nektar::SpatialDomains::TetGeom::kNtfaces
static const int kNtfaces
Definition: TetGeom.h:59
Nektar::SpatialDomains::HexGeomSharedPtr
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Nektar::SpatialDomains::MeshGraph3D::GetGeometry2D
Geometry2DSharedPtr GetGeometry2D(int gID)
Definition: MeshGraph3D.cpp:1086
Nektar::SpatialDomains::MeshGraph3D::ResolveGeomRef
void ResolveGeomRef(const std::string &prevToken, const std::string &token, Composite &composite)
Definition: MeshGraph3D.cpp:1108
Nektar::SpatialDomains::TetGeom::kNfaces
static const int kNfaces
Definition: TetGeom.h:60
Nektar::SpatialDomains::MeshGraph3D::GetSegGeom
SegGeomSharedPtr GetSegGeom(int eID)
Definition: MeshGraph3D.cpp:1077
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::ParseUtils::GenerateSeqVector
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
Nektar::SpatialDomains::MeshGraph::ReadGeometry
virtual void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph.cpp:230
Nektar::SpatialDomains::SegGeom::GetEdgeOrientation
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:293
Nektar::SpatialDomains::TriGeomMapIter
std::map< int, TriGeomSharedPtr >::iterator TriGeomMapIter
Definition: TriGeom.h:63
Nektar::SpatialDomains::MeshGraph3D::GetElementsFromFace
ElementFaceVectorSharedPtr GetElementsFromFace(Geometry2DSharedPtr face)
Return the elements (shared ptrs) that have this face.
Definition: MeshGraph3D.cpp:1339
Nektar::SpatialDomains::MeshGraph3D::ReadComposites
void ReadComposites(TiXmlDocument &doc)
Definition: MeshGraph3D.cpp:982
Nektar::StdRegions::EvaluateQuadFaceBasisKey
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:349
Nektar::SpatialDomains::SegGeomSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Nektar::SpatialDomains::DomainRangeShPtr
boost::shared_ptr< DomainRange > DomainRangeShPtr
Definition: MeshGraph.h:157
Nektar::SpatialDomains::MeshGraph3D::ReadElements
void ReadElements(TiXmlDocument &doc)
Definition: MeshGraph3D.cpp:507
MeshGraph3D.h
TriGeom.h
Nektar::SpatialDomains::Composite
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
Nektar::SpatialDomains::ElementFaceSharedPtr
boost::shared_ptr< ElementFace > ElementFaceSharedPtr
Definition: MeshGraph.h:136
Nektar::SpatialDomains::MeshGraph::m_compositesLabels
map< int, string > m_compositesLabels
Definition: MeshGraph.h:431
Nektar::SpatialDomains::TetGeom::kNqfaces
static const int kNqfaces
Definition: TetGeom.h:58
Nektar::SpatialDomains::TriGeom::kNedges
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:98
Nektar::SpatialDomains::PyrGeom::kNtfaces
static const int kNtfaces
Definition: PyrGeom.h:59
Nektar::SpatialDomains::Geometry2DSharedPtr
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::SpatialDomains::MeshGraph3D::GetFaceBasisKey
LibUtilities::BasisKey GetFaceBasisKey(Geometry2DSharedPtr face, const int flag, const std::string variable="DefaultVar")
Return the BasisKey corresponding to a face of an element.
Definition: MeshGraph3D.cpp:1353
Nektar::SpatialDomains::MeshGraph3D::ReadEdges
void ReadEdges(TiXmlDocument &doc)
Definition: MeshGraph3D.cpp:98
Nektar::SpatialDomains::PyrGeom::kNqfaces
static const int kNqfaces
Definition: PyrGeom.h:58
Nektar::SpatialDomains::Geometry3D
3D geometry information
Definition: Geometry3D.h:70
Nektar::SpatialDomains::ExpansionShPtr
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:173
Nektar::SpatialDomains::MeshGraph
Base class for a spectral/hp element mesh.
Definition: MeshGraph.h:186
Nektar::SpatialDomains::MeshGraph::ReadDomain
void ReadDomain(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1064
Nektar::SpatialDomains::PrismGeomSharedPtr
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
Nektar::SpatialDomains::MeshGraph::CheckRange
bool CheckRange(Geometry2D &geom)
Check if goemetry is in range definition if activated.
Definition: MeshGraph.cpp:2000
Nektar::SpatialDomains::MeshGraph::GetExpansion
ExpansionShPtr GetExpansion(GeometrySharedPtr geom, const std::string variable="DefaultVar")
Definition: MeshGraph.cpp:2323
Nektar::SpatialDomains::TetGeomSharedPtr
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
Nektar::SpatialDomains::TriGeomSharedPtr
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::SpatialDomains::MeshGraph::ReadExpansions
void ReadExpansions(const std::string &infilename)
Read the expansions given the XML file path.
Definition: MeshGraph.cpp:580
Nektar::SpatialDomains::HexGeom::kNqfaces
static const int kNqfaces
Definition: HexGeom.h:62
ParseUtils.hpp
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::SpatialDomains::PointGeomSharedPtr
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
Nektar::LibUtilities::CompressData::ZlibDecodeFromBase64Str
int ZlibDecodeFromBase64Str(std::string &in64, std::vector< T > &out)
Definition: CompressData.h:243
Nektar::SpatialDomains::Geometry3DSharedPtr
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52
Generated on Wed Mar 23 2016 20:14:21 for Nektar++ by   doxygen 1.8.8