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