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