Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PyrGeom.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: PyrGeom.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: Pyramidic geometry information.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <SpatialDomains/PyrGeom.h>
37 
40 #include <StdRegions/StdPyrExp.h>
41 #include <SpatialDomains/SegGeom.h>
42 #include <SpatialDomains/TriGeom.h>
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49  namespace SpatialDomains
50  {
51  PyrGeom::PyrGeom()
52  {
53  m_shapeType = LibUtilities::ePyramid;
54  }
55 
56  PyrGeom::PyrGeom(const Geometry2DSharedPtr faces[]) :
57  Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
58  {
60 
61  /// Copy the face shared pointers
62  m_faces.insert(m_faces.begin(), faces, faces+PyrGeom::kNfaces);
63 
64  /// Set up orientation vectors with correct amount of elements.
65  m_eorient.resize(kNedges);
66  m_forient.resize(kNfaces);
67 
72  SetUpXmap();
73  SetUpCoeffs(m_xmap->GetNcoeffs());
74  }
75 
77  {
78 
79  }
80 
82  {
84  {
85  int i;
86  GeomType Gtype = eRegular;
87 
88  v_FillGeom();
89 
90  // check to see if expansions are linear
91  for(i = 0; i < m_coordim; ++i)
92  {
93  if (m_xmap->GetBasisNumModes(0) != 2 ||
94  m_xmap->GetBasisNumModes(1) != 2 ||
95  m_xmap->GetBasisNumModes(2) != 2 )
96  {
97  Gtype = eDeformed;
98  }
99  }
100 
101  // check to see if all quadrilateral faces are parallelograms
102  if(Gtype == eRegular)
103  {
104  // Ensure each face is a parallelogram? Check this.
105  for (i = 0; i < m_coordim; i++)
106  {
107  if( fabs( (*m_verts[0])(i) -
108  (*m_verts[1])(i) +
109  (*m_verts[2])(i) -
110  (*m_verts[3])(i) )
112  {
113  Gtype = eDeformed;
114  break;
115  }
116  }
117  }
118 
120  Gtype, m_coordim, m_xmap, m_coeffs);
122  }
123  }
124 
126  const Array<OneD, const NekDouble> &coords,
127  Array<OneD, NekDouble> &Lcoords)
128  {
129  NekDouble ptdist = 1e6;
130 
131  v_FillGeom();
132 
133  // calculate local coordinate for coord
134  if(GetMetricInfo()->GetGtype() == eRegular)
135  { // Based on Spen's book, page 99
136 
137  // Point inside tetrahedron
138  PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
139 
140  // Edges
141  PointGeom er0, e10, e30, e40;
142  er0.Sub(r,*m_verts[0]);
143  e10.Sub(*m_verts[1],*m_verts[0]);
144  e30.Sub(*m_verts[3],*m_verts[0]);
145  e40.Sub(*m_verts[4],*m_verts[0]);
146 
147 
148  // Cross products (Normal times area)
149  PointGeom cp1030, cp3040, cp4010;
150  cp1030.Mult(e10,e30);
151  cp3040.Mult(e30,e40);
152  cp4010.Mult(e40,e10);
153 
154 
155  // Barycentric coordinates (relative volume)
156  NekDouble V = e40.dot(cp1030); // Pyramid Volume = {(e40)dot(e10)x(e30)}/4
157  NekDouble scaleFactor = 2.0/3.0;
158  NekDouble v1 = er0.dot(cp3040) / V; // volume1 = {(er0)dot(e30)x(e40)}/6
159  NekDouble v2 = er0.dot(cp4010) / V; // volume2 = {(er0)dot(e40)x(e10)}/6
160  NekDouble beta = v1 * scaleFactor;
161  NekDouble gamma = v2 * scaleFactor;
162  NekDouble delta = er0.dot(cp1030) / V; // volume3 = {(er0)dot(e10)x(e30)}/4
163 
164  // Make Pyramid bigger
165  Lcoords[0] = 2.0*beta - 1.0;
166  Lcoords[1] = 2.0*gamma - 1.0;
167  Lcoords[2] = 2.0*delta - 1.0;
168 
169  // Set ptdist to distance to nearest vertex
170  for(int i = 0; i < 5; ++i)
171  {
172  ptdist = min(ptdist,r.dist(*m_verts[i]));
173  }
174 
175  }
176  else
177  {
179  "inverse mapping must be set up to use this call");
180  }
181  return ptdist;
182  }
183 
185  {
186  return 5;
187  }
188 
190  {
191  return 8;
192  }
193 
195  {
196  return 5;
197  }
198 
199  int PyrGeom::v_GetDir(const int faceidx, const int facedir) const
200  {
201  if (faceidx == 0)
202  {
203  return facedir;
204  }
205  else if (faceidx == 1 || faceidx == 3)
206  {
207  return 2 * facedir;
208  }
209  else
210  {
211  return 1 + facedir;
212  }
213  }
214 
216  {
217  // find edge 0
218  int i,j;
219  unsigned int check;
220 
221  SegGeomSharedPtr edge;
222 
223  // First set up the 4 bottom edges
224  int f;
225  for (f = 1; f < 5; f++)
226  {
227  int nEdges = m_faces[f]->GetNumEdges();
228  check = 0;
229  for (i = 0; i < 4; i++)
230  {
231  for (j = 0; j < nEdges; j++)
232  {
233  if (m_faces[0]->GetEid(i) == m_faces[f]->GetEid(j))
234  {
235  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[0])->GetEdge(i));
236  m_edges.push_back(edge);
237  check++;
238  }
239  }
240  }
241 
242  if (check < 1)
243  {
244  std::ostringstream errstrm;
245  errstrm << "Connected faces do not share an edge. Faces ";
246  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
247  ASSERTL0(false, errstrm.str());
248  }
249  else if (check > 1)
250  {
251  std::ostringstream errstrm;
252  errstrm << "Connected faces share more than one edge. Faces ";
253  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
254  ASSERTL0(false, errstrm.str());
255  }
256  }
257 
258  // Then, set up the 4 vertical edges
259  check = 0;
260  for(i = 0; i < 3; i++) //Set up the vertical edge :face(1) and face(4)
261  {
262  for(j = 0; j < 3; j++)
263  {
264  if( (m_faces[1])->GetEid(i) == (m_faces[4])->GetEid(j) )
265  {
266  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[1])->GetEdge(i));
267  m_edges.push_back(edge);
268  check++;
269  }
270  }
271  }
272  if( check < 1 )
273  {
274  std::ostringstream errstrm;
275  errstrm << "Connected faces do not share an edge. Faces ";
276  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
277  ASSERTL0(false, errstrm.str());
278  }
279  else if( check > 1)
280  {
281  std::ostringstream errstrm;
282  errstrm << "Connected faces share more than one edge. Faces ";
283  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
284  ASSERTL0(false, errstrm.str());
285  }
286 
287  // Set up vertical edges: face(1) through face(4)
288  for (f = 1; f < 4; f++)
289  {
290  check = 0;
291  for(i = 0; i < m_faces[f]->GetNumEdges(); i++)
292  {
293  for(j = 0; j < m_faces[f+1]->GetNumEdges(); j++)
294  {
295  if( (m_faces[f])->GetEid(i) == (m_faces[f+1])->GetEid(j))
296  {
297  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[f])->GetEdge(i));
298  m_edges.push_back(edge);
299  check++;
300  }
301  }
302  }
303 
304  if( check < 1 )
305  {
306  std::ostringstream errstrm;
307  errstrm << "Connected faces do not share an edge. Faces ";
308  errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
309  ASSERTL0(false, errstrm.str());
310  }
311  else if( check > 1)
312  {
313  std::ostringstream errstrm;
314  errstrm << "Connected faces share more than one edge. Faces ";
315  errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
316  ASSERTL0(false, errstrm.str());
317  }
318  }
319  }
320 
322  {
323  // Set up the first 2 vertices (i.e. vertex 0,1)
324  if (m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0) ||
325  m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1))
326  {
327  m_verts.push_back(m_edges[0]->GetVertex(1));
328  m_verts.push_back(m_edges[0]->GetVertex(0));
329  }
330  else if (m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0) ||
331  m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1))
332  {
333  m_verts.push_back(m_edges[0]->GetVertex(0));
334  m_verts.push_back(m_edges[0]->GetVertex(1));
335  }
336  else
337  {
338  std::ostringstream errstrm;
339  errstrm << "Connected edges do not share a vertex. Edges ";
340  errstrm << m_edges[0]->GetEid() << ", " << m_edges[1]->GetEid();
341  ASSERTL0(false, errstrm.str());
342  }
343 
344  // set up the other bottom vertices (i.e. vertex 2,3)
345  for(int i = 1; i < 3; i++)
346  {
347  if (m_edges[i]->GetVid(0) == m_verts[i]->GetVid())
348  {
349  m_verts.push_back(m_edges[i]->GetVertex(1));
350  }
351  else if (m_edges[i]->GetVid(1) == m_verts[i]->GetVid())
352  {
353  m_verts.push_back(m_edges[i]->GetVertex(0));
354  }
355  else
356  {
357  std::ostringstream errstrm;
358  errstrm << "Connected edges do not share a vertex. Edges ";
359  errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
360  ASSERTL0(false, errstrm.str());
361  }
362  }
363 
364  // set up top vertex
365  if (m_edges[4]->GetVid(0) == m_verts[0]->GetVid())
366  {
367  m_verts.push_back(m_edges[4]->GetVertex(1));
368  }
369  else
370  {
371  m_verts.push_back(m_edges[4]->GetVertex(0));
372  }
373 
374  int check = 0;
375  for (int i = 5; i < 8; ++i)
376  {
377  if( (m_edges[i]->GetVid(0) == m_verts[i-4]->GetVid()
378  && m_edges[i]->GetVid(1) == m_verts[4]->GetVid())
379  ||(m_edges[i]->GetVid(1) == m_verts[i-4]->GetVid()
380  && m_edges[i]->GetVid(0) == m_verts[4]->GetVid()))
381  {
382  check++;
383  }
384  }
385  if (check != 3) {
386  std::ostringstream errstrm;
387  errstrm << "Connected edges do not share a vertex. Edges ";
388  errstrm << m_edges[3]->GetEid() << ", " << m_edges[2]->GetEid();
389  ASSERTL0(false, errstrm.str());
390  }
391  }
392 
394  {
395  // This 2D array holds the local id's of all the vertices for every
396  // edge. For every edge, they are ordered to what we define as being
397  // Forwards.
398  const unsigned int edgeVerts[kNedges][2] =
399  { {0,1}, {1,2}, {3,2}, {0,3}, {0,4}, {1,4}, {2,4}, {3,4} };
400 
401  int i;
402  for (i = 0; i < kNedges; i++)
403  {
404  if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][0]]->GetVid())
405  {
407  }
408  else if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][1]]->GetVid())
409  {
411  }
412  else
413  {
414  ASSERTL0(false,"Could not find matching vertex for the edge");
415  }
416  }
417  }
418 
420  {
421  int f,i;
422 
423  // These arrays represent the vector of the A and B coordinate of
424  // the local elemental coordinate system where A corresponds with
425  // the coordinate direction xi_i with the lowest index i (for that
426  // particular face) Coordinate 'B' then corresponds to the other
427  // local coordinate (i.e. with the highest index)
428  Array<OneD,NekDouble> elementAaxis(m_coordim);
429  Array<OneD,NekDouble> elementBaxis(m_coordim);
430 
431  // These arrays correspond to the local coordinate
432  // system of the face itself (i.e. the Geometry2D)
433  // faceAaxis correspond to the xi_0 axis
434  // faceBaxis correspond to the xi_1 axis
435  Array<OneD,NekDouble> faceAaxis(m_coordim);
436  Array<OneD,NekDouble> faceBaxis(m_coordim);
437 
438  // This is the base vertex of the face (i.e. the Geometry2D) This
439  // corresponds to thevertex with local ID 0 of the Geometry2D
440  unsigned int baseVertex;
441 
442  // The lenght of the vectors above
443  NekDouble elementAaxis_length;
444  NekDouble elementBaxis_length;
445  NekDouble faceAaxis_length;
446  NekDouble faceBaxis_length;
447 
448  // This 2D array holds the local id's of all the vertices for every
449  // face. For every face, they are ordered in such a way that the
450  // implementation below allows a unified approach for all faces.
451  const unsigned int faceVerts[kNfaces][4] = {
452  {0,1,2,3},
453  {0,1,4,0}, // Last four elements are triangles which only
454  {1,2,4,0}, // require three vertices.
455  {3,2,4,0},
456  {0,3,4,0}
457  };
458 
459  NekDouble dotproduct1 = 0.0;
460  NekDouble dotproduct2 = 0.0;
461 
462  unsigned int orientation;
463 
464  // Loop over all the faces to set up the orientation
465  for(f = 0; f < kNqfaces + kNtfaces; f++)
466  {
467  // initialisation
468  elementAaxis_length = 0.0;
469  elementBaxis_length = 0.0;
470  faceAaxis_length = 0.0;
471  faceBaxis_length = 0.0;
472 
473  dotproduct1 = 0.0;
474  dotproduct2 = 0.0;
475 
476  baseVertex = m_faces[f]->GetVid(0);
477 
478  // We are going to construct the vectors representing the A and
479  // B axis of every face. These vectors will be constructed as a
480  // vector-representation of the edges of the face. However, for
481  // both coordinate directions, we can represent the vectors by
482  // two different edges. That's why we need to make sure that we
483  // pick the edge to which the baseVertex of the
484  // Geometry2D-representation of the face belongs...
485 
486  // Compute the length of edges on a base-face
487  if (f > 0)
488  {
489  for(i = 0; i < m_coordim; i++)
490  {
491  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
492  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
493  }
494  }
495  else
496  {
497  if( baseVertex == m_verts[ faceVerts[f][0] ]->GetVid() )
498  {
499  for(i = 0; i < m_coordim; i++)
500  {
501  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
502  elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
503  }
504  }
505  else if( baseVertex == m_verts[ faceVerts[f][1] ]->GetVid() )
506  {
507  for(i = 0; i < m_coordim; i++)
508  {
509  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
510  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
511  }
512  }
513  else if( baseVertex == m_verts[ faceVerts[f][2] ]->GetVid() )
514  {
515  for(i = 0; i < m_coordim; i++)
516  {
517  elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
518  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
519  }
520  }
521  else if( baseVertex == m_verts[ faceVerts[f][3] ]->GetVid() )
522  {
523  for(i = 0; i < m_coordim; i++)
524  {
525  elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
526  elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
527  }
528  }
529  else
530  {
531  ASSERTL0(false, "Could not find matching vertex for the face");
532  }
533  }
534 
535  // Now, construct the edge-vectors of the local coordinates of
536  // the Geometry2D-representation of the face
537  for(i = 0; i < m_coordim; i++)
538  {
539  int v = m_faces[f]->GetNumVerts()-1;
540  faceAaxis[i] = (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
541  faceBaxis[i] = (*m_faces[f]->GetVertex(v))[i] - (*m_faces[f]->GetVertex(0))[i];
542 
543  elementAaxis_length += pow(elementAaxis[i],2);
544  elementBaxis_length += pow(elementBaxis[i],2);
545  faceAaxis_length += pow(faceAaxis[i],2);
546  faceBaxis_length += pow(faceBaxis[i],2);
547  }
548 
549  elementAaxis_length = sqrt(elementAaxis_length);
550  elementBaxis_length = sqrt(elementBaxis_length);
551  faceAaxis_length = sqrt(faceAaxis_length);
552  faceBaxis_length = sqrt(faceBaxis_length);
553 
554  // Calculate the inner product of both the A-axis
555  // (i.e. Elemental A axis and face A axis)
556  for(i = 0 ; i < m_coordim; i++)
557  {
558  dotproduct1 += elementAaxis[i]*faceAaxis[i];
559  }
560 
561  orientation = 0;
562  // if the innerproduct is equal to the (absolute value of the ) products of the lengths
563  // of both vectors, then, the coordinate systems will NOT be transposed
564  if( fabs(elementAaxis_length*faceAaxis_length - fabs(dotproduct1)) < NekConstants::kNekZeroTol )
565  {
566  // if the inner product is negative, both A-axis point
567  // in reverse direction
568  if(dotproduct1 < 0.0)
569  {
570  orientation += 2;
571  }
572 
573  // calculate the inner product of both B-axis
574  for(i = 0 ; i < m_coordim; i++)
575  {
576  dotproduct2 += elementBaxis[i]*faceBaxis[i];
577  }
578 
579  // if the inner product is negative, both B-axis point
580  // in reverse direction
581  if( dotproduct2 < 0.0 )
582  {
583  orientation++;
584  }
585  }
586  // The coordinate systems are transposed
587  else
588  {
589  orientation = 4;
590 
591  // Calculate the inner product between the elemental A-axis
592  // and the B-axis of the face (which are now the corresponding axis)
593  dotproduct1 = 0.0;
594  for(i = 0 ; i < m_coordim; i++)
595  {
596  dotproduct1 += elementAaxis[i]*faceBaxis[i];
597  }
598 
599  // check that both these axis are indeed parallel
600  if (fabs(elementAaxis_length*faceBaxis_length
601  - fabs(dotproduct1)) > NekConstants::kNekZeroTol)
602  {
603  cout << "Warning: Prism axes not parallel" << endl;
604  }
605 
606  // if the result is negative, both axis point in reverse
607  // directions
608  if(dotproduct1 < 0.0)
609  {
610  orientation += 2;
611  }
612 
613  // Do the same for the other two corresponding axis
614  dotproduct2 = 0.0;
615  for(i = 0 ; i < m_coordim; i++)
616  {
617  dotproduct2 += elementBaxis[i]*faceAaxis[i];
618  }
619 
620  // check that both these axis are indeed parallel
621  if (fabs(elementBaxis_length*faceAaxis_length
622  - fabs(dotproduct2)) > NekConstants::kNekZeroTol)
623  {
624  cout << "Warning: Prism axes not parallel" << endl;
625  }
626 
627  if( dotproduct2 < 0.0 )
628  {
629  orientation++;
630  }
631  }
632 
633  orientation = orientation + 5;
634 
635  // Fill the m_forient array
636  m_forient[f] = (StdRegions::Orientation) orientation;
637  }
638  }
639 
641  CurveMap &curvedEdges,
642  CurveMap &curvedFaces)
643  {
644  Geometry::v_Reset(curvedEdges, curvedFaces);
645 
646  for (int i = 0; i < 5; ++i)
647  {
648  m_faces[i]->Reset(curvedEdges, curvedFaces);
649  }
650 
651  SetUpXmap();
652  SetUpCoeffs(m_xmap->GetNcoeffs());
653  }
654 
655  /**
656  * @brief Set up the #m_xmap object by determining the order of each
657  * direction from derived faces.
658  */
660  {
661  vector<int> tmp;
662  int order0, order1;
663 
664  if (m_forient[0] < 9)
665  {
666  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
667  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
668  order0 = *max_element(tmp.begin(), tmp.end());
669  }
670  else
671  {
672  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
673  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
674  order0 = *max_element(tmp.begin(), tmp.end());
675  }
676 
677  if (m_forient[0] < 9)
678  {
679  tmp.clear();
680  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
681  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
682  tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
683  order1 = *max_element(tmp.begin(), tmp.end());
684  }
685  else
686  {
687  tmp.clear();
688  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
689  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
690  tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
691  order1 = *max_element(tmp.begin(), tmp.end());
692  }
693 
694  tmp.clear();
695  tmp.push_back(order0);
696  tmp.push_back(order1);
697  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(1));
698  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(2));
699  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(1));
700  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(2));
701  int order2 = *max_element(tmp.begin(), tmp.end());
702 
703 
704  const LibUtilities::BasisKey A1(
708  const LibUtilities::BasisKey A2(
712  const LibUtilities::BasisKey C(
716 
718  A1, A2, C);
719  }
720  }; //end of namespace
721 }; //end of namespace
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
#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
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:92
Principle Modified Functions .
Definition: BasisType.h:51
virtual int v_GetDir(const int faceidx, const int facedir) const
Definition: PyrGeom.cpp:199
static const int kNfaces
Definition: PyrGeom.h:60
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
GeomFactorsSharedPtr m_geomFactors
Definition: Geometry.h:170
NekDouble dot(PointGeom &a)
Definition: PointGeom.cpp:192
virtual int v_GetNumEdges() const
Definition: PyrGeom.cpp:189
int GetEid(int i) const
Return the ID of edge i in this element.
Definition: Geometry3D.cpp:74
Principle Modified Functions .
Definition: BasisType.h:49
void Sub(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:168
STL namespace.
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: Geometry.cpp:307
int GetFid(int i) const
Definition: Geometry.h:329
StdRegions::StdExpansionSharedPtr GetXmap() const
Definition: Geometry.h:383
virtual int v_GetNumVerts() const
Definition: PyrGeom.cpp:184
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: PyrGeom.cpp:640
static const NekDouble kNekZeroTol
static const int kNedges
Definition: PyrGeom.h:57
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Definition: PyrGeom.cpp:125
PointGeomSharedPtr GetVertex(int i) const
Definition: Geometry.h:348
static const int kNtfaces
Definition: PyrGeom.h:59
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:247
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
void Mult(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:177
const Geometry1DSharedPtr GetEdge(int i) const
static const int kNqfaces
Definition: PyrGeom.h:58
3D geometry information
Definition: Geometry3D.h:70
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
Geometry is straight-sided with constant geometric factors.
int GetVid(int i) const
Definition: Geometry.h:319
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
LibUtilities::ShapeType m_shapeType
Definition: Geometry.h:177
virtual int v_GetNumFaces() const
Definition: PyrGeom.cpp:194
NekDouble dist(PointGeom &a)
Definition: PointGeom.cpp:186
boost::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:63
Geometric information has been generated.
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
Definition: Geometry.h:484
GeomFactorsSharedPtr GetMetricInfo()
Definition: Geometry.h:299
GeomType
Indicates the type of element geometry.
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Definition: PyrGeom.cpp:659
virtual void v_GenGeomFactors()
Definition: PyrGeom.cpp:81
Geometry is curved or has non-constant factors.
int m_coordim
coordinate dimension
Definition: Geometry.h:169
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50