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