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 
36 #include <SpatialDomains/MeshGraph2D.h>
37 #include <SpatialDomains/SegGeom.h>
38 #include <SpatialDomains/TriGeom.h>
39 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
40 #include <tinyxml.h>
41 
42 namespace Nektar
43 {
44  namespace SpatialDomains
45  {
46  MeshGraph2D::MeshGraph2D()
47  {
48  }
49 
50  MeshGraph2D::~MeshGraph2D()
51  {
52  }
53 
54  MeshGraph2D::MeshGraph2D(const LibUtilities::SessionReaderSharedPtr &pSession,
55  const DomainRangeShPtr &rng)
56  : MeshGraph(pSession,rng)
57  {
58  ReadGeometry(pSession->GetDocument());
59  ReadExpansions(pSession->GetDocument());
60  }
61 
62  void MeshGraph2D::ReadGeometry(const std::string &infilename)
63  {
64  TiXmlDocument doc(infilename);
65  bool loadOkay = doc.LoadFile();
66 
67  std::stringstream errstr;
68  errstr << "Unable to load file: " << infilename << "\n";
69  errstr << doc.ErrorDesc() << " (Line " << doc.ErrorRow()
70  << ", Column " << doc.ErrorCol() << ")";
71  ASSERTL0(loadOkay, errstr.str());
72 
73  ReadGeometry(doc);
74  }
75 
76  // \brief Read segments (and general MeshGraph) given TiXmlDocument.
77  void MeshGraph2D::ReadGeometry(TiXmlDocument &doc)
78  {
79  // Read mesh first
80  MeshGraph::ReadGeometry(doc);
81  TiXmlHandle docHandle(&doc);
82 
83  TiXmlElement* mesh = NULL;
84 
85  /// Look for all geometry related data in GEOMETRY block.
86  mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
87 
88  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
89 
90  ReadCurves(doc);
91  ReadEdges(doc);
92  ReadElements(doc);
93  ReadComposites(doc);
94  ReadDomain(doc);
95  }
96 
97  void MeshGraph2D::ReadEdges(TiXmlDocument &doc)
98  {
99  /// We know we have it since we made it this far.
100  TiXmlHandle docHandle(&doc);
101  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
102  TiXmlElement* field = NULL;
103 
104  /// Look for elements in ELEMENT block.
105  field = mesh->FirstChildElement("EDGE");
106 
107  ASSERTL0(field, "Unable to find EDGE tag in file.");
108 
109  /// All elements are of the form: "<E ID="#"> ... </E>", with
110  /// ? being the element type.
111  /// Read the ID field first.
112  TiXmlElement *edge = field->FirstChildElement("E");
113 
114  /// Since all edge data is one big text block, we need to
115  /// accumulate all TEXT data and then parse it. This
116  /// approach effectively skips all comments or other node
117  /// types since we only care about the edge list. We
118  /// cannot handle missing edge numbers as we could with
119  /// missing element numbers due to the text block format.
120  std::string edgeStr;
121  int i,indx;
122  int nextEdgeNumber = -1;
123 
124  // Curved Edges
125  map<int, int> edge_curved;
126  for(i = 0; i < m_curvedEdges.size(); ++i)
127  {
128  edge_curved[m_curvedEdges[i]->m_curveID] = i;
129  }
130 
131  while(edge)
132  {
133  nextEdgeNumber++;
134 
135  int err = edge->QueryIntAttribute("ID",&indx);
136  ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
137 // ASSERTL0(indx == nextEdgeNumber, "Edge IDs must begin with zero and be sequential.");
138 
139  TiXmlNode *child = edge->FirstChild();
140  edgeStr.clear();
141  if (child->Type() == TiXmlNode::TEXT)
142  {
143  edgeStr += child->ToText()->ValueStr();
144  }
145 
146  /// Now parse out the edges, three fields at a time.
147  int vertex1, vertex2;
148  std::istringstream edgeDataStrm(edgeStr.c_str());
149 
150  try
151  {
152  while (!edgeDataStrm.fail())
153  {
154  edgeDataStrm >> vertex1 >> vertex2;
155 
156  // Must check after the read because we may be
157  // at the end and not know it. If we are at
158  // the end we will add a duplicate of the last
159  // entry if we don't check here.
160  if (!edgeDataStrm.fail())
161  {
162  PointGeomSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
163 
164  SegGeomSharedPtr edge;
165 
166  if(edge_curved.count(indx) == 0)
167  {
168  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices);
169  edge->SetGlobalID(indx); // Set global mesh id
170  }
171  else
172  {
173  edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_spaceDimension, vertices, m_curvedEdges[edge_curved.find(indx)->second]);
174  edge->SetGlobalID(indx); //Set global mesh id
175  }
176 
177  m_segGeoms[indx] = edge;
178  }
179  }
180  }
181  catch(...)
182  {
183  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
184  }
185 
186  edge = edge->NextSiblingElement("E");
187  }
188 
189  }
190 
191  void MeshGraph2D::ReadElements(TiXmlDocument &doc)
192  {
193  /// We know we have it since we made it this far.
194  TiXmlHandle docHandle(&doc);
195  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
196  TiXmlElement* field = NULL;
197 
198  /// Look for elements in ELEMENT block.
199  field = mesh->FirstChildElement("ELEMENT");
200 
201  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
202 
203  // Set up curve map for curved elements on an embedded manifold.
204  map<int, int> faceCurves;
205  map<int,int>::iterator x;
206  for (int i = 0; i < m_curvedFaces.size(); ++i)
207  {
208  faceCurves[m_curvedFaces[i]->m_curveID] = i;
209  }
210 
211  int nextElementNumber = -1;
212 
213  /// All elements are of the form: "<? ID="#"> ... </?>", with
214  /// ? being the element type.
215 
216  TiXmlElement *element = field->FirstChildElement();
217 
218  while (element)
219  {
220  std::string elementType(element->ValueStr());
221 
222  ASSERTL0(elementType == "Q" || elementType == "T",
223  (std::string("Unknown 2D element type: ") + elementType).c_str());
224 
225  /// These should be ordered.
226  nextElementNumber++;
227 
228  /// Read id attribute.
229  int indx;
230  int err = element->QueryIntAttribute("ID", &indx);
231  ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
232 // ASSERTL0(indx == nextElementNumber, "Element IDs must begin with zero and be sequential.");
233 
234  /// Read text element description.
235  TiXmlNode* elementChild = element->FirstChild();
236  std::string elementStr;
237  while(elementChild)
238  {
239  if (elementChild->Type() == TiXmlNode::TEXT)
240  {
241  elementStr += elementChild->ToText()->ValueStr();
242  }
243  elementChild = elementChild->NextSibling();
244  }
245 
246  ASSERTL0(!elementStr.empty(), "Unable to read element description body.");
247 
248  /// Parse out the element components corresponding to type of element.
249  if (elementType == "T")
250  {
251  // Read three edge numbers
252  int edge1, edge2, edge3;
253  std::istringstream elementDataStrm(elementStr.c_str());
254 
255  try
256  {
257  elementDataStrm >> edge1;
258  elementDataStrm >> edge2;
259  elementDataStrm >> edge3;
260 
261  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
262 
263  /// Create a TriGeom to hold the new definition.
264  SegGeomSharedPtr edges[TriGeom::kNedges] =
265  {
266  GetSegGeom(edge1),
267  GetSegGeom(edge2),
268  GetSegGeom(edge3)
269  };
270 
271  StdRegions::Orientation edgeorient[TriGeom::kNedges] =
272  {
273  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
274  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
275  SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
276  };
277 
278  TriGeomSharedPtr trigeom;
279  if ((x = faceCurves.find(indx)) == faceCurves.end())
280  {
281  trigeom = MemoryManager<TriGeom>
282  ::AllocateSharedPtr(indx,
283  edges,
284  edgeorient);
285  }
286  else
287  {
288  trigeom = MemoryManager<TriGeom>
289  ::AllocateSharedPtr(indx,
290  edges,
291  edgeorient,
292  m_curvedFaces[x->second]);
293  }
294  trigeom->SetGlobalID(indx);
295 
296  m_triGeoms[indx] = trigeom;
297  }
298  catch(...)
299  {
300  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
301  }
302  }
303  else if (elementType == "Q")
304  {
305  // Read four edge numbers
306  int edge1, edge2, edge3, edge4;
307  std::istringstream elementDataStrm(elementStr.c_str());
308 
309  try
310  {
311  elementDataStrm >> edge1;
312  elementDataStrm >> edge2;
313  elementDataStrm >> edge3;
314  elementDataStrm >> edge4;
315 
316  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
317 
318  /// Create a QuadGeom to hold the new definition.
319  SegGeomSharedPtr edges[QuadGeom::kNedges] =
320  {GetSegGeom(edge1),GetSegGeom(edge2),
321  GetSegGeom(edge3),GetSegGeom(edge4)};
322 
323  StdRegions::Orientation edgeorient[QuadGeom::kNedges] =
324  {
325  SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
326  SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
327  SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
328  SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
329  };
330 
331  QuadGeomSharedPtr quadgeom;
332  if ((x = faceCurves.find(indx)) == faceCurves.end())
333  {
334  quadgeom = MemoryManager<QuadGeom>
335  ::AllocateSharedPtr(indx,
336  edges,
337  edgeorient);
338  }
339  else
340  {
341  quadgeom = MemoryManager<QuadGeom>
342  ::AllocateSharedPtr(indx,
343  edges,
344  edgeorient,
345  m_curvedFaces[x->second]);
346  }
347  quadgeom->SetGlobalID(indx);
348 
349  m_quadGeoms[indx] = quadgeom;
350 
351  }
352  catch(...)
353  {
354  NEKERROR(ErrorUtil::efatal,(std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
355  }
356  }
357 
358  /// Keep looking
359  element = element->NextSiblingElement();
360  }
361  }
362 
363  void MeshGraph2D::ReadComposites(TiXmlDocument &doc)
364  {
365  TiXmlHandle docHandle(&doc);
366 
367  /// We know we have it since we made it this far.
368  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
369  TiXmlElement* field = NULL;
370 
371  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
372 
373  /// Look for elements in ELEMENT block.
374  field = mesh->FirstChildElement("COMPOSITE");
375 
376  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
377 
378  int nextCompositeNumber = -1;
379 
380  /// All elements are of the form: "<C ID = "N"> ... </C>".
381  /// Read the ID field first.
382  TiXmlElement *composite = field->FirstChildElement("C");
383 
384  while (composite)
385  {
386  nextCompositeNumber++;
387 
388  int indx;
389  int err = composite->QueryIntAttribute("ID", &indx);
390  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
391 // ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with zero and be sequential.");
392 
393  TiXmlNode* compositeChild = composite->FirstChild();
394  // This is primarily to skip comments that may be present.
395  // Comments appear as nodes just like elements.
396  // We are specifically looking for text in the body
397  // of the definition.
398  while(compositeChild && compositeChild->Type() != TiXmlNode::TEXT)
399  {
400  compositeChild = compositeChild->NextSibling();
401  }
402 
403  ASSERTL0(compositeChild, "Unable to read composite definition body.");
404  std::string compositeStr = compositeChild->ToText()->ValueStr();
405 
406  /// Parse out the element components corresponding to type of element.
407 
408  std::istringstream compositeDataStrm(compositeStr.c_str());
409 
410  try
411  {
412  bool first = true;
413  std::string prevCompositeElementStr;
414 
415  while (!compositeDataStrm.fail())
416  {
417  std::string compositeElementStr;
418  compositeDataStrm >> compositeElementStr;
419 
420  if (!compositeDataStrm.fail())
421  {
422  if (first)
423  {
424  first = false;
425 
426  Composite curVector(MemoryManager<GeometryVector>::AllocateSharedPtr());
427  m_meshComposites[indx] = curVector;
428  }
429 
430  if (compositeElementStr.length() > 0)
431  {
432  ResolveGeomRef(prevCompositeElementStr, compositeElementStr, m_meshComposites[indx]);
433  }
434  prevCompositeElementStr = compositeElementStr;
435  }
436  }
437  }
438  catch(...)
439  {
440  NEKERROR(ErrorUtil::efatal,
441  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
442  }
443 
444  /// Keep looking
445  composite = composite->NextSiblingElement("C");
446  }
447  }
448 
449 
450  SegGeomSharedPtr MeshGraph2D::GetSegGeom(int eID)
451  {
452  SegGeomSharedPtr returnval;
453  SegGeomMap::iterator x = m_segGeoms.find(eID);
454  ASSERTL0(x != m_segGeoms.end(), "Segment not found.");
455  return x->second;
456  };
457 
458 
459  // Take the string that is the composite reference and find the
460  // pointer to the Geometry object corresponding to it.
461 
462  // The only allowable combinations of previous and current items
463  // are V (0D); E (1D); and T and Q (2D). Only elements of the same
464  // dimension are allowed to be grouped.
465  void MeshGraph2D::ResolveGeomRef(const std::string &prevToken, const std::string &token,
466  Composite& composite)
467  {
468  try
469  {
470  std::istringstream tokenStream(token);
471  std::istringstream prevTokenStream(prevToken);
472 
473  char type;
474  char prevType;
475 
476  tokenStream >> type;
477 
478  std::string::size_type indxBeg = token.find_first_of('[') + 1;
479  std::string::size_type indxEnd = token.find_last_of(']') - 1;
480 
481  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());
482 
483  std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
484  std::vector<unsigned int> seqVector;
485  std::vector<unsigned int>::iterator seqIter;
486 
487  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
488 
489  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
490 
491  prevTokenStream >> prevType;
492 
493  // All composites must be of the same dimension.
494  bool validSequence = (prevToken.empty() || // No previous, then current is just fine.
495  (type == 'V' && prevType == 'V') ||
496  (type == 'E' && prevType == 'E') ||
497  ((type == 'T' || type == 'Q') &&
498  (prevType == 'T' || prevType == 'Q')));
499 
500  ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
501  + type + " and " + prevType + ".").c_str());
502 
503 
504  switch(type)
505  {
506  case 'E': // Edge
507  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
508  {
509  if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
510  {
511  char errStr[16] = "";
512  ::sprintf(errStr, "%d", *seqIter);
513  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
514  }
515  else
516  {
517  composite->push_back(m_segGeoms[*seqIter]);
518  }
519  }
520  break;
521 
522  case 'T': // Triangle
523  {
524  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
525  {
526  if (m_triGeoms.count(*seqIter) == 0 )
527  {
528  char errStr[16] = "";
529  ::sprintf(errStr, "%d", *seqIter);
530  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
531  }
532  else
533  {
534  if(CheckRange(*m_triGeoms[*seqIter]))
535  {
536  composite->push_back(m_triGeoms[*seqIter]);
537  }
538  }
539  }
540  }
541  break;
542 
543  case 'Q': // Quad
544  {
545  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
546  {
547  if (m_quadGeoms.count(*seqIter) == 0)
548  {
549  char errStr[16] = "";
550  ::sprintf(errStr, "%d", *seqIter);
551  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr +std::string(" in Composite section")).c_str());
552  }
553  else
554  {
555  if(CheckRange(*m_quadGeoms[*seqIter]))
556  {
557  composite->push_back(m_quadGeoms[*seqIter]);
558  }
559  }
560  }
561  }
562  break;
563 
564  case 'V': // Vertex
565  for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
566  {
567  if (*seqIter >= m_vertSet.size())
568  {
569  char errStr[16] = "";
570  ::sprintf(errStr, "%d", *seqIter);
571  NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
572  }
573  else
574  {
575  composite->push_back(m_vertSet[*seqIter]);
576  }
577  }
578  break;
579 
580  default:
581  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
582  }
583  }
584  catch(...)
585  {
586  NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
587  }
588 
589  return;
590  }
591 
592  ElementEdgeVectorSharedPtr MeshGraph2D::GetElementsFromEdge(Geometry1DSharedPtr edge)
593  {
594  SegGeomSharedPtr Sedge;
595 
596  if(!(Sedge = boost::dynamic_pointer_cast<SegGeom>(edge)))
597  {
598  ASSERTL0(false,"Dynamics cast failed");
599 
600  }
601  return GetElementsFromEdge(Sedge);
602 
603  }
604  ElementEdgeVectorSharedPtr MeshGraph2D::GetElementsFromEdge(SegGeomSharedPtr edge)
605  {
606  // Search tris and quads
607  // Need to iterate through vectors because there may be multiple
608  // occurrences.
609  ElementEdgeSharedPtr elementEdge;
610  //TriGeomVectorIter triIter;
611 
612  ElementEdgeVectorSharedPtr returnval = MemoryManager<ElementEdgeVector>::AllocateSharedPtr();
613 
614  CompositeMapIter compIter;
615  TriGeomSharedPtr triGeomShPtr;
616  QuadGeomSharedPtr quadGeomShPtr;
617 
618  GeometryVectorIter geomIter;
619 
620  for(int d = 0; d < m_domain.size(); ++d)
621  {
622  for (compIter = m_domain[d].begin(); compIter != m_domain[d].end(); ++compIter)
623  {
624  for (geomIter = (compIter->second)->begin(); geomIter != (compIter->second)->end(); ++geomIter)
625  {
626  triGeomShPtr = boost::dynamic_pointer_cast<TriGeom>(*geomIter);
627  quadGeomShPtr = boost::dynamic_pointer_cast<QuadGeom>(*geomIter);
628 
629  if (triGeomShPtr || quadGeomShPtr)
630  {
631  int edgeNum;
632  if (triGeomShPtr)
633  {
634  if ((edgeNum = triGeomShPtr->WhichEdge(edge)) > -1)
635  {
636  elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
637  elementEdge->m_Element = triGeomShPtr;
638  elementEdge->m_EdgeIndx = edgeNum;
639  returnval->push_back(elementEdge);
640  }
641  }
642  else if (quadGeomShPtr)
643  {
644  if ((edgeNum = quadGeomShPtr->WhichEdge(edge)) > -1)
645  {
646  elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
647  elementEdge->m_Element = quadGeomShPtr;
648  elementEdge->m_EdgeIndx = edgeNum;
649  returnval->push_back(elementEdge);
650  }
651  }
652  }
653  }
654  }
655  }
656 
657  return returnval;
658  }
659 
660  LibUtilities::BasisKey MeshGraph2D::GetEdgeBasisKey(SegGeomSharedPtr edge, const std::string variable)
661  {
662  ElementEdgeVectorSharedPtr elements = GetElementsFromEdge(edge);
663  // Perhaps, a check should be done here to ensure that
664  // in case elements->size!=1, all elements to which
665  // the edge belongs have the same type and order of
666  // expansion such that no confusion can arise.
667  ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element, variable);
668 
669  int edge_id = (*elements)[0]->m_EdgeIndx;
670 
671  if((*elements)[0]->m_Element->GetShapeType() == LibUtilities::eTriangle)
672  {
673  edge_id = (edge_id)? 1:0;
674  }
675  else
676  {
677  edge_id = edge_id%2;
678  }
679 
680  int nummodes = expansion->m_basisKeyVector[edge_id].GetNumModes();
681  int numpoints = expansion->m_basisKeyVector[edge_id].GetNumPoints();
682 
683  if((*elements)[0]->m_Element->GetShapeType() == LibUtilities::eTriangle)
684  {
685  // Use edge 0 to define basis of order relevant to edge
686  switch(expansion->m_basisKeyVector[edge_id].GetBasisType())
687  {
688  case LibUtilities::eGLL_Lagrange:
689  {
690  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
691  {
692  case LibUtilities::eGaussLobattoLegendre:
693  {
694  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
695  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
696  }
697  break;
698 
699  default:
700  ASSERTL0(false,"Unexpected points distribution");
701 
702  // It doesn't matter what we return
703  // here since the ASSERT will stop
704  // execution. Just return something
705  // to prevent warnings messages.
706  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
707  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
708  break;
709  }
710  }
711  break;
712  case LibUtilities::eOrtho_B: // Assume this is called from nodal triangular basis
713  {
714  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
715  {
716  case LibUtilities::eGaussRadauMAlpha1Beta0:
717  {
718  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
719  return LibUtilities::BasisKey(LibUtilities::eGLL_Lagrange,nummodes,pkey);
720  }
721  break;
722 
723  default:
724  ASSERTL0(false,"Unexpected points distribution");
725 
726  // It doesn't matter what we return
727  // here since the ASSERT will stop
728  // execution. Just return something
729  // to prevent warnings messages.
730  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
731  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
732  break;
733  }
734  }
735  break;
736  case LibUtilities::eModified_B:
737  {
738  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
739  {
740  case LibUtilities::eGaussRadauMAlpha1Beta0:
741  {
742  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
743  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
744  }
745  break;
746 
747  default:
748  ASSERTL0(false,"Unexpected points distribution");
749 
750  // It doesn't matter what we return
751  // here since the ASSERT will stop
752  // execution. Just return something
753  // to prevent warnings messages.
754  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
755  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
756  break;
757  }
758  }
759  break;
760  case LibUtilities::eModified_A:
761  {
762  switch(expansion->m_basisKeyVector[edge_id].GetPointsType())
763  {
764  case LibUtilities::eGaussLobattoLegendre:
765  {
766  const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussLobattoLegendre);
767  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
768  }
769  break;
770  default:
771  ASSERTL0(false,"Unexpected points distribution");
772  // It doesn't matter what we return here
773  // since the ASSERT will stop execution.
774  // Just return something to prevent
775  // warnings messages.
776  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
777  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
778  break;
779  }
780  }
781  break;
782 
783  default:
784  ASSERTL0(false,"Unexpected basis distribution");
785  // It doesn't matter what we return here since the
786  // ASSERT will stop execution. Just return
787  // something to prevent warnings messages.
788  const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
789  return LibUtilities::BasisKey(expansion->m_basisKeyVector[0].GetBasisType(),nummodes,pkey);
790  }
791  }
792  else
793  {
794  // Quadrilateral
795  const LibUtilities::PointsKey pkey(numpoints,expansion->m_basisKeyVector[edge_id].GetPointsType());
796  return LibUtilities::BasisKey(expansion->m_basisKeyVector[edge_id].GetBasisType(),nummodes,pkey);
797  }
798 
799  ASSERTL0(false, "Unable to determine edge points type.");
800  return LibUtilities::NullBasisKey;
801  }
802  }; //end of namespace
803 }; //end of namespace
804 
805 //
806 // $Log: MeshGraph2D.cpp,v $
807 // Revision 1.41 2009/12/17 02:08:04 bnelson
808 // Fixed visual studio compiler warning.
809 //
810 // Revision 1.40 2009/10/22 17:41:47 cbiotto
811 // Update for variable order expansion
812 //
813 // Revision 1.39 2009/07/02 13:32:24 sehunchun
814 // *** empty log message ***
815 //
816 // Revision 1.38 2009/07/02 11:39:44 sehunchun
817 // Modification for 2D gemoetry embedded in 3D
818 //
819 // Revision 1.37 2009/05/01 13:23:21 pvos
820 // Fixed various bugs
821 //
822 // Revision 1.36 2009/04/20 16:13:23 sherwin
823 // Modified Import and Write functions and redefined how Expansion is used
824 //
825 // Revision 1.35 2009/01/21 16:59:03 pvos
826 // Added additional geometric factors to improve efficiency
827 //
828 // Revision 1.34 2009/01/12 10:26:59 pvos
829 // Added input tags for nodal expansions
830 //
831 // Revision 1.33 2008/09/12 11:26:19 pvos
832 // Updates for mappings in 3D
833 //
834 // Revision 1.32 2008/09/09 14:21:44 sherwin
835 // Updates for first working version of curved edges
836 //
837 // Revision 1.31 2008/08/14 22:11:03 sherwin
838 // Mods for HDG update
839 //
840 // Revision 1.30 2008/07/29 22:23:36 sherwin
841 // various mods for DG advection solver in Multiregions. Added virtual calls to Geometry, Geometry1D, 2D and 3D
842 //
843 // Revision 1.29 2008/05/29 19:07:39 delisi
844 // 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.
845 //
846 // Revision 1.28 2008/05/28 21:42:18 jfrazier
847 // Minor comment spelling change.
848 //
849 // Revision 1.27 2008/03/18 14:14:49 pvos
850 // Update for nodal triangular helmholtz solver
851 //
852 // Revision 1.26 2008/03/11 17:02:24 jfrazier
853 // Modified the GetElementsFromEdge to use the domain.
854 //
855 // Revision 1.25 2008/01/21 19:58:14 sherwin
856 // Updated so that QuadGeom and TriGeom have SegGeoms instead of EdgeComponents
857 //
858 // Revision 1.24 2007/12/17 20:27:24 sherwin
859 // Added normals to GeomFactors
860 //
861 // Revision 1.23 2007/12/14 17:39:18 jfrazier
862 // Fixed composite reader to handle ranges and comma separated lists.
863 //
864 // Revision 1.22 2007/12/11 21:51:53 jfrazier
865 // Updated 2d components so elements could be retrieved from edges.
866 //
867 // Revision 1.21 2007/12/04 03:02:26 jfrazier
868 // Changed to stringstream.
869 //
870 // Revision 1.20 2007/09/20 22:25:06 jfrazier
871 // Added expansion information to meshgraph class.
872 //
873 // Revision 1.19 2007/07/26 01:38:33 jfrazier
874 // Cleanup of some attribute reading code.
875 //
876 // Revision 1.18 2007/07/24 16:52:09 jfrazier
877 // Added domain code.
878 //
879 // Revision 1.17 2007/07/05 04:21:10 jfrazier
880 // Changed id format and propagated from 1d to 2d.
881 //
882 // Revision 1.16 2007/06/10 02:27:11 jfrazier
883 // Another checkin with an incremental completion of the boundary conditions reader
884 //
885 // Revision 1.15 2007/06/07 23:55:24 jfrazier
886 // Intermediate revisions to add parsing for boundary conditions file.
887 //
888 // Revision 1.14 2007/05/28 21:48:42 sherwin
889 // Update for 2D functionality
890 //
891 // Revision 1.13 2006/10/17 22:26:01 jfrazier
892 // Added capability to specify ranges in composite definition.
893 //
894 // Revision 1.12 2006/10/17 18:42:54 jfrazier
895 // Removed "NUMBER" attribute in items.
896 //
897 // Revision 1.11 2006/09/26 23:41:53 jfrazier
898 // Updated to account for highest level NEKTAR tag and changed the geometry tag to GEOMETRY.
899 //
900 // Revision 1.10 2006/08/24 18:50:00 jfrazier
901 // Completed error checking on permissable composite item combinations.
902 //
903 // Revision 1.9 2006/08/18 19:37:17 jfrazier
904 // *** empty log message ***
905 //
906 // Revision 1.8 2006/08/17 22:55:00 jfrazier
907 // Continued adding code to process composites in the mesh2d.
908 //
909 // Revision 1.7 2006/08/16 23:34:42 jfrazier
910 // *** empty log message ***
911 //
912 // Revision 1.6 2006/06/02 18:48:40 sherwin
913 // Modifications to make ProjectLoc2D run bit there are bus errors for order > 3
914 //
915 // Revision 1.5 2006/06/01 14:15:30 sherwin
916 // Added typdef of boost wrappers and made GeoFac a boost shared pointer.
917 //
918 // Revision 1.4 2006/05/30 14:00:04 sherwin
919 // Updates to make MultiRegions and its Demos work
920 //
921 // Revision 1.3 2006/05/23 19:56:33 jfrazier
922 // These build and run, but the expansion pieces are commented out
923 // because they would not run.
924 //
925 // Revision 1.2 2006/05/09 13:37:01 jfrazier
926 // Removed duplicate definition of shared vertex pointer.
927 //
928 // Revision 1.1 2006/05/04 18:59:01 kirby
929 // *** empty log message ***
930 //
931 // Revision 1.11 2006/04/11 23:18:11 jfrazier
932 // Completed MeshGraph2D for tri's and quads. Not thoroughly tested.
933 //
934 // Revision 1.10 2006/04/09 02:08:35 jfrazier
935 // Added precompiled header.
936 //
937 // Revision 1.9 2006/04/04 23:12:37 jfrazier
938 // More updates to readers. Still work to do on MeshGraph2D to store tris and quads.
939 //
940 // Revision 1.8 2006/03/25 00:58:29 jfrazier
941 // Many changes dealing with fundamental structure and reading/writing.
942 //
943 // Revision 1.7 2006/03/12 14:20:43 sherwin
944 //
945 // First compiling version of SpatialDomains and associated modifications
946 //
947 // Revision 1.6 2006/03/12 07:42:03 sherwin
948 //
949 // Updated member names and StdRegions call. Still has not been compiled
950 //
951 // Revision 1.5 2006/02/26 21:19:43 bnelson
952 // Fixed a variety of compiler errors caused by updates to the coding standard.
953 //
954 // Revision 1.4 2006/02/19 01:37:34 jfrazier
955 // Initial attempt at bringing into conformance with the coding standard. Still more work to be done. Has not been compiled.
956 //
957 //
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2