Reference documentation for deal.II version GIT relicensing-399-g79d89019c5 2024-04-16 15:00:02+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
20
23
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,
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 {
117 break;
119 xi_scale = 0.5;
120 xi_translation = subface_no % 2 * 0.5;
121 break;
123 eta_scale = 0.5;
124 eta_translation = subface_no % 2 * 0.5;
125 break;
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:
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:
183 }
184 return {};
185 }
186
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]
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 +
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
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 =
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
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;
353 AssertDimension(q_points.size(), 1);
354
355 q_points[0] = Point<dim>(static_cast<double>(face_no));
356}
357
358
359
360template <>
361void
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;
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:
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:
412 }
413 }
414 else
415 {
417 }
418}
419
420
421
422template <>
423void
425 const Quadrature<2> &quadrature,
426 const unsigned int face_no,
427 std::vector<Point<3>> &q_points)
428{
430 (void)reference_cell;
431
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>
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 <>
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{
466
467 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
468 quadrature,
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
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;
492 AssertDimension(q_points.size(), 1);
493
494 q_points[0] = Point<dim>(static_cast<double>(face_no));
495}
496
497
498
499template <>
500void
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;
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:
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:
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:
563 }
564 break;
565 default:
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:
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:
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:
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:
628 }
629 break;
630
631 default:
633 }
634 }
635 else
636 {
638 }
639}
640
641
642
643template <>
644void
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{
653 (void)reference_cell;
654
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>
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,
678{
680 reference_cell,
681 quadrature,
682 face_no,
683 subface_no,
685}
686
687
688
689template <>
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{
702
703 const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
704 quadrature,
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
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 <>
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 <>
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 =
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
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
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 {
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 <>
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(
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 // We will use linear polynomials to map the reference quadrature
915 // points correctly to on faces. There are as many linear shape
916 // functions as there are vertices in the face.
917 const unsigned int n_linear_shape_functions = faces[face_no].size();
918 std::vector<Tensor<1, 2>> shape_derivatives;
919
920 const auto &poly =
921 (n_linear_shape_functions == 3 ?
922 static_cast<const ScalarPolynomialsBase<2> &>(poly_tri) :
923 static_cast<const ScalarPolynomialsBase<2> &>(poly_quad));
924
925 // ... and over all possible orientations
926 for (unsigned char orientation = 0;
927 orientation < reference_cell.n_face_orientations(face_no);
928 ++orientation)
929 {
930 const auto &face = faces[face_no];
931
932 const boost::container::small_vector<Point<3>, 8> support_points =
933 reference_cell.face_reference_cell(face_no)
934 .permute_by_combined_orientation<Point<3>>(face, orientation);
935
936 // the quadrature rule to be projected ...
937 const auto &sub_quadrature_points =
938 quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
939 const auto &sub_quadrature_weights =
940 quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
941
942 // loop over all quadrature points
943 for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
944 {
945 Point<3> mapped_point;
946
947 // map reference quadrature point
948 for (unsigned int i = 0; i < n_linear_shape_functions; ++i)
949 mapped_point +=
950 support_points[i] *
951 poly.compute_value(i, sub_quadrature_points[j]);
952
953 points.push_back(mapped_point);
954
955 // scale quadrature weight
956 const double scaling = [&]() {
957 const unsigned int dim_ = 2;
958 const unsigned int spacedim = 3;
959
961
962 shape_derivatives.resize(n_linear_shape_functions);
963
964 for (unsigned int i = 0; i < n_linear_shape_functions; ++i)
965 shape_derivatives[i] =
966 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
967
968 for (unsigned int k = 0; k < n_linear_shape_functions; ++k)
969 for (unsigned int i = 0; i < spacedim; ++i)
970 for (unsigned int j = 0; j < dim_; ++j)
971 DX_t[j][i] +=
972 shape_derivatives[k][j] * support_points[k][i];
973
975 for (unsigned int i = 0; i < dim_; ++i)
976 for (unsigned int j = 0; j < dim_; ++j)
977 G[i][j] = DX_t[i] * DX_t[j];
978
979 return std::sqrt(determinant(G));
980 }();
981
982 weights.push_back(sub_quadrature_weights[j] * scaling);
983 }
984 }
985 }
986
987 // construct new quadrature rule
988 return Quadrature<3>(std::move(points), std::move(weights));
989 };
990
991 if ((reference_cell == ReferenceCells::Tetrahedron) ||
992 (reference_cell == ReferenceCells::Wedge) ||
993 (reference_cell == ReferenceCells::Pyramid))
994 {
995 std::vector<std::vector<Point<3>>> face_vertex_locations(
996 reference_cell.n_faces());
997 for (const unsigned int f : reference_cell.face_indices())
998 {
999 face_vertex_locations[f].resize(
1000 reference_cell.face_reference_cell(f).n_vertices());
1001 for (const unsigned int v :
1002 reference_cell.face_reference_cell(f).vertex_indices())
1003 face_vertex_locations[f][v] =
1004 reference_cell.face_vertex_location<3>(f, v);
1005 }
1006
1007 return process(face_vertex_locations);
1008 }
1009 else
1010 {
1012
1013 const unsigned int dim = 3;
1014
1015 unsigned int n_points_total = 0;
1016
1017 if (quadrature.size() == 1)
1018 n_points_total =
1019 quadrature[0].size() * GeometryInfo<dim>::faces_per_cell;
1020 else
1021 {
1023 for (const auto &q : quadrature)
1024 n_points_total += q.size();
1025 }
1026
1027 n_points_total *= 8;
1028
1029 // first fix quadrature points
1030 std::vector<Point<dim>> q_points;
1031 q_points.reserve(n_points_total);
1032
1033 std::vector<double> weights;
1034 weights.reserve(n_points_total);
1035
1036 for (unsigned char offset = 0; offset < 8; ++offset)
1037 {
1038 const auto mutation = internal::QProjector::mutate_points_with_offset(
1039 quadrature[0].get_points(), offset);
1040 // project to each face and append results
1041 for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
1042 ++face)
1043 {
1044 const unsigned int q_index = quadrature.size() == 1 ? 0 : face;
1045
1046 internal::QProjector::project_to_hex_face_and_append(
1047 q_index > 0 ? internal::QProjector::mutate_points_with_offset(
1048 quadrature[face].get_points(), offset) :
1049 mutation,
1050 face,
1051 q_points);
1052
1053 std::copy(quadrature[q_index].get_weights().begin(),
1054 quadrature[q_index].get_weights().end(),
1055 std::back_inserter(weights));
1056 }
1057 }
1058
1059 Assert(q_points.size() == n_points_total, ExcInternalError());
1060 Assert(weights.size() == n_points_total, ExcInternalError());
1061
1062 return Quadrature<dim>(std::move(q_points), std::move(weights));
1063 }
1064}
1065
1066
1067
1068template <>
1071 const Quadrature<0> &quadrature)
1072{
1073 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1074 (void)reference_cell;
1075
1076 const unsigned int dim = 1;
1077
1078 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
1079 subfaces_per_face =
1081
1082 // first fix quadrature points
1083 std::vector<Point<dim>> q_points;
1084 q_points.reserve(n_points * n_faces * subfaces_per_face);
1085 std::vector<Point<dim>> help(n_points);
1086
1087 // project to each face and copy
1088 // results
1089 for (unsigned int face = 0; face < n_faces; ++face)
1090 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1091 {
1092 project_to_subface(reference_cell, quadrature, face, subface, help);
1093 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1094 }
1095
1096 // next copy over weights
1097 std::vector<double> weights;
1098 weights.reserve(n_points * n_faces * subfaces_per_face);
1099 for (unsigned int face = 0; face < n_faces; ++face)
1100 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1101 std::copy(quadrature.get_weights().begin(),
1102 quadrature.get_weights().end(),
1103 std::back_inserter(weights));
1104
1105 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1107 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1109
1110 return Quadrature<dim>(std::move(q_points), std::move(weights));
1111}
1112
1113
1114
1115template <>
1118 const SubQuadrature &quadrature)
1119{
1120 if (reference_cell == ReferenceCells::Triangle ||
1121 reference_cell == ReferenceCells::Tetrahedron)
1122 return Quadrature<2>(); // nothing to do
1123
1125
1126 const unsigned int dim = 2;
1127
1128 const unsigned int n_points = quadrature.size(),
1130 subfaces_per_face =
1132
1133 // first fix quadrature points
1134 std::vector<Point<dim>> q_points;
1135 q_points.reserve(n_points * n_faces * subfaces_per_face);
1136 std::vector<Point<dim>> help(n_points);
1137
1138 // project to each face and copy
1139 // results
1140 for (unsigned int face = 0; face < n_faces; ++face)
1141 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1142 {
1143 project_to_subface(reference_cell, quadrature, face, subface, help);
1144 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1145 }
1146
1147 // next copy over weights
1148 std::vector<double> weights;
1149 weights.reserve(n_points * n_faces * subfaces_per_face);
1150 for (unsigned int face = 0; face < n_faces; ++face)
1151 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1152 std::copy(quadrature.get_weights().begin(),
1153 quadrature.get_weights().end(),
1154 std::back_inserter(weights));
1155
1156 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1158 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1160
1161 return Quadrature<dim>(std::move(q_points), std::move(weights));
1162}
1163
1164
1165
1166template <>
1169 const SubQuadrature &quadrature)
1170{
1171 if (reference_cell == ReferenceCells::Triangle ||
1172 reference_cell == ReferenceCells::Tetrahedron)
1173 return Quadrature<3>(); // nothing to do
1174
1176
1177 const unsigned int dim = 3;
1178
1179 const unsigned int n_points = quadrature.size(),
1181 total_subfaces_per_face = 2 + 2 + 4;
1182
1183 // first fix quadrature points
1184 std::vector<Point<dim>> q_points;
1185 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1186
1187 std::vector<double> weights;
1188 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1189
1190 // do the following for all possible mutations of a face (mutation==0
1191 // corresponds to a face with standard orientation, no flip and no rotation)
1192 for (unsigned char offset = 0; offset < 8; ++offset)
1193 {
1194 const auto mutation =
1195 internal::QProjector::mutate_points_with_offset(quadrature.get_points(),
1196 offset);
1197
1198 // project to each face and copy results
1199 for (unsigned int face = 0; face < n_faces; ++face)
1200 for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1201 ref_case >= RefinementCase<dim - 1>::cut_x;
1202 --ref_case)
1203 for (unsigned int subface = 0;
1205 RefinementCase<dim - 1>(ref_case));
1206 ++subface)
1207 {
1208 internal::QProjector::project_to_hex_face_and_append(
1209 mutation,
1210 face,
1211 q_points,
1212 RefinementCase<dim - 1>(ref_case),
1213 subface);
1214
1215 // next copy over weights
1216 std::copy(quadrature.get_weights().begin(),
1217 quadrature.get_weights().end(),
1218 std::back_inserter(weights));
1219 }
1220 }
1221
1222 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1224 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1226
1227 return Quadrature<dim>(std::move(q_points), std::move(weights));
1228}
1229
1230
1231
1232template <int dim>
1235 const Quadrature<dim> &quadrature,
1236 const unsigned int child_no)
1237{
1238 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1240 (void)reference_cell;
1241
1243
1244 const unsigned int n_q_points = quadrature.size();
1245
1246 std::vector<Point<dim>> q_points(n_q_points);
1247 for (unsigned int i = 0; i < n_q_points; ++i)
1248 q_points[i] =
1250 child_no);
1251
1252 // for the weights, things are
1253 // equally simple: copy them and
1254 // scale them
1255 std::vector<double> weights = quadrature.get_weights();
1256 for (unsigned int i = 0; i < n_q_points; ++i)
1257 weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1258
1259 return Quadrature<dim>(q_points, weights);
1260}
1261
1262
1263
1264template <int dim>
1267 const Quadrature<dim> &quadrature)
1268{
1269 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1271 (void)reference_cell;
1272
1273 const unsigned int n_points = quadrature.size(),
1275
1276 std::vector<Point<dim>> q_points(n_points * n_children);
1277 std::vector<double> weights(n_points * n_children);
1278
1279 // project to each child and copy
1280 // results
1281 for (unsigned int child = 0; child < n_children; ++child)
1282 {
1283 Quadrature<dim> help =
1284 project_to_child(reference_cell, quadrature, child);
1285 for (unsigned int i = 0; i < n_points; ++i)
1286 {
1287 q_points[child * n_points + i] = help.point(i);
1288 weights[child * n_points + i] = help.weight(i);
1289 }
1290 }
1291 return Quadrature<dim>(q_points, weights);
1292}
1293
1294
1295
1296template <int dim>
1299 const Quadrature<1> &quadrature,
1300 const Point<dim> &p1,
1301 const Point<dim> &p2)
1302{
1303 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1305 (void)reference_cell;
1306
1307 const unsigned int n = quadrature.size();
1308 std::vector<Point<dim>> points(n);
1309 std::vector<double> weights(n);
1310 const double length = p1.distance(p2);
1311
1312 for (unsigned int k = 0; k < n; ++k)
1313 {
1314 const double alpha = quadrature.point(k)[0];
1315 points[k] = alpha * p2;
1316 points[k] += (1. - alpha) * p1;
1317 weights[k] = length * quadrature.weight(k);
1318 }
1319 return Quadrature<dim>(points, weights);
1320}
1321
1322
1323
1324template <int dim>
1327 const unsigned int face_no,
1328 const bool face_orientation,
1329 const bool face_flip,
1330 const bool face_rotation,
1331 const unsigned int n_quadrature_points)
1332{
1333 return face(reference_cell,
1334 face_no,
1335 internal::combined_face_orientation(face_orientation,
1336 face_rotation,
1337 face_flip),
1338 n_quadrature_points);
1339}
1340
1341
1342
1343template <int dim>
1346 const ReferenceCell &reference_cell,
1347 const unsigned int face_no,
1348 const unsigned char combined_orientation,
1349 const unsigned int n_quadrature_points)
1350{
1351 // TODO: once we move to representing the default orientation as 0 (instead of
1352 // 1) we can get rid of the dim = 1 check
1353 Assert(dim == 1 ||
1354 (combined_orientation < reference_cell.n_face_orientations(face_no)),
1356 if (reference_cell == ReferenceCells::Triangle ||
1357 reference_cell == ReferenceCells::Tetrahedron)
1358 {
1359 if (dim == 2)
1360 return {(2 * face_no +
1361 (combined_orientation ==
1363 1 :
1364 0)) *
1365 n_quadrature_points};
1366 else if (dim == 3)
1367 {
1368 return {(reference_cell.n_face_orientations(face_no) * face_no +
1369 combined_orientation) *
1370 n_quadrature_points};
1371 }
1372 }
1373
1374 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1376
1378
1379 switch (dim)
1380 {
1381 case 1:
1382 case 2:
1383 return face_no * n_quadrature_points;
1384 case 3:
1385 return (face_no +
1386 GeometryInfo<dim>::faces_per_cell * combined_orientation) *
1387 n_quadrature_points;
1388 default:
1390 }
1392}
1393
1394
1395
1396template <int dim>
1399 const ReferenceCell &reference_cell,
1400 const unsigned int face_no,
1401 const bool face_orientation,
1402 const bool face_flip,
1403 const bool face_rotation,
1404 const hp::QCollection<dim - 1> &quadrature)
1405{
1406 return face(reference_cell,
1407 face_no,
1408 internal::combined_face_orientation(face_orientation,
1409 face_rotation,
1410 face_flip),
1411 quadrature);
1412}
1413
1414
1415
1416template <int dim>
1419 const ReferenceCell &reference_cell,
1420 const unsigned int face_no,
1421 const unsigned char combined_orientation,
1422 const hp::QCollection<dim - 1> &quadrature)
1423{
1424 if (reference_cell == ReferenceCells::Triangle ||
1425 reference_cell == ReferenceCells::Tetrahedron ||
1426 reference_cell == ReferenceCells::Wedge ||
1427 reference_cell == ReferenceCells::Pyramid)
1428 {
1429 unsigned int offset = 0;
1430 Assert(combined_orientation < reference_cell.n_face_orientations(face_no),
1432
1433 static const unsigned int X = numbers::invalid_unsigned_int;
1434 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1435 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1436 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1437 static const std::array<unsigned int, 5> scale_pyramid = {
1438 {8, 6, 6, 6, 6}};
1439
1440 const auto &scale =
1441 (reference_cell == ReferenceCells::Triangle) ?
1442 scale_tri :
1443 ((reference_cell == ReferenceCells::Tetrahedron) ?
1444 scale_tet :
1445 ((reference_cell == ReferenceCells::Wedge) ? scale_wedge :
1446 scale_pyramid));
1447
1448 if (quadrature.size() == 1)
1449 offset = scale[0] * quadrature[0].size() * face_no;
1450 else
1451 for (unsigned int i = 0; i < face_no; ++i)
1452 offset += scale[i] * quadrature[i].size();
1453
1454 if (dim == 2)
1455 return {offset +
1456 (combined_orientation ==
1458 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1459 else if (dim == 3)
1460 {
1461 return {offset +
1462 combined_orientation *
1463 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1464 }
1465 }
1466
1467 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1469
1471
1472 switch (dim)
1473 {
1474 case 1:
1475 case 2:
1476 {
1477 if (quadrature.size() == 1)
1478 return quadrature[0].size() * face_no;
1479 else
1480 {
1481 unsigned int result = 0;
1482 for (unsigned int i = 0; i < face_no; ++i)
1483 result += quadrature[i].size();
1484 return result;
1485 }
1486 }
1487 case 3:
1488 {
1489 if (quadrature.size() == 1)
1490 return (face_no +
1491 combined_orientation * GeometryInfo<dim>::faces_per_cell) *
1492 quadrature[0].size();
1493 else
1494 {
1495 unsigned int n_points_i = 0;
1496 for (unsigned int i = 0; i < face_no; ++i)
1497 n_points_i += quadrature[i].size();
1498
1499 unsigned int n_points = 0;
1500 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1501 ++i)
1502 n_points += quadrature[i].size();
1503
1504 return n_points_i + combined_orientation * n_points;
1505 }
1506 }
1507
1508 default:
1510 }
1512}
1513
1514
1515
1516template <int dim>
1519 const ReferenceCell &reference_cell,
1520 const unsigned int face_no,
1521 const unsigned int subface_no,
1522 const bool face_orientation,
1523 const bool face_flip,
1524 const bool face_rotation,
1525 const unsigned int n_quadrature_points,
1526 const internal::SubfaceCase<dim> ref_case)
1527{
1529 reference_cell,
1530 face_no,
1531 subface_no,
1532 internal::combined_face_orientation(face_orientation,
1533 face_rotation,
1534 face_flip),
1535 n_quadrature_points,
1536 ref_case);
1537}
1538
1539
1540
1541template <>
1544 const ReferenceCell &reference_cell,
1545 const unsigned int face_no,
1546 const unsigned int subface_no,
1547 const unsigned char /*combined_orientation*/,
1548 const unsigned int n_quadrature_points,
1550{
1551 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1552 (void)reference_cell;
1553
1557
1558 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1559 n_quadrature_points);
1560}
1561
1562
1563
1564template <>
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,
1573{
1575 (void)reference_cell;
1576
1580
1581 return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1582 n_quadrature_points);
1583}
1584
1585
1586
1587template <>
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<3> ref_case)
1596{
1597 const unsigned int dim = 3;
1598
1599 const auto [face_orientation, face_rotation, face_flip] =
1600 internal::split_face_orientation(combined_orientation);
1601
1603 (void)reference_cell;
1604
1608
1609 // in 3d, we have to account for faces that
1610 // have non-standard orientation. thus, we
1611 // have to store _eight_ data sets per face
1612 // or subface already for the isotropic
1613 // case. Additionally, we have three
1614 // different refinement cases, resulting in
1615 // <tt>4 + 2 + 2 = 8</tt> different subfaces
1616 // for each face.
1617 const unsigned int total_subfaces_per_face = 8;
1618
1619 // set up a table with the offsets for a
1620 // given refinement case respecting the
1621 // corresponding number of subfaces. the
1622 // index corresponds to (RefineCase::Type - 1)
1623
1624 // note, that normally we should use the
1625 // obvious offsets 0,2,6. However, prior to
1626 // the implementation of anisotropic
1627 // refinement, in many places of the library
1628 // the convention was used, that the first
1629 // dataset with offset 0 corresponds to a
1630 // standard (isotropic) face
1631 // refinement. therefore we use the offsets
1632 // 6,4,0 here to stick to that (implicit)
1633 // convention
1634 static const unsigned int ref_case_offset[3] = {
1635 6, // cut_x
1636 4, // cut_y
1637 0 // cut_xy
1638 };
1639
1640 const std::pair<unsigned int, RefinementCase<2>>
1641 final_subface_no_and_ref_case =
1642 internal::QProjector::select_subface_no_and_refinement_case(
1643 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1644
1645 return ((face_no * total_subfaces_per_face +
1646 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1647 final_subface_no_and_ref_case.first) +
1648 GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face *
1649 combined_orientation) *
1650 n_quadrature_points;
1651}
1652
1653
1654
1655template <int dim>
1658 const SubQuadrature &quadrature,
1659 const unsigned int face_no)
1660{
1661 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1663 (void)reference_cell;
1664
1665 std::vector<Point<dim>> points(quadrature.size());
1666 project_to_face(reference_cell, quadrature, face_no, points);
1667 return Quadrature<dim>(points, quadrature.get_weights());
1668}
1669
1670
1671
1672template <int dim>
1675 const SubQuadrature &quadrature,
1676 const unsigned int face_no,
1677 const unsigned int subface_no,
1678 const RefinementCase<dim - 1> &ref_case)
1679{
1680 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1682 (void)reference_cell;
1683
1684 std::vector<Point<dim>> points(quadrature.size());
1686 reference_cell, quadrature, face_no, subface_no, points, ref_case);
1687 return Quadrature<dim>(points, quadrature.get_weights());
1688}
1689
1690
1691// explicit instantiations; note: we need them all for all dimensions
1692template class QProjector<1>;
1693template class QProjector<2>;
1694template class QProjector<3>;
1695
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
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)
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)
static Quadrature< dim > project_to_all_subfaces(const ReferenceCell &reference_cell, const SubQuadrature &quadrature)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
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)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
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)
static Quadrature< dim > project_to_line(const ReferenceCell &reference_cell, const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
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)
static Quadrature< dim > project_to_all_children(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const unsigned char orientation) const
static constexpr unsigned char default_combined_face_orientation()
CollectionIterator< T > begin() const
Definition collection.h:283
unsigned int size() const
Definition collection.h:264
CollectionIterator< T > end() const
Definition collection.h:292
unsigned int max_n_quadrature_points() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_ASSERT_UNREACHABLE()
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Line
std::tuple< bool, bool, bool > split_face_orientation(const unsigned char combined_face_orientation)
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
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)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)