Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MeshGraph2D.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////////
2 // File: MeshGraph2D.cpp
3 //
4 // For more information, please see: http://www.nektar.info/
5 //
6 // The MIT License
7 //
8 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
9 // Department of Aeronautics, Imperial College London (UK), and Scientific
10 // Computing and Imaging Institute, University of Utah (USA).
11 //
12 // License for the specific language governing rights and limitations under
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description:
32 //
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 #include <SpatialDomains/SegGeom.h>
38 #include <SpatialDomains/TriGeom.h>
40 #include <tinyxml.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace SpatialDomains
47  {
48  MeshGraph2D::MeshGraph2D()
49  {
50  }
51 
52  MeshGraph2D::~MeshGraph2D()
53  {
54  }
55 
56  MeshGraph2D::MeshGraph2D(const LibUtilities::SessionReaderSharedPtr &pSession,
57  const DomainRangeShPtr &rng)
58  : MeshGraph(pSession,rng)
59  {
60  ReadGeometry(pSession->GetDocument());
61  ReadExpansions(pSession->GetDocument());
62  }
63 
64  void MeshGraph2D::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 MeshGraph2D::ReadGeometry(TiXmlDocument &doc)
80  {
81  // Read mesh first
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  ReadElements(doc);
95  ReadComposites(doc);
96  ReadDomain(doc);
97  }
98 
99  void MeshGraph2D::ReadEdges(TiXmlDocument &doc)
100  {
101  /// We know we have it since we made it this far.
102  TiXmlHandle docHandle(&doc);
103  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
104  TiXmlElement* field = NULL;
106 
107  /// Look for elements in ELEMENT block.
108  field = mesh->FirstChildElement("EDGE");
109 
110  ASSERTL0(field, "Unable to find EDGE tag in file.");
111 
112  string IsCompressed;
113  field->QueryStringAttribute("COMPRESSED",&IsCompressed);
114 
115  if(IsCompressed.size())
116  {
117  ASSERTL0(boost::iequals(IsCompressed,
119  "Compressed formats do not match. Expected :"
121  + " but got " + std::string(IsCompressed));
122 
123  // Extract the edge body
124  TiXmlNode* edgeChild = field->FirstChild();
125  ASSERTL0(edgeChild, "Unable to extract the data from "
126  "the compressed edge tag.");
127 
128  std::string edgeStr;
129  if (edgeChild->Type() == TiXmlNode::TINYXML_TEXT)
130  {
131  edgeStr += edgeChild->ToText()->ValueStr();
132  }
133 
134  std::vector<LibUtilities::MeshEdge> edgeData;
136  edgeData);
137  int indx;
138  for(int i = 0; i < edgeData.size(); ++i)
139  {
140  indx = edgeData[i].id;
141 
142  PointGeomSharedPtr vertices[2] =
143  {GetVertex(edgeData[i].v0), GetVertex(edgeData[i].v1)};
144  SegGeomSharedPtr edge;
145 
146  it = m_curvedEdges.find(indx);
147  if (it == m_curvedEdges.end())
148  {
150  indx, m_spaceDimension, vertices);
151  }
152  else
153  {
155  indx, m_spaceDimension, vertices,
156  it->second);
157  }
158  m_segGeoms[indx] = edge;
159  }
160  }
161  else
162  {
163  /// All elements are of the form: "<E ID="#"> ... </E>", with
164  /// ? being the element type.
165  /// Read the ID field first.
166  TiXmlElement *edge = field->FirstChildElement("E");
167 
168  /// Since all edge data is one big text block, we need to
169  /// accumulate all TINYXML_TEXT data and then parse it. This
170  /// approach effectively skips all comments or other node
171  /// types since we only care about the edge list. We
172  /// cannot handle missing edge numbers as we could with
173  /// missing element numbers due to the text block format.
174  std::string edgeStr;
175  int indx;
176 
177  while(edge)
178  {
179  int err = edge->QueryIntAttribute("ID",&indx);
180  ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
181 
182  TiXmlNode *child = edge->FirstChild();
183  edgeStr.clear();
184  if (child->Type() == TiXmlNode::TINYXML_TEXT)
185  {
186  edgeStr += child->ToText()->ValueStr();
187  }
188 
189  /// Now parse out the edges, three fields at a time.
190  int vertex1, vertex2;
191  std::istringstream edgeDataStrm(edgeStr.c_str());
192 
193  try
194  {
195  while (!edgeDataStrm.fail())
196  {
197  edgeDataStrm >> vertex1 >> vertex2;
198 
199  // Must check after the read because we may be
200  // at the end and not know it. If we are at
201  // the end we will add a duplicate of the last
202  // entry if we don't check here.
203  if (!edgeDataStrm.fail())
204  {
205  PointGeomSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
206 
207  SegGeomSharedPtr edge;
208 
209  it = m_curvedEdges.find(indx);
210 
211  if(it == m_curvedEdges.end())
212  {
214  edge->SetGlobalID(indx); // Set global mesh id
215  }
216  else
217  {
218  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices, it->second);
219  edge->SetGlobalID(indx); //Set global mesh id
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 MeshGraph2D::ReadElements(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 ELEMENT block.
244  field = mesh->FirstChildElement("ELEMENT");
245 
246  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
247 
248  // Set up curve map for curved elements on an embedded manifold.
250 
251  /// All elements are of the form: "<? ID="#"> ... </?>", with
252  /// ? being the element type.
253 
254  TiXmlElement *element = field->FirstChildElement();
255 
256  while (element)
257  {
258  std::string elementType(element->ValueStr());
259 
260  ASSERTL0(elementType == "Q" || elementType == "T",
261  (std::string("Unknown 2D element type: ") + elementType).c_str());
262 
263  string IsCompressed;
264  element->QueryStringAttribute("COMPRESSED",&IsCompressed);
265 
266  if(IsCompressed.size())
267  {
268  ASSERTL0(boost::iequals(IsCompressed,
270  "Compressed formats do not match. Expected :"
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;
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.
301  {
302  GetSegGeom(faceData[i].e[0]),
303  GetSegGeom(faceData[i].e[1]),
304  GetSegGeom(faceData[i].e[2])
305  };
306 
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 =
319  indx, edges, edgeorient);
320  }
321  else
322  {
323  trigeom =
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;
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.
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 
353  {
354  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
355  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
356  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
357  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
358  };
359 
360  QuadGeomSharedPtr quadgeom;
361  if (it == m_curvedEdges.end())
362  {
363  quadgeom =
365  indx, edges, edgeorient);
366  }
367  else
368  {
369  quadgeom =
371  indx, edges, edgeorient,
372  it->second);
373  }
374  quadgeom->SetGlobalID(indx);
375  m_quadGeoms[indx] = quadgeom;
376  }
377  }
378  }
379  else
380  {
381  /// Read id attribute.
382  int indx;
383  int err = element->QueryIntAttribute("ID", &indx);
384  ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
385 
386  it = m_curvedFaces.find(indx);
387 
388  /// Read text element description.
389  TiXmlNode* elementChild = element->FirstChild();
390  std::string elementStr;
391  while(elementChild)
392  {
393  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
394  {
395  elementStr += elementChild->ToText()->ValueStr();
396  }
397  elementChild = elementChild->NextSibling();
398  }
399 
400  ASSERTL0(!elementStr.empty(), "Unable to read element description body.");
401 
402  /// Parse out the element components corresponding to type of element.
403  if (elementType == "T")
404  {
405  // Read three edge numbers
406  int edge1, edge2, edge3;
407  std::istringstream elementDataStrm(elementStr.c_str());
408 
409  try
410  {
411  elementDataStrm >> edge1;
412  elementDataStrm >> edge2;
413  elementDataStrm >> edge3;
414 
415  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
416 
417  /// Create a TriGeom to hold the new definition.
419  {
420  GetSegGeom(edge1),
421  GetSegGeom(edge2),
422  GetSegGeom(edge3)
423  };
424 
426  {
427  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
428  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
429  SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
430  };
431 
432  TriGeomSharedPtr trigeom;
433  if (it == m_curvedFaces.end())
434  {
435  trigeom = MemoryManager<TriGeom>
437  edges,
438  edgeorient);
439  }
440  else
441  {
442  trigeom = MemoryManager<TriGeom>
444  edges,
445  edgeorient,
446  it->second);
447  }
448  trigeom->SetGlobalID(indx);
449 
450  m_triGeoms[indx] = trigeom;
451  }
452  catch(...)
453  {
454  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
455  }
456  }
457  else if (elementType == "Q")
458  {
459  // Read four edge numbers
460  int edge1, edge2, edge3, edge4;
461  std::istringstream elementDataStrm(elementStr.c_str());
462 
463  try
464  {
465  elementDataStrm >> edge1;
466  elementDataStrm >> edge2;
467  elementDataStrm >> edge3;
468  elementDataStrm >> edge4;
469 
470  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
471 
472  /// Create a QuadGeom to hold the new definition.
474  {GetSegGeom(edge1),GetSegGeom(edge2),
475  GetSegGeom(edge3),GetSegGeom(edge4)};
476 
478  {
479  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
480  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
481  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
482  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
483  };
484 
485  QuadGeomSharedPtr quadgeom;
486  if (it == m_curvedFaces.end())
487  {
488  quadgeom = MemoryManager<QuadGeom>
490  edges,
491  edgeorient);
492  }
493  else
494  {
495  quadgeom = MemoryManager<QuadGeom>
497  edges,
498  edgeorient,
499  it->second);
500  }
501  quadgeom->SetGlobalID(indx);
502 
503  m_quadGeoms[indx] = quadgeom;
504 
505  }
506  catch(...)
507  {
508  NEKERROR(ErrorUtil::efatal,(std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
509  }
510  }
511  }
512  /// Keep looking
513  element = element->NextSiblingElement();
514  }
515  }
516 
517  void MeshGraph2D::ReadComposites(TiXmlDocument &doc)
518  {
519  TiXmlHandle docHandle(&doc);
520 
521  /// We know we have it since we made it this far.
522  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
523  TiXmlElement* field = NULL;
524 
525  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
526 
527  /// Look for elements in ELEMENT block.
528  field = mesh->FirstChildElement("COMPOSITE");
529 
530  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
531 
532  int nextCompositeNumber = -1;
533 
534  /// All elements are of the form: "<C ID = "N"> ... </C>".
535  /// Read the ID field first.
536  TiXmlElement *composite = field->FirstChildElement("C");
537 
538  while (composite)
539  {
540  nextCompositeNumber++;
541 
542  int indx;
543  int err = composite->QueryIntAttribute("ID", &indx);
544  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
545 
546  // read and store label if they exist
547  string labelstr;
548  err = composite->QueryStringAttribute("LABEL", &labelstr);
549  if(err == TIXML_SUCCESS)
550  {
551  m_compositesLabels[indx] = labelstr;
552  }
553 
554  TiXmlNode* compositeChild = composite->FirstChild();
555  // This is primarily to skip comments that may be present.
556  // Comments appear as nodes just like elements.
557  // We are specifically looking for text in the body
558  // of the definition.
559  while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
560  {
561  compositeChild = compositeChild->NextSibling();
562  }
563 
564  ASSERTL0(compositeChild, "Unable to read composite definition body.");
565  std::string compositeStr = compositeChild->ToText()->ValueStr();
566 
567  /// Parse out the element components corresponding to type of element.
568 
569  std::istringstream compositeDataStrm(compositeStr.c_str());
570 
571  try
572  {
573  bool first = true;
574  std::string prevCompositeElementStr;
575 
576  while (!compositeDataStrm.fail())
577  {
578  std::string compositeElementStr;
579  compositeDataStrm >> compositeElementStr;
580 
581  if (!compositeDataStrm.fail())
582  {
583  if (first)
584  {
585  first = false;
586 
588  m_meshComposites[indx] = curVector;
589  }
590 
591  if (compositeElementStr.length() > 0)
592  {
593  ResolveGeomRef(prevCompositeElementStr, compositeElementStr, m_meshComposites[indx]);
594  }
595  prevCompositeElementStr = compositeElementStr;
596  }
597  }
598  }
599  catch(...)
600  {
602  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
603  }
604 
605  /// Keep looking
606  composite = composite->NextSiblingElement("C");
607  }
608  }
609 
610 
612  {
613  SegGeomSharedPtr returnval;
614  SegGeomMap::iterator x = m_segGeoms.find(eID);
615  ASSERTL0(x != m_segGeoms.end(), "Segment not found.");
616  return x->second;
617  };
618 
619 
620  // Take the string that is the composite reference and find the
621  // pointer to the Geometry object corresponding to it.
622 
623  // The only allowable combinations of previous and current items
624  // are V (0D); E (1D); and T and Q (2D). Only elements of the same
625  // dimension are allowed to be grouped.
626  void MeshGraph2D::ResolveGeomRef(const std::string &prevToken, const std::string &token,
627  Composite& composite)
628  {
629  try
630  {
631  std::istringstream tokenStream(token);
632  std::istringstream prevTokenStream(prevToken);
633 
634  char type;
635  char prevType;
636 
637  tokenStream >> type;
638 
639  std::string::size_type indxBeg = token.find_first_of('[') + 1;
640  std::string::size_type indxEnd = token.find_last_of(']') - 1;
641 
642  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());
643 
644  std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
645  std::vector<unsigned int> seqVector;
647 
648  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
649 
650  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
651 
652  prevTokenStream >> prevType;
653 
654  // All composites must be of the same dimension.
655  bool validSequence = (prevToken.empty() || // No previous, then current is just fine.
656  (type == 'V' && prevType == 'V') ||
657  (type == 'E' && prevType == 'E') ||
658  ((type == 'T' || type == 'Q') &&
659  (prevType == 'T' || prevType == 'Q')));
660 
661  ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
662  + type + " and " + prevType + ".").c_str());
663 
664 
665  switch(type)
666  {
667  case 'E': // Edge
668  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
669  {
670  if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
671  {
672  char errStr[16] = "";
673  ::sprintf(errStr, "%d", *seqIter);
674  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
675  }
676  else
677  {
678  composite->push_back(m_segGeoms[*seqIter]);
679  }
680  }
681  break;
682 
683  case 'T': // Triangle
684  {
685  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
686  {
687  if (m_triGeoms.count(*seqIter) == 0 )
688  {
689  char errStr[16] = "";
690  ::sprintf(errStr, "%d", *seqIter);
691  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
692  }
693  else
694  {
695  if(CheckRange(*m_triGeoms[*seqIter]))
696  {
697  composite->push_back(m_triGeoms[*seqIter]);
698  }
699  }
700  }
701  }
702  break;
703 
704  case 'Q': // Quad
705  {
706  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
707  {
708  if (m_quadGeoms.count(*seqIter) == 0)
709  {
710  char errStr[16] = "";
711  ::sprintf(errStr, "%d", *seqIter);
712  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr +std::string(" in Composite section")).c_str());
713  }
714  else
715  {
716  if(CheckRange(*m_quadGeoms[*seqIter]))
717  {
718  composite->push_back(m_quadGeoms[*seqIter]);
719  }
720  }
721  }
722  }
723  break;
724 
725  case 'V': // Vertex
726  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
727  {
728  if (*seqIter >= m_vertSet.size())
729  {
730  char errStr[16] = "";
731  ::sprintf(errStr, "%d", *seqIter);
732  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
733  }
734  else
735  {
736  composite->push_back(m_vertSet[*seqIter]);
737  }
738  }
739  break;
740 
741  default:
742  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
743  }
744  }
745  catch(...)
746  {
747  NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
748  }
749 
750  return;
751  }
752 
754  {
755  SegGeomSharedPtr Sedge;
756 
757  if(!(Sedge = boost::dynamic_pointer_cast<SegGeom>(edge)))
758  {
759  ASSERTL0(false,"Dynamics cast failed");
760 
761  }
762  return GetElementsFromEdge(Sedge);
763 
764  }
766  {
767  // Search tris and quads
768  // Need to iterate through vectors because there may be multiple
769  // occurrences.
770  ElementEdgeSharedPtr elementEdge;
771  //TriGeomVectorIter triIter;
772 
774 
775  CompositeMapIter compIter;
776  TriGeomSharedPtr triGeomShPtr;
777  QuadGeomSharedPtr quadGeomShPtr;
778 
779  GeometryVectorIter geomIter;
780 
781  for(int d = 0; d < m_domain.size(); ++d)
782  {
783  for (compIter = m_domain[d].begin(); compIter != m_domain[d].end(); ++compIter)
784  {
785  for (geomIter = (compIter->second)->begin(); geomIter != (compIter->second)->end(); ++geomIter)
786  {
787  triGeomShPtr = boost::dynamic_pointer_cast<TriGeom>(*geomIter);
788  quadGeomShPtr = boost::dynamic_pointer_cast<QuadGeom>(*geomIter);
789 
790  if (triGeomShPtr || quadGeomShPtr)
791  {
792  int edgeNum;
793  if (triGeomShPtr)
794  {
795  if ((edgeNum = triGeomShPtr->WhichEdge(edge)) > -1)
796  {
798  elementEdge->m_Element = triGeomShPtr;
799  elementEdge->m_EdgeIndx = edgeNum;
800  returnval->push_back(elementEdge);
801  }
802  }
803  else if (quadGeomShPtr)
804  {
805  if ((edgeNum = quadGeomShPtr->WhichEdge(edge)) > -1)
806  {
808  elementEdge->m_Element = quadGeomShPtr;
809  elementEdge->m_EdgeIndx = edgeNum;
810  returnval->push_back(elementEdge);
811  }
812  }
813  }
814  }
815  }
816  }
817 
818  return returnval;
819  }
820 
822  {
824  // Perhaps, a check should be done here to ensure that
825  // in case elements->size!=1, all elements to which
826  // the edge belongs have the same type and order of
827  // expansion such that no confusion can arise.
828  ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element, variable);
829 
830  int edge_id = (*elements)[0]->m_EdgeIndx;
831 
832  if((*elements)[0]->m_Element->GetShapeType() == LibUtilities::eTriangle)
833  {
834  edge_id = (edge_id)? 1:0;
835  }
836  else
837  {
838  edge_id = edge_id%2;
839  }
840 
841  int nummodes = expansion->m_basisKeyVector[edge_id].GetNumModes();
842  int numpoints = expansion->m_basisKeyVector[edge_id].GetNumPoints();
843 
844  if((*elements)[0]->m_Element->GetShapeType() == LibUtilities::eTriangle)
845  {
846  // Use edge 0 to define basis of order relevant to edge
847  switch(expansion->m_basisKeyVector[edge_id].GetBasisType())
848  {
850  {
851  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
852  {
854  {
856  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
857  }
858  break;
859 
860  default:
861  ASSERTL0(false,"Unexpected points distribution");
862 
863  // It doesn't matter what we return
864  // here since the ASSERT will stop
865  // execution. Just return something
866  // to prevent warnings messages.
868  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
869  break;
870  }
871  }
872  break;
873  case LibUtilities::eOrtho_B: // Assume this is called from nodal triangular basis
874  {
875  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
876  {
878  {
881  }
882  break;
883 
884  default:
885  ASSERTL0(false,"Unexpected points distribution");
886 
887  // It doesn't matter what we return
888  // here since the ASSERT will stop
889  // execution. Just return something
890  // to prevent warnings messages.
892  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
893  break;
894  }
895  }
896  break;
898  {
899  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
900  {
902  {
904  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
905  }
906  break;
908  {
910  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
911  }
912  break;
913 
914  default:
915  ASSERTL0(false,"Unexpected points distribution");
916 
917  // It doesn't matter what we return
918  // here since the ASSERT will stop
919  // execution. Just return something
920  // to prevent warnings messages.
922  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
923  break;
924  }
925  }
926  break;
928  {
929  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
930  {
932  {
934  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
935  }
936  break;
937  default:
938  ASSERTL0(false,"Unexpected points distribution");
939  // It doesn't matter what we return here
940  // since the ASSERT will stop execution.
941  // Just return something to prevent
942  // warnings messages.
944  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
945  break;
946  }
947  }
948  break;
949 
950  default:
951  ASSERTL0(false,"Unexpected basis distribution");
952  // It doesn't matter what we return here since the
953  // ASSERT will stop execution. Just return
954  // something to prevent warnings messages.
956  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
957  }
958  }
959  else
960  {
961  // Quadrilateral
962  const LibUtilities::PointsKey pkey(numpoints,expansion->m_basisKeyVector[edge_id].GetPointsType());
963  return LibUtilities::BasisKey(expansion->m_basisKeyVector[edge_id].GetBasisType(),nummodes,pkey);
964  }
965 
966  ASSERTL0(false, "Unable to determine edge points type.");
968  }
969  }; //end of namespace
970 }; //end of namespace
971 
972 //
973 // $Log: MeshGraph2D.cpp,v $
974 // Revision 1.41 2009/12/17 02:08:04 bnelson
975 // Fixed visual studio compiler warning.
976 //
977 // Revision 1.40 2009/10/22 17:41:47 cbiotto
978 // Update for variable order expansion
979 //
980 // Revision 1.39 2009/07/02 13:32:24 sehunchun
981 // *** empty log message ***
982 //
983 // Revision 1.38 2009/07/02 11:39:44 sehunchun
984 // Modification for 2D gemoetry embedded in 3D
985 //
986 // Revision 1.37 2009/05/01 13:23:21 pvos
987 // Fixed various bugs
988 //
989 // Revision 1.36 2009/04/20 16:13:23 sherwin
990 // Modified Import and Write functions and redefined how Expansion is used
991 //
992 // Revision 1.35 2009/01/21 16:59:03 pvos
993 // Added additional geometric factors to improve efficiency
994 //
995 // Revision 1.34 2009/01/12 10:26:59 pvos
996 // Added input tags for nodal expansions
997 //
998 // Revision 1.33 2008/09/12 11:26:19 pvos
999 // Updates for mappings in 3D
1000 //
1001 // Revision 1.32 2008/09/09 14:21:44 sherwin
1002 // Updates for first working version of curved edges
1003 //
1004 // Revision 1.31 2008/08/14 22:11:03 sherwin
1005 // Mods for HDG update
1006 //
1007 // Revision 1.30 2008/07/29 22:23:36 sherwin
1008 // various mods for DG advection solver in Multiregions. Added virtual calls to Geometry, Geometry1D, 2D and 3D
1009 //
1010 // Revision 1.29 2008/05/29 19:07:39 delisi
1011 // Removed the Write(..) methods, so it is only in the base MeshGraph class. Also, added a line to set the global ID of the geometry object for every element read in.
1012 //
1013 // Revision 1.28 2008/05/28 21:42:18 jfrazier
1014 // Minor comment spelling change.
1015 //
1016 // Revision 1.27 2008/03/18 14:14:49 pvos
1017 // Update for nodal triangular helmholtz solver
1018 //
1019 // Revision 1.26 2008/03/11 17:02:24 jfrazier
1020 // Modified the GetElementsFromEdge to use the domain.
1021 //
1022 // Revision 1.25 2008/01/21 19:58:14 sherwin
1023 // Updated so that QuadGeom and TriGeom have SegGeoms instead of EdgeComponents
1024 //
1025 // Revision 1.24 2007/12/17 20:27:24 sherwin
1026 // Added normals to GeomFactors
1027 //
1028 // Revision 1.23 2007/12/14 17:39:18 jfrazier
1029 // Fixed composite reader to handle ranges and comma separated lists.
1030 //
1031 // Revision 1.22 2007/12/11 21:51:53 jfrazier
1032 // Updated 2d components so elements could be retrieved from edges.
1033 //
1034 // Revision 1.21 2007/12/04 03:02:26 jfrazier
1035 // Changed to stringstream.
1036 //
1037 // Revision 1.20 2007/09/20 22:25:06 jfrazier
1038 // Added expansion information to meshgraph class.
1039 //
1040 // Revision 1.19 2007/07/26 01:38:33 jfrazier
1041 // Cleanup of some attribute reading code.
1042 //
1043 // Revision 1.18 2007/07/24 16:52:09 jfrazier
1044 // Added domain code.
1045 //
1046 // Revision 1.17 2007/07/05 04:21:10 jfrazier
1047 // Changed id format and propagated from 1d to 2d.
1048 //
1049 // Revision 1.16 2007/06/10 02:27:11 jfrazier
1050 // Another checkin with an incremental completion of the boundary conditions reader
1051 //
1052 // Revision 1.15 2007/06/07 23:55:24 jfrazier
1053 // Intermediate revisions to add parsing for boundary conditions file.
1054 //
1055 // Revision 1.14 2007/05/28 21:48:42 sherwin
1056 // Update for 2D functionality
1057 //
1058 // Revision 1.13 2006/10/17 22:26:01 jfrazier
1059 // Added capability to specify ranges in composite definition.
1060 //
1061 // Revision 1.12 2006/10/17 18:42:54 jfrazier
1062 // Removed "NUMBER" attribute in items.
1063 //
1064 // Revision 1.11 2006/09/26 23:41:53 jfrazier
1065 // Updated to account for highest level NEKTAR tag and changed the geometry tag to GEOMETRY.
1066 //
1067 // Revision 1.10 2006/08/24 18:50:00 jfrazier
1068 // Completed error checking on permissable composite item combinations.
1069 //
1070 // Revision 1.9 2006/08/18 19:37:17 jfrazier
1071 // *** empty log message ***
1072 //
1073 // Revision 1.8 2006/08/17 22:55:00 jfrazier
1074 // Continued adding code to process composites in the mesh2d.
1075 //
1076 // Revision 1.7 2006/08/16 23:34:42 jfrazier
1077 // *** empty log message ***
1078 //
1079 // Revision 1.6 2006/06/02 18:48:40 sherwin
1080 // Modifications to make ProjectLoc2D run bit there are bus errors for order > 3
1081 //
1082 // Revision 1.5 2006/06/01 14:15:30 sherwin
1083 // Added typdef of boost wrappers and made GeoFac a boost shared pointer.
1084 //
1085 // Revision 1.4 2006/05/30 14:00:04 sherwin
1086 // Updates to make MultiRegions and its Demos work
1087 //
1088 // Revision 1.3 2006/05/23 19:56:33 jfrazier
1089 // These build and run, but the expansion pieces are commented out
1090 // because they would not run.
1091 //
1092 // Revision 1.2 2006/05/09 13:37:01 jfrazier
1093 // Removed duplicate definition of shared vertex pointer.
1094 //
1095 // Revision 1.1 2006/05/04 18:59:01 kirby
1096 // *** empty log message ***
1097 //
1098 // Revision 1.11 2006/04/11 23:18:11 jfrazier
1099 // Completed MeshGraph2D for tri's and quads. Not thoroughly tested.
1100 //
1101 // Revision 1.10 2006/04/09 02:08:35 jfrazier
1102 // Added precompiled header.
1103 //
1104 // Revision 1.9 2006/04/04 23:12:37 jfrazier
1105 // More updates to readers. Still work to do on MeshGraph2D to store tris and quads.
1106 //
1107 // Revision 1.8 2006/03/25 00:58:29 jfrazier
1108 // Many changes dealing with fundamental structure and reading/writing.
1109 //
1110 // Revision 1.7 2006/03/12 14:20:43 sherwin
1111 //
1112 // First compiling version of SpatialDomains and associated modifications
1113 //
1114 // Revision 1.6 2006/03/12 07:42:03 sherwin
1115 //
1116 // Updated member names and StdRegions call. Still has not been compiled
1117 //
1118 // Revision 1.5 2006/02/26 21:19:43 bnelson
1119 // Fixed a variety of compiler errors caused by updates to the coding standard.
1120 //
1121 // Revision 1.4 2006/02/19 01:37:34 jfrazier
1122 // Initial attempt at bringing into conformance with the coding standard. Still more work to be done. Has not been compiled.
1123 //
1124 //
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
void ReadEdges(TiXmlDocument &doc)
Definition: MeshGraph2D.cpp:99
void ReadComposites(TiXmlDocument &doc)
std::map< int, std::string > m_compositesLabels
Definition: MeshGraph.h:431
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void ReadCurves(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1160
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
void ReadElements(TiXmlDocument &doc)
LibUtilities::BasisKey GetEdgeBasisKey(SegGeomSharedPtr edge, const std::string variable="DefaultVar")
Return the BasisKey corresponding to an edge of an element If the expansion is a triangle the Modifie...
Principle Modified Functions .
Definition: BasisType.h:49
PointGeomSharedPtr GetVertex(int id)
Definition: MeshGraph.h:587
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
virtual void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph.cpp:232
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:293
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
std::vector< GeometrySharedPtr >::iterator GeometryVectorIter
Definition: Geometry.h:58
Principle Orthogonal Functions .
Definition: BasisType.h:47
boost::shared_ptr< ElementEdge > ElementEdgeSharedPtr
Definition: MeshGraph.h:132
void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph2D.cpp:64
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Principle Modified Functions .
Definition: BasisType.h:50
boost::shared_ptr< DomainRange > DomainRangeShPtr
Definition: MeshGraph.h:157
std::map< int, Composite >::iterator CompositeMapIter
Definition: MeshGraph.h:116
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
ElementEdgeVectorSharedPtr GetElementsFromEdge(SegGeomSharedPtr edge)
Return the elements (shared ptrs) that have this edge.
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:98
boost::shared_ptr< ElementEdgeVector > ElementEdgeVectorSharedPtr
Definition: MeshGraph.h:134
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::vector< CompositeMap > m_domain
Definition: MeshGraph.h:432
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:173
Base class for a spectral/hp element mesh.
Definition: MeshGraph.h:186
void ReadDomain(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1066
bool CheckRange(Geometry2D &geom)
Check if goemetry is in range definition if activated.
Definition: MeshGraph.cpp:2002
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
ExpansionShPtr GetExpansion(GeometrySharedPtr geom, const std::string variable="DefaultVar")
Definition: MeshGraph.cpp:2325
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
void ReadExpansions(const std::string &infilename)
Read the expansions given the XML file path.
Definition: MeshGraph.cpp:582
Lagrange for SEM basis .
Definition: BasisType.h:53
void ResolveGeomRef(const std::string &prevToken, const std::string &token, Composite &composite)
Describes the specification for a Basis.
Definition: Basis.h:50
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
int ZlibDecodeFromBase64Str(std::string &in64, std::vector< T > &out)
Definition: CompressData.h:243
SegGeomSharedPtr GetSegGeom(int eID)