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