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