deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
qprojector.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/derivative_form.h>
16#include <deal.II/base/geometry_info.h>
17#include <deal.II/base/polynomials_barycentric.h>
18#include <deal.II/base/qprojector.h>
19#include <deal.II/base/tensor_product_polynomials.h>
20
21#include <deal.II/grid/reference_cell.h>
22#include <deal.II/grid/tria_orientation.h>
23
24DEAL_II_NAMESPACE_OPEN
25
26
27namespace internal
28{
29 namespace QProjector
30 {
31 namespace
32 {
33 // Reflect points across the y = x line.
34 std::vector<Point<2>>
35 reflect(const std::vector<Point<2>> &points)
36 {
37 std::vector<Point<2>> q_points;
38 q_points.reserve(points.size());
39 for (const Point<2> &p : points)
40 q_points.emplace_back(p[1], p[0]);
41
42 return q_points;
43 }
44
45
46 // Rotate points in the plane around the positive z axis (i.e.,
47 // counter-clockwise).
48 std::vector<Point<2>>
49 rotate(const std::vector<Point<2>> &points, const unsigned int n_times)
50 {
51 std::vector<Point<2>> q_points;
52 q_points.reserve(points.size());
53 switch (n_times % 4)
54 {
55 case 0:
56 // 0 degree. the point remains as it is.
57 for (const Point<2> &p : points)
58 q_points.push_back(p);
59 break;
60 case 1:
61 // 90 degree counterclockwise
62 for (const Point<2> &p : points)
63 q_points.emplace_back(1.0 - p[1], p[0]);
64 break;
65 case 2:
66 // 180 degree counterclockwise
67 for (const Point<2> &p : points)
68 q_points.emplace_back(1.0 - p[0], 1.0 - p[1]);
69 break;
70 case 3:
71 // 270 degree counterclockwise
72 for (const Point<2> &p : points)
73 q_points.emplace_back(p[1], 1.0 - p[0]);
74 break;
75 }
76
77 return q_points;
78 }
79
86 void
87 project_to_hex_face_and_append(
88 const std::vector<Point<2>> &points,
89 const unsigned int face_no,
90 std::vector<Point<3>> &q_points,
91 const RefinementCase<2> &ref_case = RefinementCase<2>::no_refinement,
92 const unsigned int subface_no = 0)
93 {
94 // one coordinate is at a const value. for faces 0, 2 and 4 this value
95 // is 0.0, for faces 1, 3 and 5 it is 1.0
96 const double const_value = face_no % 2;
97
98 // local 2d coordinates are xi and eta, global 3d coordinates are x, y
99 // and z. those have to be mapped. the following indices tell, which
100 // global coordinate (0->x, 1->y, 2->z) corresponds to which local one
101 const unsigned int xi_index = (1 + face_no / 2) % 3,
102 eta_index = (2 + face_no / 2) % 3,
103 const_index = face_no / 2;
104
105 // for a standard face (no refinement), we use the default values of
106 // the xi and eta scales and translations, otherwise the xi and eta
107 // values will be scaled (by factor 0.5 or factor 1.0) depending on
108 // the refinement case and translated (by 0.0 or 0.5) depending on the
109 // refinement case and subface_no
110 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
111 eta_translation = 0.0;
112
113 // set the scale and translation parameter for individual subfaces
114 switch (ref_case)
115 {
116 case RefinementCase<2>::no_refinement:
117 break;
118 case RefinementCase<2>::cut_x:
119 xi_scale = 0.5;
120 xi_translation = subface_no % 2 * 0.5;
121 break;
122 case RefinementCase<2>::cut_y:
123 eta_scale = 0.5;
124 eta_translation = subface_no % 2 * 0.5;
125 break;
126 case RefinementCase<2>::cut_xy:
127 xi_scale = 0.5;
128 eta_scale = 0.5;
129 xi_translation = int(subface_no % 2) * 0.5;
130 eta_translation = int(subface_no / 2) * 0.5;
131 break;
132 default:
133 DEAL_II_ASSERT_UNREACHABLE();
134 break;
135 }
136
137 // finally, compute the scaled, translated, projected quadrature
138 // points
139 for (const Point<2> &p : points)
140 {
141 Point<3> cell_point;
142 cell_point[xi_index] = xi_scale * p[0] + xi_translation;
143 cell_point[eta_index] = eta_scale * p[1] + eta_translation;
144 cell_point[const_index] = const_value;
145 q_points.push_back(cell_point);
146 }
147 }
148
149 std::vector<Point<2>>
150 mutate_points_with_offset(const std::vector<Point<2>> &points,
151 const unsigned char combined_orientation)
152 {
153 // These rotations are backwards (relative to the standard notion of,
154 // e.g., what rotation index 7 means) since they are rotations about the
155 // positive z axis in 2d: i.e., they are done from the perspective of
156 // 'inside' a cell instead of the perspective of an abutting cell.
157 //
158 // For example: consider points on face 4 of a hexahedron with
159 // orientation 3. In 2d, rotating such points clockwise is the same as
160 // rotating them counter-clockwise from the perspective of the abutting
161 // face. Hence, such points must be rotated 90 degrees
162 // counter-clockwise.
163 switch (combined_orientation)
164 {
165 case 0:
166 return reflect(points);
167 case 2:
168 return rotate(reflect(points), 3);
169 case 4:
170 return rotate(reflect(points), 2);
171 case 6:
172 return rotate(reflect(points), 1);
173 case 1:
174 return points;
175 case 3:
176 return rotate(points, 1);
177 case 5:
178 return rotate(points, 2);
179 case 7:
180 return rotate(points, 3);
181 default:
182 DEAL_II_ASSERT_UNREACHABLE();
183 }
184 return {};
185 }
186
187 Quadrature<2>
188 mutate_quadrature(const Quadrature<2> &quadrature,
189 const unsigned char combined_orientation)
190 {
191 return Quadrature<2>(mutate_points_with_offset(quadrature.get_points(),
192 combined_orientation),
193 quadrature.get_weights());
194 }
195
196 std::pair<unsigned int, RefinementCase<2>>
197 select_subface_no_and_refinement_case(
198 const unsigned int subface_no,
199 const bool face_orientation,
200 const bool face_flip,
201 const bool face_rotation,
202 const internal::SubfaceCase<3> ref_case)
203 {
204 constexpr int dim = 3;
205 // for each subface of a given FaceRefineCase
206 // there is a corresponding equivalent
207 // subface number of one of the "standard"
208 // RefineCases (cut_x, cut_y, cut_xy). Map
209 // the given values to those equivalent
210 // ones.
211
212 // first, define an invalid number
213 static const unsigned int e = numbers::invalid_unsigned_int;
214
215 static const RefinementCase<dim - 1>
216 equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
217 [GeometryInfo<3>::max_children_per_face] = {
218 // case_none. there should be only
219 // invalid values here. However, as
220 // this function is also called (in
221 // tests) for cells which have no
222 // refined faces, use isotropic
223 // refinement instead
224 {RefinementCase<dim - 1>::cut_xy,
225 RefinementCase<dim - 1>::cut_xy,
226 RefinementCase<dim - 1>::cut_xy,
227 RefinementCase<dim - 1>::cut_xy},
228 // case_x
229 {RefinementCase<dim - 1>::cut_x,
230 RefinementCase<dim - 1>::cut_x,
231 RefinementCase<dim - 1>::no_refinement,
232 RefinementCase<dim - 1>::no_refinement},
233 // case_x1y
234 {RefinementCase<dim - 1>::cut_xy,
235 RefinementCase<dim - 1>::cut_xy,
236 RefinementCase<dim - 1>::cut_x,
237 RefinementCase<dim - 1>::no_refinement},
238 // case_x2y
239 {RefinementCase<dim - 1>::cut_x,
240 RefinementCase<dim - 1>::cut_xy,
241 RefinementCase<dim - 1>::cut_xy,
242 RefinementCase<dim - 1>::no_refinement},
243 // case_x1y2y
244 {RefinementCase<dim - 1>::cut_xy,
245 RefinementCase<dim - 1>::cut_xy,
246 RefinementCase<dim - 1>::cut_xy,
247 RefinementCase<dim - 1>::cut_xy},
248 // case_y
249 {RefinementCase<dim - 1>::cut_y,
250 RefinementCase<dim - 1>::cut_y,
251 RefinementCase<dim - 1>::no_refinement,
252 RefinementCase<dim - 1>::no_refinement},
253 // case_y1x
254 {RefinementCase<dim - 1>::cut_xy,
255 RefinementCase<dim - 1>::cut_xy,
256 RefinementCase<dim - 1>::cut_y,
257 RefinementCase<dim - 1>::no_refinement},
258 // case_y2x
259 {RefinementCase<dim - 1>::cut_y,
260 RefinementCase<dim - 1>::cut_xy,
261 RefinementCase<dim - 1>::cut_xy,
262 RefinementCase<dim - 1>::no_refinement},
263 // case_y1x2x
264 {RefinementCase<dim - 1>::cut_xy,
265 RefinementCase<dim - 1>::cut_xy,
266 RefinementCase<dim - 1>::cut_xy,
267 RefinementCase<dim - 1>::cut_xy},
268 // case_xy (case_isotropic)
269 {RefinementCase<dim - 1>::cut_xy,
270 RefinementCase<dim - 1>::cut_xy,
271 RefinementCase<dim - 1>::cut_xy,
272 RefinementCase<dim - 1>::cut_xy}};
273
274 static const unsigned int
275 equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic +
276 1][GeometryInfo<3>::max_children_per_face] =
277 {// case_none, see above
278 {0, 1, 2, 3},
279 // case_x
280 {0, 1, e, e},
281 // case_x1y
282 {0, 2, 1, e},
283 // case_x2y
284 {0, 1, 3, e},
285 // case_x1y2y
286 {0, 2, 1, 3},
287 // case_y
288 {0, 1, e, e},
289 // case_y1x
290 {0, 1, 1, e},
291 // case_y2x
292 {0, 2, 3, e},
293 // case_y1x2x
294 {0, 1, 2, 3},
295 // case_xy (case_isotropic)
296 {0, 1, 2, 3}};
297
298 // If face-orientation or face_rotation are
299 // non-standard, cut_x and cut_y have to be
300 // exchanged.
301 static const RefinementCase<dim - 1> ref_case_permutation[4] = {
302 RefinementCase<dim - 1>::no_refinement,
303 RefinementCase<dim - 1>::cut_y,
304 RefinementCase<dim - 1>::cut_x,
305 RefinementCase<dim - 1>::cut_xy};
306
307 // set a corresponding (equivalent)
308 // RefineCase and subface number
309 const RefinementCase<dim - 1> equ_ref_case =
310 equivalent_refine_case[ref_case][subface_no];
311 const unsigned int equ_subface_no =
312 equivalent_subface_number[ref_case][subface_no];
313 // make sure, that we got a valid subface and RefineCase
314 Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
315 ExcInternalError());
316 Assert(equ_subface_no != e, ExcInternalError());
317 // now, finally respect non-standard faces
318 const RefinementCase<dim - 1> final_ref_case =
319 (face_orientation == face_rotation ?
320 ref_case_permutation[equ_ref_case] :
321 equ_ref_case);
322
323 const unsigned int final_subface_no =
324 GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(
325 final_ref_case),
326 4,
327 equ_subface_no,
328 face_orientation,
329 face_flip,
330 face_rotation,
331 equ_ref_case);
332
333 return std::make_pair(final_subface_no, final_ref_case);
334 }
335 } // namespace
336 } // namespace QProjector
337} // namespace internal
338
339
340
341template <>
342void
343QProjector<1>::project_to_face(const ReferenceCell &reference_cell,
344 const Quadrature<0> &,
345 const unsigned int face_no,
346 std::vector<Point<1>> &q_points)
347{
348 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
349 (void)reference_cell;
350
351 const unsigned int dim = 1;
352 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
353 AssertDimension(q_points.size(), 1);
354
355 q_points[0] = Point<dim>(static_cast<double>(face_no));
356}
357
358
359
360template <>
361void
362QProjector<2>::project_to_face(const ReferenceCell &reference_cell,
363 const Quadrature<1> &quadrature,
364 const unsigned int face_no,
365 std::vector<Point<2>> &q_points)
366{
367 const unsigned int dim = 2;
368 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
369 Assert(q_points.size() == quadrature.size(),
370 ExcDimensionMismatch(q_points.size(), quadrature.size()));
371
372 if (reference_cell == ReferenceCells::Triangle)
373 {
374 // use linear polynomial to map the reference quadrature points correctly
375 // on faces, i.e., BarycentricPolynomials<1>(1)
376 for (unsigned int p = 0; p < quadrature.size(); ++p)
377 switch (face_no)
378 {
379 case 0:
380 q_points[p] = Point<dim>(quadrature.point(p)[0], 0);
381 break;
382 case 1:
383 q_points[p] =
384 Point<dim>(1 - quadrature.point(p)[0], quadrature.point(p)[0]);
385 break;
386 case 2:
387 q_points[p] = Point<dim>(0, 1 - quadrature.point(p)[0]);
388 break;
389 default:
390 DEAL_II_ASSERT_UNREACHABLE();
391 }
392 }
393 else if (reference_cell == ReferenceCells::Quadrilateral)
394 {
395 for (unsigned int p = 0; p < quadrature.size(); ++p)
396 switch (face_no)
397 {
398 case 0:
399 q_points[p] = Point<dim>(0, quadrature.point(p)[0]);
400 break;
401 case 1:
402 q_points[p] = Point<dim>(1, quadrature.point(p)[0]);
403 break;
404 case 2:
405 q_points[p] = Point<dim>(quadrature.point(p)[0], 0);
406 break;
407 case 3:
408 q_points[p] = Point<dim>(quadrature.point(p)[0], 1);
409 break;
410 default:
411 DEAL_II_ASSERT_UNREACHABLE();
412 }
413 }
414 else
415 {
416 DEAL_II_ASSERT_UNREACHABLE();
417 }
418}
419
420
421
422template <>
423void
424QProjector<3>::project_to_face(const ReferenceCell &reference_cell,
425 const Quadrature<2> &quadrature,
426 const unsigned int face_no,
427 std::vector<Point<3>> &q_points)
428{
429 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
430 (void)reference_cell;
431
432 AssertIndexRange(face_no, GeometryInfo<3>::faces_per_cell);
433 Assert(q_points.size() == quadrature.size(),
434 ExcDimensionMismatch(q_points.size(), quadrature.size()));
435 q_points.clear();
436 internal::QProjector::project_to_hex_face_and_append(quadrature.get_points(),
437 face_no,
438 q_points);
439}
440
441
442template <int dim>
443Quadrature<dim>
444QProjector<dim>::project_to_oriented_face(const ReferenceCell &reference_cell,
445 const Quadrature<dim - 1> &quadrature,
446 const unsigned int face_no,
447 const bool,
448 const bool,
449 const bool)
450{
451 return QProjector<dim>::project_to_face(reference_cell, quadrature, face_no);
452}
453
454
455
456template <>
457Quadrature<3>
458QProjector<3>::project_to_oriented_face(const ReferenceCell &reference_cell,
459 const Quadrature<2> &quadrature,
460 const unsigned int face_no,
461 const bool face_orientation,
462 const bool face_flip,
463 const bool face_rotation)
464{
465 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
466
467 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
468 quadrature,
469 internal::combined_face_orientation(face_orientation,
470 face_rotation,
471 face_flip));
472
473 return QProjector<3>::project_to_face(reference_cell, mutation, face_no);
474}
475
476
477
478template <>
479void
480QProjector<1>::project_to_subface(const ReferenceCell &reference_cell,
481 const Quadrature<0> &,
482 const unsigned int face_no,
483 const unsigned int,
484 std::vector<Point<1>> &q_points,
485 const RefinementCase<0> &)
486{
487 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
488 (void)reference_cell;
489
490 const unsigned int dim = 1;
491 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
492 AssertDimension(q_points.size(), 1);
493
494 q_points[0] = Point<dim>(static_cast<double>(face_no));
495}
496
497
498
499template <>
500void
501QProjector<2>::project_to_subface(const ReferenceCell &reference_cell,
502 const Quadrature<1> &quadrature,
503 const unsigned int face_no,
504 const unsigned int subface_no,
505 std::vector<Point<2>> &q_points,
506 const RefinementCase<1> &)
507{
508 const unsigned int dim = 2;
509 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
510 AssertIndexRange(subface_no, GeometryInfo<dim>::max_children_per_face);
511
512 Assert(q_points.size() == quadrature.size(),
513 ExcDimensionMismatch(q_points.size(), quadrature.size()));
514
515 if (reference_cell == ReferenceCells::Triangle)
516 {
517 // use linear polynomial to map the reference quadrature points correctly
518 // on faces, i.e., BarycentricPolynomials<1>(1)
519 for (unsigned int p = 0; p < quadrature.size(); ++p)
520 switch (face_no)
521 {
522 case 0:
523 switch (subface_no)
524 {
525 case 0:
526 q_points[p] = Point<dim>(quadrature.point(p)[0] / 2, 0);
527 break;
528 case 1:
529 q_points[p] =
530 Point<dim>(0.5 + quadrature.point(p)[0] / 2, 0);
531 break;
532 default:
533 DEAL_II_ASSERT_UNREACHABLE();
534 }
535 break;
536 case 1:
537 switch (subface_no)
538 {
539 case 0:
540 q_points[p] = Point<dim>(1 - quadrature.point(p)[0] / 2,
541 quadrature.point(p)[0] / 2);
542 break;
543 case 1:
544 q_points[p] = Point<dim>(0.5 - quadrature.point(p)[0] / 2,
545 0.5 + quadrature.point(p)[0] / 2);
546 break;
547 default:
548 DEAL_II_ASSERT_UNREACHABLE();
549 }
550 break;
551 case 2:
552 switch (subface_no)
553 {
554 case 0:
555 q_points[p] = Point<dim>(0, 1 - quadrature.point(p)[0] / 2);
556 break;
557 case 1:
558 q_points[p] =
559 Point<dim>(0, 0.5 - quadrature.point(p)[0] / 2);
560 break;
561 default:
562 DEAL_II_ASSERT_UNREACHABLE();
563 }
564 break;
565 default:
566 DEAL_II_ASSERT_UNREACHABLE();
567 }
568 }
569 else if (reference_cell == ReferenceCells::Quadrilateral)
570 {
571 for (unsigned int p = 0; p < quadrature.size(); ++p)
572 switch (face_no)
573 {
574 case 0:
575 switch (subface_no)
576 {
577 case 0:
578 q_points[p] = Point<dim>(0, quadrature.point(p)[0] / 2);
579 break;
580 case 1:
581 q_points[p] =
582 Point<dim>(0, quadrature.point(p)[0] / 2 + 0.5);
583 break;
584 default:
585 DEAL_II_ASSERT_UNREACHABLE();
586 }
587 break;
588 case 1:
589 switch (subface_no)
590 {
591 case 0:
592 q_points[p] = Point<dim>(1, quadrature.point(p)[0] / 2);
593 break;
594 case 1:
595 q_points[p] =
596 Point<dim>(1, quadrature.point(p)[0] / 2 + 0.5);
597 break;
598 default:
599 DEAL_II_ASSERT_UNREACHABLE();
600 }
601 break;
602 case 2:
603 switch (subface_no)
604 {
605 case 0:
606 q_points[p] = Point<dim>(quadrature.point(p)[0] / 2, 0);
607 break;
608 case 1:
609 q_points[p] =
610 Point<dim>(quadrature.point(p)[0] / 2 + 0.5, 0);
611 break;
612 default:
613 DEAL_II_ASSERT_UNREACHABLE();
614 }
615 break;
616 case 3:
617 switch (subface_no)
618 {
619 case 0:
620 q_points[p] = Point<dim>(quadrature.point(p)[0] / 2, 1);
621 break;
622 case 1:
623 q_points[p] =
624 Point<dim>(quadrature.point(p)[0] / 2 + 0.5, 1);
625 break;
626 default:
627 DEAL_II_ASSERT_UNREACHABLE();
628 }
629 break;
630
631 default:
632 DEAL_II_ASSERT_UNREACHABLE();
633 }
634 }
635 else
636 {
637 DEAL_II_ASSERT_UNREACHABLE();
638 }
639}
640
641
642
643template <>
644void
645QProjector<3>::project_to_subface(const ReferenceCell &reference_cell,
646 const Quadrature<2> &quadrature,
647 const unsigned int face_no,
648 const unsigned int subface_no,
649 std::vector<Point<3>> &q_points,
650 const RefinementCase<2> &ref_case)
651{
652 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
653 (void)reference_cell;
654
655 AssertIndexRange(face_no, GeometryInfo<3>::faces_per_cell);
656 AssertIndexRange(subface_no, GeometryInfo<3>::max_children_per_face);
657 Assert(q_points.size() == quadrature.size(),
658 ExcDimensionMismatch(q_points.size(), quadrature.size()));
659
660 q_points.clear();
661 internal::QProjector::project_to_hex_face_and_append(
662 quadrature.get_points(), face_no, q_points, ref_case, subface_no);
663}
664
665
666
667template <int dim>
668Quadrature<dim>
669QProjector<dim>::project_to_oriented_subface(
670 const ReferenceCell &reference_cell,
671 const Quadrature<dim - 1> &quadrature,
672 const unsigned int face_no,
673 const unsigned int subface_no,
674 const bool,
675 const bool,
676 const bool,
677 const internal::SubfaceCase<dim>)
678{
679 return QProjector<dim>::project_to_subface(
680 reference_cell,
681 quadrature,
682 face_no,
683 subface_no,
684 RefinementCase<dim - 1>::isotropic_refinement);
685}
686
687
688
689template <>
690Quadrature<3>
691QProjector<3>::project_to_oriented_subface(
692 const ReferenceCell &reference_cell,
693 const Quadrature<2> &quadrature,
694 const unsigned int face_no,
695 const unsigned int subface_no,
696 const bool face_orientation,
697 const bool face_flip,
698 const bool face_rotation,
699 const internal::SubfaceCase<3> ref_case)
700{
701 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
702
703 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
704 quadrature,
705 internal::combined_face_orientation(face_orientation,
706 face_rotation,
707 face_flip));
708
709 const std::pair<unsigned int, RefinementCase<2>>
710 final_subface_no_and_ref_case =
711 internal::QProjector::select_subface_no_and_refinement_case(
712 subface_no, face_orientation, face_flip, face_rotation, ref_case);
713
714 return QProjector<3>::project_to_subface(
715 reference_cell,
716 mutation,
717 face_no,
718 final_subface_no_and_ref_case.first,
719 final_subface_no_and_ref_case.second);
720}
721
722
723
724template <>
725Quadrature<1>
726QProjector<1>::project_to_all_faces(const ReferenceCell &reference_cell,
727 const hp::QCollection<0> &quadrature)
728{
729 AssertDimension(quadrature.size(), 1);
730 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
731 (void)reference_cell;
732
733 const unsigned int dim = 1;
734
735 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;
736
737 // first fix quadrature points
738 std::vector<Point<dim>> q_points;
739 q_points.reserve(n_points * n_faces);
740 std::vector<Point<dim>> help(n_points);
741
742
743 // project to each face and append
744 // results
745 for (unsigned int face = 0; face < n_faces; ++face)
746 {
747 project_to_face(reference_cell,
748 quadrature[quadrature.size() == 1 ? 0 : face],
749 face,
750 help);
751 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
752 }
753
754 // next copy over weights
755 std::vector<double> weights;
756 weights.reserve(n_points * n_faces);
757 for (unsigned int face = 0; face < n_faces; ++face)
758 std::copy(
759 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().begin(),
760 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().end(),
761 std::back_inserter(weights));
762
763 Assert(q_points.size() == n_points * n_faces, ExcInternalError());
764 Assert(weights.size() == n_points * n_faces, ExcInternalError());
765
766 return Quadrature<dim>(std::move(q_points), std::move(weights));
767}
768
769
770
771template <>
772Quadrature<2>
773QProjector<2>::project_to_all_faces(const ReferenceCell &reference_cell,
774 const hp::QCollection<1> &quadrature)
775{
776 if (reference_cell == ReferenceCells::Triangle)
777 {
778 const auto support_points_line =
779 [](const auto &face, const auto &orientation) -> std::vector<Point<2>> {
780 // MSVC struggles when using face.first.begin()
781 const Point<2, double> *vertices_ptr = &face.first[0];
782 ArrayView<const Point<2>> vertices(vertices_ptr, face.first.size());
783 const auto temp =
784 ReferenceCells::Line.permute_by_combined_orientation(vertices,
785 orientation);
786 return std::vector<Point<2>>(temp.begin(),
787 temp.begin() + face.first.size());
788 };
789
790 // reference faces (defined by its support points and arc length)
791 const std::array<std::pair<std::array<Point<2>, 2>, double>, 3> faces = {
792 {{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0},
793 {{{Point<2>(1.0, 0.0), Point<2>(0.0, 1.0)}}, std::sqrt(2.0)},
794 {{{Point<2>(0.0, 1.0), Point<2>(0.0, 0.0)}}, 1.0}}};
795
796 // linear polynomial to map the reference quadrature points correctly
797 // on faces
798 const auto poly = BarycentricPolynomials<1>::get_fe_p_basis(1);
799
800 // new (projected) quadrature points and weights
801 std::vector<Point<2>> points;
802 std::vector<double> weights;
803
804 // loop over all faces (lines) ...
805 for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
806 // ... and over all possible orientations
807 for (unsigned int orientation = 0; orientation < 2; ++orientation)
808 {
809 const auto &face = faces[face_no];
810
811 // determine support point of the current line with the correct
812 // orientation
813 std::vector<Point<2>> support_points =
814 support_points_line(face, orientation);
815
816 // the quadrature rule to be projected ...
817 const auto &sub_quadrature_points =
818 quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
819 const auto &sub_quadrature_weights =
820 quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
821
822 // loop over all quadrature points
823 for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
824 {
825 Point<2> mapped_point;
826
827 // map reference quadrature point
828 for (unsigned int i = 0; i < 2; ++i)
829 mapped_point +=
830 support_points[i] *
831 poly.compute_value(i, sub_quadrature_points[j]);
832
833 points.emplace_back(mapped_point);
834
835 // scale weight by arc length
836 weights.emplace_back(sub_quadrature_weights[j] * face.second);
837 }
838 }
839
840 // construct new quadrature rule
841 return Quadrature<2>(std::move(points), std::move(weights));
842 }
843
844 Assert(reference_cell == ReferenceCells::Quadrilateral, ExcNotImplemented());
845
846 const unsigned int dim = 2;
847
848 const unsigned int n_faces = GeometryInfo<dim>::faces_per_cell;
849
850 unsigned int n_points_total = 0;
851
852 if (quadrature.size() == 1)
853 n_points_total = quadrature[0].size() * GeometryInfo<dim>::faces_per_cell;
854 else
855 {
856 AssertDimension(quadrature.size(), GeometryInfo<dim>::faces_per_cell);
857 for (const auto &q : quadrature)
858 n_points_total += q.size();
859 }
860
861 // first fix quadrature points
862 std::vector<Point<dim>> q_points;
863 q_points.reserve(n_points_total);
864 std::vector<Point<dim>> help;
865 help.reserve(quadrature.max_n_quadrature_points());
866
867 // project to each face and append
868 // results
869 for (unsigned int face = 0; face < n_faces; ++face)
870 {
871 help.resize(quadrature[quadrature.size() == 1 ? 0 : face].size());
872 project_to_face(reference_cell,
873 quadrature[quadrature.size() == 1 ? 0 : face],
874 face,
875 help);
876 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
877 }
878
879 // next copy over weights
880 std::vector<double> weights;
881 weights.reserve(n_points_total);
882 for (unsigned int face = 0; face < n_faces; ++face)
883 std::copy(
884 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().begin(),
885 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().end(),
886 std::back_inserter(weights));
887
888 Assert(q_points.size() == n_points_total, ExcInternalError());
889 Assert(weights.size() == n_points_total, ExcInternalError());
890
891 return Quadrature<dim>(std::move(q_points), std::move(weights));
892}
893
894
895
896template <>
897Quadrature<3>
898QProjector<3>::project_to_all_faces(const ReferenceCell &reference_cell,
899 const hp::QCollection<2> &quadrature)
900{
901 const auto process = [&](const std::vector<std::vector<Point<3>>> &faces) {
902 // new (projected) quadrature points and weights
903 std::vector<Point<3>> points;
904 std::vector<double> weights;
905
906 const auto poly_tri = BarycentricPolynomials<2>::get_fe_p_basis(1);
907 const TensorProductPolynomials<2> poly_quad(
908 Polynomials::generate_complete_Lagrange_basis(
909 {Point<1>(0.0), Point<1>(1.0)}));
910
911 // loop over all faces (triangles) ...
912 for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
913 {
914 const ReferenceCell face_reference_cell =
915 reference_cell.face_reference_cell(face_no);
916 // We will use linear polynomials to map the reference quadrature
917 // points correctly to on faces. There are as many linear shape
918 // functions as there are vertices in the face.
919 const unsigned int n_linear_shape_functions = faces[face_no].size();
920 std::vector<Tensor<1, 2>> shape_derivatives;
921
922 const auto &poly =
923 (n_linear_shape_functions == 3 ?
924 static_cast<const ScalarPolynomialsBase<2> &>(poly_tri) :
925 static_cast<const ScalarPolynomialsBase<2> &>(poly_quad));
926
927 // ... and over all possible orientations
928 for (unsigned char orientation = 0;
929 orientation < reference_cell.n_face_orientations(face_no);
930 ++orientation)
931 {
932 const auto &face = faces[face_no];
933
934 // The goal of this function is to compute identical sets of
935 // quadrature points on the common face of two abutting cells. Our
936 // orientation convention is that, given such a pair of abutting
937 // cells:
938 //
939 // 1. The shared face, from the perspective of the first cell, is
940 // in the default orientation.
941 // 2. The shared face, from the perspective of the second cell, has
942 // its orientation computed relative to the first cell: i.e.,
943 // 'orientation' is the vertex permutation applied to the first
944 // cell's face to get the second cell's face.
945 //
946 // The first case is trivial since points do not need to be
947 // oriented. However, in the second case, we need to use the
948 // *reverse* of the stored orientation (i.e., the permutation
949 // applied to the second cell's face which yields the first cell's
950 // face) so that we get identical quadrature points.
951 //
952 // For more information see connectivity.h.
953 const boost::container::small_vector<Point<3>, 8> support_points =
954 face_reference_cell.permute_by_combined_orientation<Point<3>>(
955 face,
956 face_reference_cell.get_inverse_combined_orientation(
957 orientation));
958
959 // the quadrature rule to be projected ...
960 const auto &sub_quadrature_points =
961 quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
962 const auto &sub_quadrature_weights =
963 quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
964
965 // loop over all quadrature points
966 for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
967 {
968 Point<3> mapped_point;
969
970 // map reference quadrature point
971 for (unsigned int i = 0; i < n_linear_shape_functions; ++i)
972 mapped_point +=
973 support_points[i] *
974 poly.compute_value(i, sub_quadrature_points[j]);
975
976 points.push_back(mapped_point);
977
978 // scale quadrature weight
979 const double scaling = [&]() {
980 const unsigned int dim_ = 2;
981 const unsigned int spacedim = 3;
982
983 Tensor<1, dim_, Tensor<1, spacedim>> DX_t;
984
985 shape_derivatives.resize(n_linear_shape_functions);
986
987 for (unsigned int i = 0; i < n_linear_shape_functions; ++i)
988 shape_derivatives[i] =
989 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
990
991 for (unsigned int k = 0; k < n_linear_shape_functions; ++k)
992 for (unsigned int i = 0; i < spacedim; ++i)
993 for (unsigned int j = 0; j < dim_; ++j)
994 DX_t[j][i] +=
995 shape_derivatives[k][j] * support_points[k][i];
996
997 Tensor<2, dim_> G;
998 for (unsigned int i = 0; i < dim_; ++i)
999 for (unsigned int j = 0; j < dim_; ++j)
1000 G[i][j] = DX_t[i] * DX_t[j];
1001
1002 return std::sqrt(determinant(G));
1003 }();
1004
1005 weights.push_back(sub_quadrature_weights[j] * scaling);
1006 }
1007 }
1008 }
1009
1010 // construct new quadrature rule
1011 return Quadrature<3>(std::move(points), std::move(weights));
1012 };
1013
1014 if ((reference_cell == ReferenceCells::Tetrahedron) ||
1015 (reference_cell == ReferenceCells::Wedge) ||
1016 (reference_cell == ReferenceCells::Pyramid))
1017 {
1018 std::vector<std::vector<Point<3>>> face_vertex_locations(
1019 reference_cell.n_faces());
1020 for (const unsigned int f : reference_cell.face_indices())
1021 {
1022 face_vertex_locations[f].resize(
1023 reference_cell.face_reference_cell(f).n_vertices());
1024 for (const unsigned int v :
1025 reference_cell.face_reference_cell(f).vertex_indices())
1026 face_vertex_locations[f][v] =
1027 reference_cell.face_vertex_location<3>(f, v);
1028 }
1029
1030 return process(face_vertex_locations);
1031 }
1032 else
1033 {
1034 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
1035
1036 const unsigned int dim = 3;
1037
1038 unsigned int n_points_total = 0;
1039
1040 if (quadrature.size() == 1)
1041 n_points_total =
1042 quadrature[0].size() * GeometryInfo<dim>::faces_per_cell;
1043 else
1044 {
1045 AssertDimension(quadrature.size(), GeometryInfo<dim>::faces_per_cell);
1046 for (const auto &q : quadrature)
1047 n_points_total += q.size();
1048 }
1049
1050 n_points_total *= 8;
1051
1052 // first fix quadrature points
1053 std::vector<Point<dim>> q_points;
1054 q_points.reserve(n_points_total);
1055
1056 std::vector<double> weights;
1057 weights.reserve(n_points_total);
1058
1059 for (unsigned char offset = 0; offset < 8; ++offset)
1060 {
1061 const auto mutation = internal::QProjector::mutate_points_with_offset(
1062 quadrature[0].get_points(), offset);
1063 // project to each face and append results
1064 for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
1065 ++face)
1066 {
1067 const unsigned int q_index = quadrature.size() == 1 ? 0 : face;
1068
1069 internal::QProjector::project_to_hex_face_and_append(
1070 q_index > 0 ? internal::QProjector::mutate_points_with_offset(
1071 quadrature[face].get_points(), offset) :
1072 mutation,
1073 face,
1074 q_points);
1075
1076 std::copy(quadrature[q_index].get_weights().begin(),
1077 quadrature[q_index].get_weights().end(),
1078 std::back_inserter(weights));
1079 }
1080 }
1081
1082 Assert(q_points.size() == n_points_total, ExcInternalError());
1083 Assert(weights.size() == n_points_total, ExcInternalError());
1084
1085 return Quadrature<dim>(std::move(q_points), std::move(weights));
1086 }
1087}
1088
1089
1090
1091template <>
1092Quadrature<1>
1093QProjector<1>::project_to_all_subfaces(const ReferenceCell &reference_cell,
1094 const Quadrature<0> &quadrature)
1095{
1096 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1097 (void)reference_cell;
1098
1099 const unsigned int dim = 1;
1100
1101 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
1102 subfaces_per_face =
1103 GeometryInfo<dim>::max_children_per_face;
1104
1105 // first fix quadrature points
1106 std::vector<Point<dim>> q_points;
1107 q_points.reserve(n_points * n_faces * subfaces_per_face);
1108 std::vector<Point<dim>> help(n_points);
1109
1110 // project to each face and copy
1111 // results
1112 for (unsigned int face = 0; face < n_faces; ++face)
1113 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1114 {
1115 project_to_subface(reference_cell, quadrature, face, subface, help);
1116 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1117 }
1118
1119 // next copy over weights
1120 std::vector<double> weights;
1121 weights.reserve(n_points * n_faces * subfaces_per_face);
1122 for (unsigned int face = 0; face < n_faces; ++face)
1123 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1124 std::copy(quadrature.get_weights().begin(),
1125 quadrature.get_weights().end(),
1126 std::back_inserter(weights));
1127
1128 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1129 ExcInternalError());
1130 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1131 ExcInternalError());
1132
1133 return Quadrature<dim>(std::move(q_points), std::move(weights));
1134}
1135
1136
1137
1138template <>
1139Quadrature<2>
1140QProjector<2>::project_to_all_subfaces(const ReferenceCell &reference_cell,
1141 const SubQuadrature &quadrature)
1142{
1143 if (reference_cell == ReferenceCells::Triangle ||
1144 reference_cell == ReferenceCells::Tetrahedron)
1145 return Quadrature<2>(); // nothing to do
1146
1147 Assert(reference_cell == ReferenceCells::Quadrilateral, ExcNotImplemented());
1148
1149 const unsigned int dim = 2;
1150
1151 const unsigned int n_points = quadrature.size(),
1152 n_faces = GeometryInfo<dim>::faces_per_cell,
1153 subfaces_per_face =
1154 GeometryInfo<dim>::max_children_per_face;
1155
1156 // first fix quadrature points
1157 std::vector<Point<dim>> q_points;
1158 q_points.reserve(n_points * n_faces * subfaces_per_face);
1159 std::vector<Point<dim>> help(n_points);
1160
1161 // project to each face and copy
1162 // results
1163 for (unsigned int face = 0; face < n_faces; ++face)
1164 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1165 {
1166 project_to_subface(reference_cell, quadrature, face, subface, help);
1167 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1168 }
1169
1170 // next copy over weights
1171 std::vector<double> weights;
1172 weights.reserve(n_points * n_faces * subfaces_per_face);
1173 for (unsigned int face = 0; face < n_faces; ++face)
1174 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1175 std::copy(quadrature.get_weights().begin(),
1176 quadrature.get_weights().end(),
1177 std::back_inserter(weights));
1178
1179 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1180 ExcInternalError());
1181 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1182 ExcInternalError());
1183
1184 return Quadrature<dim>(std::move(q_points), std::move(weights));
1185}
1186
1187
1188
1189template <>
1190Quadrature<3>
1191QProjector<3>::project_to_all_subfaces(const ReferenceCell &reference_cell,
1192 const SubQuadrature &quadrature)
1193{
1194 if (reference_cell == ReferenceCells::Triangle ||
1195 reference_cell == ReferenceCells::Tetrahedron)
1196 return Quadrature<3>(); // nothing to do
1197
1198 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
1199
1200 const unsigned int dim = 3;
1201
1202 const unsigned int n_points = quadrature.size(),
1203 n_faces = GeometryInfo<dim>::faces_per_cell,
1204 total_subfaces_per_face = 2 + 2 + 4;
1205
1206 // first fix quadrature points
1207 std::vector<Point<dim>> q_points;
1208 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1209
1210 std::vector<double> weights;
1211 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1212
1213 // do the following for all possible mutations of a face (mutation==0
1214 // corresponds to a face with standard orientation, no flip and no rotation)
1215 for (unsigned char offset = 0; offset < 8; ++offset)
1216 {
1217 const auto mutation =
1218 internal::QProjector::mutate_points_with_offset(quadrature.get_points(),
1219 offset);
1220
1221 // project to each face and copy results
1222 for (unsigned int face = 0; face < n_faces; ++face)
1223 for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1224 ref_case >= RefinementCase<dim - 1>::cut_x;
1225 --ref_case)
1226 for (unsigned int subface = 0;
1227 subface < GeometryInfo<dim - 1>::n_children(
1228 RefinementCase<dim - 1>(ref_case));
1229 ++subface)
1230 {
1231 internal::QProjector::project_to_hex_face_and_append(
1232 mutation,
1233 face,
1234 q_points,
1235 RefinementCase<dim - 1>(ref_case),
1236 subface);
1237
1238 // next copy over weights
1239 std::copy(quadrature.get_weights().begin(),
1240 quadrature.get_weights().end(),
1241 std::back_inserter(weights));
1242 }
1243 }
1244
1245 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1246 ExcInternalError());
1247 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1248 ExcInternalError());
1249
1250 return Quadrature<dim>(std::move(q_points), std::move(weights));
1251}
1252
1253
1254
1255template <int dim>
1256Quadrature<dim>
1257QProjector<dim>::project_to_child(const ReferenceCell &reference_cell,
1258 const Quadrature<dim> &quadrature,
1259 const unsigned int child_no)
1260{
1261 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1262 ExcNotImplemented());
1263 (void)reference_cell;
1264
1265 AssertIndexRange(child_no, GeometryInfo<dim>::max_children_per_cell);
1266
1267 const unsigned int n_q_points = quadrature.size();
1268
1269 std::vector<Point<dim>> q_points(n_q_points);
1270 for (unsigned int i = 0; i < n_q_points; ++i)
1271 q_points[i] =
1272 GeometryInfo<dim>::child_to_cell_coordinates(quadrature.point(i),
1273 child_no);
1274
1275 // for the weights, things are
1276 // equally simple: copy them and
1277 // scale them
1278 std::vector<double> weights = quadrature.get_weights();
1279 for (unsigned int i = 0; i < n_q_points; ++i)
1280 weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1281
1282 return Quadrature<dim>(q_points, weights);
1283}
1284
1285
1286
1287template <int dim>
1288Quadrature<dim>
1289QProjector<dim>::project_to_all_children(const ReferenceCell &reference_cell,
1290 const Quadrature<dim> &quadrature)
1291{
1292 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1293 ExcNotImplemented());
1294 (void)reference_cell;
1295
1296 const unsigned int n_points = quadrature.size(),
1297 n_children = GeometryInfo<dim>::max_children_per_cell;
1298
1299 std::vector<Point<dim>> q_points(n_points * n_children);
1300 std::vector<double> weights(n_points * n_children);
1301
1302 // project to each child and copy
1303 // results
1304 for (unsigned int child = 0; child < n_children; ++child)
1305 {
1306 Quadrature<dim> help =
1307 project_to_child(reference_cell, quadrature, child);
1308 for (unsigned int i = 0; i < n_points; ++i)
1309 {
1310 q_points[child * n_points + i] = help.point(i);
1311 weights[child * n_points + i] = help.weight(i);
1312 }
1313 }
1314 return Quadrature<dim>(q_points, weights);
1315}
1316
1317
1318
1319template <int dim>
1320Quadrature<dim>
1321QProjector<dim>::project_to_line(const ReferenceCell &reference_cell,
1322 const Quadrature<1> &quadrature,
1323 const Point<dim> &p1,
1324 const Point<dim> &p2)
1325{
1326 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1327 ExcNotImplemented());
1328 (void)reference_cell;
1329
1330 const unsigned int n = quadrature.size();
1331 std::vector<Point<dim>> points(n);
1332 std::vector<double> weights(n);
1333 const double length = p1.distance(p2);
1334
1335 for (unsigned int k = 0; k < n; ++k)
1336 {
1337 const double alpha = quadrature.point(k)[0];
1338 points[k] = alpha * p2;
1339 points[k] += (1. - alpha) * p1;
1340 weights[k] = length * quadrature.weight(k);
1341 }
1342 return Quadrature<dim>(points, weights);
1343}
1344
1345
1346
1347template <int dim>
1348typename QProjector<dim>::DataSetDescriptor
1349QProjector<dim>::DataSetDescriptor::face(const ReferenceCell &reference_cell,
1350 const unsigned int face_no,
1351 const bool face_orientation,
1352 const bool face_flip,
1353 const bool face_rotation,
1354 const unsigned int n_quadrature_points)
1355{
1356 return face(reference_cell,
1357 face_no,
1358 internal::combined_face_orientation(face_orientation,
1359 face_rotation,
1360 face_flip),
1361 n_quadrature_points);
1362}
1363
1364
1365
1366template <int dim>
1367typename QProjector<dim>::DataSetDescriptor
1368QProjector<dim>::DataSetDescriptor::face(
1369 const ReferenceCell &reference_cell,
1370 const unsigned int face_no,
1371 const unsigned char combined_orientation,
1372 const unsigned int n_quadrature_points)
1373{
1374 // TODO: once we move to representing the default orientation as 0 (instead of
1375 // 1) we can get rid of the dim = 1 check
1376 Assert(dim == 1 ||
1377 (combined_orientation < reference_cell.n_face_orientations(face_no)),
1378 ExcInternalError());
1379 if (reference_cell == ReferenceCells::Triangle ||
1380 reference_cell == ReferenceCells::Tetrahedron)
1381 {
1382 if (dim == 2)
1383 return {(2 * face_no +
1384 (combined_orientation ==
1385 ReferenceCell::default_combined_face_orientation() ?
1386 1 :
1387 0)) *
1388 n_quadrature_points};
1389 else if (dim == 3)
1390 {
1391 return {(reference_cell.n_face_orientations(face_no) * face_no +
1392 combined_orientation) *
1393 n_quadrature_points};
1394 }
1395 }
1396
1397 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1398 ExcNotImplemented());
1399
1400 Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
1401
1402 switch (dim)
1403 {
1404 case 1:
1405 case 2:
1406 return face_no * n_quadrature_points;
1407 case 3:
1408 return (face_no +
1409 GeometryInfo<dim>::faces_per_cell * combined_orientation) *
1410 n_quadrature_points;
1411 default:
1412 DEAL_II_ASSERT_UNREACHABLE();
1413 }
1414 return numbers::invalid_unsigned_int;
1415}
1416
1417
1418
1419template <int dim>
1420typename QProjector<dim>::DataSetDescriptor
1421QProjector<dim>::DataSetDescriptor::face(
1422 const ReferenceCell &reference_cell,
1423 const unsigned int face_no,
1424 const bool face_orientation,
1425 const bool face_flip,
1426 const bool face_rotation,
1427 const hp::QCollection<dim - 1> &quadrature)
1428{
1429 return face(reference_cell,
1430 face_no,
1431 internal::combined_face_orientation(face_orientation,
1432 face_rotation,
1433 face_flip),
1434 quadrature);
1435}
1436
1437
1438
1439template <int dim>
1440typename QProjector<dim>::DataSetDescriptor
1441QProjector<dim>::DataSetDescriptor::face(
1442 const ReferenceCell &reference_cell,
1443 const unsigned int face_no,
1444 const unsigned char combined_orientation,
1445 const hp::QCollection<dim - 1> &quadrature)
1446{
1447 if (reference_cell == ReferenceCells::Triangle ||
1448 reference_cell == ReferenceCells::Tetrahedron ||
1449 reference_cell == ReferenceCells::Wedge ||
1450 reference_cell == ReferenceCells::Pyramid)
1451 {
1452 unsigned int offset = 0;
1453 Assert(combined_orientation < reference_cell.n_face_orientations(face_no),
1454 ExcInternalError());
1455
1456 static const unsigned int X = numbers::invalid_unsigned_int;
1457 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1458 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1459 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1460 static const std::array<unsigned int, 5> scale_pyramid = {
1461 {8, 6, 6, 6, 6}};
1462
1463 const auto &scale =
1464 (reference_cell == ReferenceCells::Triangle) ?
1465 scale_tri :
1466 ((reference_cell == ReferenceCells::Tetrahedron) ?
1467 scale_tet :
1468 ((reference_cell == ReferenceCells::Wedge) ? scale_wedge :
1469 scale_pyramid));
1470
1471 if (quadrature.size() == 1)
1472 offset = scale[0] * quadrature[0].size() * face_no;
1473 else
1474 for (unsigned int i = 0; i < face_no; ++i)
1475 offset += scale[i] * quadrature[i].size();
1476
1477 if (dim == 2)
1478 return {offset +
1479 (combined_orientation ==
1480 ReferenceCell::default_combined_face_orientation()) *
1481 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1482 else if (dim == 3)
1483 {
1484 return {offset +
1485 combined_orientation *
1486 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1487 }
1488 }
1489
1490 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1491 ExcNotImplemented());
1492
1493 Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
1494
1495 switch (dim)
1496 {
1497 case 1:
1498 case 2:
1499 {
1500 if (quadrature.size() == 1)
1501 return quadrature[0].size() * face_no;
1502 else
1503 {
1504 unsigned int result = 0;
1505 for (unsigned int i = 0; i < face_no; ++i)
1506 result += quadrature[i].size();
1507 return result;
1508 }
1509 }
1510 case 3:
1511 {
1512 if (quadrature.size() == 1)
1513 return (face_no +
1514 combined_orientation * GeometryInfo<dim>::faces_per_cell) *
1515 quadrature[0].size();
1516 else
1517 {
1518 unsigned int n_points_i = 0;
1519 for (unsigned int i = 0; i < face_no; ++i)
1520 n_points_i += quadrature[i].size();
1521
1522 unsigned int n_points = 0;
1523 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1524 ++i)
1525 n_points += quadrature[i].size();
1526
1527 return n_points_i + combined_orientation * n_points;
1528 }
1529 }
1530
1531 default:
1532 DEAL_II_ASSERT_UNREACHABLE();
1533 }
1534 return numbers::invalid_unsigned_int;
1535}
1536
1537
1538
1539template <int dim>
1540typename QProjector<dim>::DataSetDescriptor
1541QProjector<dim>::DataSetDescriptor::subface(
1542 const ReferenceCell &reference_cell,
1543 const unsigned int face_no,
1544 const unsigned int subface_no,
1545 const bool face_orientation,
1546 const bool face_flip,
1547 const bool face_rotation,
1548 const unsigned int n_quadrature_points,
1549 const internal::SubfaceCase<dim> ref_case)
1550{
1551 return QProjector<dim>::DataSetDescriptor::subface(
1552 reference_cell,
1553 face_no,
1554 subface_no,
1555 internal::combined_face_orientation(face_orientation,
1556 face_rotation,
1557 face_flip),
1558 n_quadrature_points,
1559 ref_case);
1560}
1561
1562
1563
1564template <>
1565QProjector<1>::DataSetDescriptor
1566QProjector<1>::DataSetDescriptor::subface(
1567 const ReferenceCell &reference_cell,
1568 const unsigned int face_no,
1569 const unsigned int subface_no,
1570 const unsigned char /*combined_orientation*/,
1571 const unsigned int n_quadrature_points,
1572 const internal::SubfaceCase<1>)
1573{
1574 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1575 (void)reference_cell;
1576
1577 Assert(face_no < GeometryInfo<1>::faces_per_cell, ExcInternalError());
1578 Assert(subface_no < GeometryInfo<1>::max_children_per_face,
1579 ExcInternalError());
1580
1581 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1582 n_quadrature_points);
1583}
1584
1585
1586
1587template <>
1588QProjector<2>::DataSetDescriptor
1589QProjector<2>::DataSetDescriptor::subface(
1590 const ReferenceCell &reference_cell,
1591 const unsigned int face_no,
1592 const unsigned int subface_no,
1593 const unsigned char /*combined_orientation*/,
1594 const unsigned int n_quadrature_points,
1595 const internal::SubfaceCase<2>)
1596{
1597 Assert(reference_cell == ReferenceCells::Quadrilateral, ExcNotImplemented());
1598 (void)reference_cell;
1599
1600 Assert(face_no < GeometryInfo<2>::faces_per_cell, ExcInternalError());
1601 Assert(subface_no < GeometryInfo<2>::max_children_per_face,
1602 ExcInternalError());
1603
1604 return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1605 n_quadrature_points);
1606}
1607
1608
1609
1610template <>
1611QProjector<3>::DataSetDescriptor
1612QProjector<3>::DataSetDescriptor::subface(
1613 const ReferenceCell &reference_cell,
1614 const unsigned int face_no,
1615 const unsigned int subface_no,
1616 const unsigned char combined_orientation,
1617 const unsigned int n_quadrature_points,
1618 const internal::SubfaceCase<3> ref_case)
1619{
1620 const unsigned int dim = 3;
1621
1622 const auto [face_orientation, face_rotation, face_flip] =
1623 internal::split_face_orientation(combined_orientation);
1624
1625 Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
1626 (void)reference_cell;
1627
1628 Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
1629 Assert(subface_no < GeometryInfo<dim>::max_children_per_face,
1630 ExcInternalError());
1631
1632 // in 3d, we have to account for faces that
1633 // have non-standard orientation. thus, we
1634 // have to store _eight_ data sets per face
1635 // or subface already for the isotropic
1636 // case. Additionally, we have three
1637 // different refinement cases, resulting in
1638 // <tt>4 + 2 + 2 = 8</tt> different subfaces
1639 // for each face.
1640 const unsigned int total_subfaces_per_face = 8;
1641
1642 // set up a table with the offsets for a
1643 // given refinement case respecting the
1644 // corresponding number of subfaces. the
1645 // index corresponds to (RefineCase::Type - 1)
1646
1647 // note, that normally we should use the
1648 // obvious offsets 0,2,6. However, prior to
1649 // the implementation of anisotropic
1650 // refinement, in many places of the library
1651 // the convention was used, that the first
1652 // dataset with offset 0 corresponds to a
1653 // standard (isotropic) face
1654 // refinement. therefore we use the offsets
1655 // 6,4,0 here to stick to that (implicit)
1656 // convention
1657 static const unsigned int ref_case_offset[3] = {
1658 6, // cut_x
1659 4, // cut_y
1660 0 // cut_xy
1661 };
1662
1663 const std::pair<unsigned int, RefinementCase<2>>
1664 final_subface_no_and_ref_case =
1665 internal::QProjector::select_subface_no_and_refinement_case(
1666 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1667
1668 return ((face_no * total_subfaces_per_face +
1669 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1670 final_subface_no_and_ref_case.first) +
1671 GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face *
1672 combined_orientation) *
1673 n_quadrature_points;
1674}
1675
1676
1677
1678template <int dim>
1679Quadrature<dim>
1680QProjector<dim>::project_to_face(const ReferenceCell &reference_cell,
1681 const SubQuadrature &quadrature,
1682 const unsigned int face_no)
1683{
1684 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1685 ExcNotImplemented());
1686 (void)reference_cell;
1687
1688 std::vector<Point<dim>> points(quadrature.size());
1689 project_to_face(reference_cell, quadrature, face_no, points);
1690 return Quadrature<dim>(points, quadrature.get_weights());
1691}
1692
1693
1694
1695template <int dim>
1696Quadrature<dim>
1697QProjector<dim>::project_to_subface(const ReferenceCell &reference_cell,
1698 const SubQuadrature &quadrature,
1699 const unsigned int face_no,
1700 const unsigned int subface_no,
1701 const RefinementCase<dim - 1> &ref_case)
1702{
1703 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1704 ExcNotImplemented());
1705 (void)reference_cell;
1706
1707 std::vector<Point<dim>> points(quadrature.size());
1708 project_to_subface(
1709 reference_cell, quadrature, face_no, subface_no, points, ref_case);
1710 return Quadrature<dim>(points, quadrature.get_weights());
1711}
1712
1713
1714// explicit instantiations; note: we need them all for all dimensions
1715template class QProjector<1>;
1716template class QProjector<2>;
1717template class QProjector<3>;
1718
1719DEAL_II_NAMESPACE_CLOSE
BarycentricPolynomials::get_fe_p_basis
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition polynomials_barycentric.cc:59
Point
Definition point.h:111
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
QProjector::DataSetDescriptor::subface
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
Definition qprojector.cc:1541
QProjector::DataSetDescriptor::face
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition qprojector.cc:1349
QProjector::project_to_all_subfaces
static Quadrature< dim > project_to_all_subfaces(const ReferenceCell &reference_cell, const SubQuadrature &quadrature)
QProjector::project_to_child
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition qprojector.cc:1257
QProjector::project_to_oriented_subface
static Quadrature< dim > project_to_oriented_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const internal::SubfaceCase< dim > ref_case)
Definition qprojector.cc:669
QProjector::project_to_all_faces
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
QProjector::project_to_oriented_face
static Quadrature< dim > project_to_oriented_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation)
Definition qprojector.cc:444
QProjector::project_to_line
static Quadrature< dim > project_to_line(const ReferenceCell &reference_cell, const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
Definition qprojector.cc:1321
QProjector::project_to_subface
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
QProjector::project_to_all_children
static Quadrature< dim > project_to_all_children(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature)
Definition qprojector.cc:1289
QProjector::project_to_face
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
Quadrature::point
const Point< dim > & point(const unsigned int i) const
Quadrature::weight
double weight(const unsigned int i) const
Quadrature::get_weights
const std::vector< double > & get_weights() const
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
Quadrature::size
unsigned int size() const
ReferenceCell::permute_by_combined_orientation
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const unsigned char orientation) const
Definition reference_cell.h:3286
ReferenceCell::get_inverse_combined_orientation
unsigned char get_inverse_combined_orientation(const unsigned char orientation) const
Definition reference_cell.h:3325
ReferenceCell::default_combined_face_orientation
static constexpr unsigned char default_combined_face_orientation()
Definition reference_cell.h:1915
hp::Collection::begin
CollectionIterator< T > begin() const
Definition collection.h:343
hp::Collection::size
unsigned int size() const
Definition collection.h:315
hp::Collection::end
CollectionIterator< T > end() const
Definition collection.h:352
hp::QCollection::max_n_quadrature_points
unsigned int max_n_quadrature_points() const
Definition q_collection.h:173
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
derivative_form.h
geometry_info.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exceptions.h:1938
size
std::size_t size
Definition mpi.cc:734
Polynomials::generate_complete_Lagrange_basis
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition polynomial.cc:686
ReferenceCells::Tetrahedron
constexpr const ReferenceCell Tetrahedron
Definition reference_cell.h:1090
ReferenceCells::Quadrilateral
constexpr const ReferenceCell Quadrilateral
Definition reference_cell.h:1088
ReferenceCells::Wedge
constexpr const ReferenceCell Wedge
Definition reference_cell.h:1094
ReferenceCells::Pyramid
constexpr const ReferenceCell Pyramid
Definition reference_cell.h:1092
ReferenceCells::Triangle
constexpr const ReferenceCell Triangle
Definition reference_cell.h:1086
ReferenceCells::Hexahedron
constexpr const ReferenceCell Hexahedron
Definition reference_cell.h:1096
ReferenceCells::Line
constexpr const ReferenceCell Line
Definition reference_cell.h:1084
internal::split_face_orientation
std::tuple< bool, bool, bool > split_face_orientation(const unsigned char combined_face_orientation)
Definition tria_orientation.h:45
internal::combined_face_orientation
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
Definition tria_orientation.h:32
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6869
polynomials_barycentric.h
reference_cell.h
GeometryInfo::child_cell_on_face
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
determinant
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:2639
tensor_product_polynomials.h
tria_orientation.h
Generated by doxygen 1.9.8