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