deal.II version GIT relicensing-3509-g790ffced1c 2025-06-14 07:20: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
18
21
22#include <boost/container/small_vector.hpp>
23
24
26
27
28namespace internal
29{
30 namespace QProjector
31 {
32 namespace
33 {
34 // Reflect points across the y = x line.
35 std::vector<Point<2>>
36 reflect(const std::vector<Point<2>> &points)
37 {
38 std::vector<Point<2>> q_points;
39 q_points.reserve(points.size());
40 for (const Point<2> &p : points)
41 q_points.emplace_back(p[1], p[0]);
42
43 return q_points;
44 }
45
46
47 // Rotate points in the plane around the positive z axis (i.e.,
48 // counter-clockwise).
49 std::vector<Point<2>>
50 rotate(const std::vector<Point<2>> &points, const unsigned int n_times)
51 {
52 std::vector<Point<2>> q_points;
53 q_points.reserve(points.size());
54 switch (n_times % 4)
55 {
56 case 0:
57 // 0 degree. the point remains as it is.
58 for (const Point<2> &p : points)
59 q_points.push_back(p);
60 break;
61 case 1:
62 // 90 degree counterclockwise
63 for (const Point<2> &p : points)
64 q_points.emplace_back(1.0 - p[1], p[0]);
65 break;
66 case 2:
67 // 180 degree counterclockwise
68 for (const Point<2> &p : points)
69 q_points.emplace_back(1.0 - p[0], 1.0 - p[1]);
70 break;
71 case 3:
72 // 270 degree counterclockwise
73 for (const Point<2> &p : points)
74 q_points.emplace_back(p[1], 1.0 - p[0]);
75 break;
76 }
77
78 return q_points;
79 }
80
87 void
88 project_to_hex_face_and_append(
89 const std::vector<Point<2>> &points,
90 const unsigned int face_no,
91 std::vector<Point<3>> &q_points,
93 const unsigned int subface_no = 0)
94 {
95 // one coordinate is at a const value. for faces 0, 2 and 4 this value
96 // is 0.0, for faces 1, 3 and 5 it is 1.0
97 const double const_value = face_no % 2;
98
99 // local 2d coordinates are xi and eta, global 3d coordinates are x, y
100 // and z. those have to be mapped. the following indices tell, which
101 // global coordinate (0->x, 1->y, 2->z) corresponds to which local one
102 const unsigned int xi_index = (1 + face_no / 2) % 3,
103 eta_index = (2 + face_no / 2) % 3,
104 const_index = face_no / 2;
105
106 // for a standard face (no refinement), we use the default values of
107 // the xi and eta scales and translations, otherwise the xi and eta
108 // values will be scaled (by factor 0.5 or factor 1.0) depending on
109 // the refinement case and translated (by 0.0 or 0.5) depending on the
110 // refinement case and subface_no
111 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
112 eta_translation = 0.0;
113
114 // set the scale and translation parameter for individual subfaces
115 switch (ref_case)
116 {
118 break;
120 xi_scale = 0.5;
121 xi_translation = subface_no % 2 * 0.5;
122 break;
124 eta_scale = 0.5;
125 eta_translation = subface_no % 2 * 0.5;
126 break;
128 xi_scale = 0.5;
129 eta_scale = 0.5;
130 xi_translation = int(subface_no % 2) * 0.5;
131 eta_translation = int(subface_no / 2) * 0.5;
132 break;
133 default:
135 break;
136 }
137
138 // finally, compute the scaled, translated, projected quadrature
139 // points
140 for (const Point<2> &p : points)
141 {
142 Point<3> cell_point;
143 cell_point[xi_index] = xi_scale * p[0] + xi_translation;
144 cell_point[eta_index] = eta_scale * p[1] + eta_translation;
145 cell_point[const_index] = const_value;
146 q_points.push_back(cell_point);
147 }
148 }
149
150 std::vector<Point<2>>
151 mutate_points_with_offset(
152 const std::vector<Point<2>> &points,
153 const types::geometric_orientation combined_orientation)
154 {
155 // These rotations are backwards (relative to the standard notion of,
156 // e.g., what rotation index 7 means) since they are rotations about the
157 // positive z axis in 2d: i.e., they are done from the perspective of
158 // 'inside' a cell instead of the perspective of an abutting cell.
159 //
160 // For example: consider points on face 4 of a hexahedron with
161 // orientation 3. In 2d, rotating such points clockwise is the same as
162 // rotating them counter-clockwise from the perspective of the abutting
163 // face. Hence, such points must be rotated 90 degrees
164 // counter-clockwise.
165 switch (combined_orientation)
166 {
167 case 1:
168 return reflect(points);
169 case 3:
170 return rotate(reflect(points), 3);
171 case 5:
172 return rotate(reflect(points), 2);
173 case 7:
174 return rotate(reflect(points), 1);
175 case 0:
176 return points;
177 case 2:
178 return rotate(points, 1);
179 case 4:
180 return rotate(points, 2);
181 case 6:
182 return rotate(points, 3);
183 default:
185 }
186 return {};
187 }
188
189
190
191 std::pair<unsigned int, RefinementCase<2>>
192 select_subface_no_and_refinement_case(
193 const unsigned int subface_no,
194 const bool face_orientation,
195 const bool face_flip,
196 const bool face_rotation,
197 const internal::SubfaceCase<3> ref_case)
198 {
199 constexpr int dim = 3;
200 // for each subface of a given FaceRefineCase
201 // there is a corresponding equivalent
202 // subface number of one of the "standard"
203 // RefineCases (cut_x, cut_y, cut_xy). Map
204 // the given values to those equivalent
205 // ones.
206
207 // first, define an invalid number
208 static const unsigned int e = numbers::invalid_unsigned_int;
209
210 static const RefinementCase<dim - 1>
211 equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
213 // case_none. there should be only
214 // invalid values here. However, as
215 // this function is also called (in
216 // tests) for cells which have no
217 // refined faces, use isotropic
218 // refinement instead
219 {RefinementCase<dim - 1>::cut_xy,
220 RefinementCase<dim - 1>::cut_xy,
221 RefinementCase<dim - 1>::cut_xy,
222 RefinementCase<dim - 1>::cut_xy},
223 // case_x
224 {RefinementCase<dim - 1>::cut_x,
225 RefinementCase<dim - 1>::cut_x,
226 RefinementCase<dim - 1>::no_refinement,
227 RefinementCase<dim - 1>::no_refinement},
228 // case_x1y
229 {RefinementCase<dim - 1>::cut_xy,
230 RefinementCase<dim - 1>::cut_xy,
231 RefinementCase<dim - 1>::cut_x,
232 RefinementCase<dim - 1>::no_refinement},
233 // case_x2y
234 {RefinementCase<dim - 1>::cut_x,
235 RefinementCase<dim - 1>::cut_xy,
236 RefinementCase<dim - 1>::cut_xy,
237 RefinementCase<dim - 1>::no_refinement},
238 // case_x1y2y
239 {RefinementCase<dim - 1>::cut_xy,
240 RefinementCase<dim - 1>::cut_xy,
241 RefinementCase<dim - 1>::cut_xy,
242 RefinementCase<dim - 1>::cut_xy},
243 // case_y
244 {RefinementCase<dim - 1>::cut_y,
245 RefinementCase<dim - 1>::cut_y,
246 RefinementCase<dim - 1>::no_refinement,
247 RefinementCase<dim - 1>::no_refinement},
248 // case_y1x
249 {RefinementCase<dim - 1>::cut_xy,
250 RefinementCase<dim - 1>::cut_xy,
251 RefinementCase<dim - 1>::cut_y,
252 RefinementCase<dim - 1>::no_refinement},
253 // case_y2x
254 {RefinementCase<dim - 1>::cut_y,
255 RefinementCase<dim - 1>::cut_xy,
256 RefinementCase<dim - 1>::cut_xy,
257 RefinementCase<dim - 1>::no_refinement},
258 // case_y1x2x
259 {RefinementCase<dim - 1>::cut_xy,
260 RefinementCase<dim - 1>::cut_xy,
261 RefinementCase<dim - 1>::cut_xy,
262 RefinementCase<dim - 1>::cut_xy},
263 // case_xy (case_isotropic)
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
269 static const unsigned int
270 equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic +
272 {// case_none, see above
273 {0, 1, 2, 3},
274 // case_x
275 {0, 1, e, e},
276 // case_x1y
277 {0, 2, 1, e},
278 // case_x2y
279 {0, 1, 3, e},
280 // case_x1y2y
281 {0, 2, 1, 3},
282 // case_y
283 {0, 1, e, e},
284 // case_y1x
285 {0, 1, 1, e},
286 // case_y2x
287 {0, 2, 3, e},
288 // case_y1x2x
289 {0, 1, 2, 3},
290 // case_xy (case_isotropic)
291 {0, 1, 2, 3}};
292
293 // If face-orientation or face_rotation are
294 // non-standard, cut_x and cut_y have to be
295 // exchanged.
296 static const RefinementCase<dim - 1> ref_case_permutation[4] = {
297 RefinementCase<dim - 1>::no_refinement,
298 RefinementCase<dim - 1>::cut_y,
299 RefinementCase<dim - 1>::cut_x,
300 RefinementCase<dim - 1>::cut_xy};
301
302 // set a corresponding (equivalent)
303 // RefineCase and subface number
304 const RefinementCase<dim - 1> equ_ref_case =
305 equivalent_refine_case[ref_case][subface_no];
306 const unsigned int equ_subface_no =
307 equivalent_subface_number[ref_case][subface_no];
308 // make sure, that we got a valid subface and RefineCase
311 Assert(equ_subface_no != e, ExcInternalError());
312 // now, finally respect non-standard faces
313 const RefinementCase<dim - 1> final_ref_case =
314 (face_orientation == face_rotation ?
315 ref_case_permutation[equ_ref_case] :
316 equ_ref_case);
317
318 const unsigned int final_subface_no =
320 final_ref_case),
321 4,
322 equ_subface_no,
323 face_orientation,
324 face_flip,
325 face_rotation,
326 equ_ref_case);
327
328 return std::make_pair(final_subface_no, final_ref_case);
329 }
330
357 template <int dim>
358 void
359 append_subobject_rule(
360 const ReferenceCell &face_reference_cell,
361 const Quadrature<dim - 1> &quadrature,
362 const std::vector<Point<dim>> &vertices,
363 const double measure,
364 const types::geometric_orientation combined_orientation,
365 std::vector<Point<dim>> &points,
366 std::vector<double> &weights)
367 {
368 const auto support_points =
369 face_reference_cell.permute_by_combined_orientation(
370 make_const_array_view(vertices),
371 face_reference_cell.get_inverse_combined_orientation(
372 combined_orientation));
373
374 for (unsigned int j = 0; j < quadrature.size(); ++j)
375 {
376 Point<dim> mapped_point;
377
378 // map reference quadrature point
379 for (const unsigned int vertex_no :
380 face_reference_cell.vertex_indices())
381 mapped_point +=
382 support_points[vertex_no] *
383 face_reference_cell.d_linear_shape_function(quadrature.point(j),
384 vertex_no);
385
386 points.push_back(mapped_point);
387
388 // rescale quadrature weights so that the sum of the weights on
389 // each face equals the measure of that face.
390 weights.push_back(quadrature.weight(j) * measure /
391 face_reference_cell.volume());
392 }
393 }
394 } // namespace
395 } // namespace QProjector
396} // namespace internal
397
398
399
400template <>
401void
403 const Quadrature<0> &quadrature,
404 const unsigned int face_no,
405 std::vector<Point<1>> &q_points)
406{
407 AssertDimension(quadrature.size(), q_points.size());
408 const auto face_quadrature =
409 QProjector<1>::project_to_face(reference_cell,
410 quadrature,
411 face_no,
413 q_points = face_quadrature.get_points();
414}
415
416
417
418template <>
419void
421 const Quadrature<1> &quadrature,
422 const unsigned int face_no,
423 std::vector<Point<2>> &q_points)
424{
425 AssertDimension(quadrature.size(), q_points.size());
426 const auto face_quadrature =
427 QProjector<2>::project_to_face(reference_cell,
428 quadrature,
429 face_no,
431 q_points = face_quadrature.get_points();
432}
433
434
435
436template <>
437void
439 const Quadrature<2> &quadrature,
440 const unsigned int face_no,
441 std::vector<Point<3>> &q_points)
442{
444 (void)reference_cell;
445
447 Assert(q_points.size() == quadrature.size(),
448 ExcDimensionMismatch(q_points.size(), quadrature.size()));
449 q_points.clear();
450 internal::QProjector::project_to_hex_face_and_append(quadrature.get_points(),
451 face_no,
452 q_points);
453}
454
455
456
457template <int dim>
460 const Quadrature<dim - 1> &quadrature,
461 const unsigned int face_no,
462 const bool face_orientation,
463 const bool face_flip,
464 const bool face_rotation)
465{
467 reference_cell,
468 quadrature,
469 face_no,
471 face_rotation,
472 face_flip));
473}
474
475
476
477template <int dim>
480 const ReferenceCell &reference_cell,
481 const Quadrature<dim - 1> &quadrature,
482 const unsigned int face_no,
483 const types::geometric_orientation combined_orientation)
484{
485 AssertIndexRange(face_no, reference_cell.n_faces());
486 AssertIndexRange(combined_orientation,
487 reference_cell.n_face_orientations(face_no));
488 AssertDimension(reference_cell.get_dimension(), dim);
489
490 std::vector<Point<dim>> points;
491 std::vector<double> weights;
492
493 const ReferenceCell face_reference_cell =
494 reference_cell.face_reference_cell(face_no);
495 std::vector<Point<dim>> face_vertices(face_reference_cell.n_vertices());
496 for (const unsigned int vertex_no : face_reference_cell.vertex_indices())
497 face_vertices[vertex_no] =
498 reference_cell.face_vertex_location<dim>(face_no, vertex_no);
499 internal::QProjector::append_subobject_rule(face_reference_cell,
500 quadrature,
501 face_vertices,
502 reference_cell.face_measure(
503 face_no),
504 combined_orientation,
505 points,
506 weights);
507
508 return Quadrature<dim>(std::move(points), std::move(weights));
509}
510
511
512
513template <>
514void
516 const Quadrature<0> &,
517 const unsigned int face_no,
518 const unsigned int,
519 std::vector<Point<1>> &q_points,
520 const RefinementCase<0> &)
521{
522 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
523 (void)reference_cell;
524
525 const unsigned int dim = 1;
527 AssertDimension(q_points.size(), 1);
528
529 q_points[0] = Point<dim>(static_cast<double>(face_no));
530}
531
532
533
534template <>
535void
537 const Quadrature<1> &quadrature,
538 const unsigned int face_no,
539 const unsigned int subface_no,
540 std::vector<Point<2>> &q_points,
541 const RefinementCase<1> &ref_case)
542{
543 AssertDimension(quadrature.size(), q_points.size());
544 const auto face_quadrature =
545 project_to_subface(reference_cell,
546 quadrature,
547 face_no,
548 subface_no,
550 ref_case);
551 q_points = face_quadrature.get_points();
552}
553
554
555
556template <>
557void
559 const Quadrature<2> &quadrature,
560 const unsigned int face_no,
561 const unsigned int subface_no,
562 std::vector<Point<3>> &q_points,
563 const RefinementCase<2> &ref_case)
564{
566 (void)reference_cell;
567 AssertDimension(quadrature.size(), q_points.size());
568 const auto face_quadrature =
569 project_to_subface(reference_cell,
570 quadrature,
571 face_no,
572 subface_no,
574 ref_case);
575 q_points = face_quadrature.get_points();
576}
577
578
579
580template <int dim>
583 const ReferenceCell &reference_cell,
584 const Quadrature<dim - 1> &quadrature,
585 const unsigned int face_no,
586 const unsigned int subface_no,
587 const bool,
588 const bool,
589 const bool,
591{
593 reference_cell,
594 quadrature,
595 face_no,
596 subface_no,
598}
599
600
601
602template <int dim>
605 const ReferenceCell &reference_cell,
606 const SubQuadrature &quadrature,
607 const unsigned int face_no,
608 const unsigned int subface_no,
609 const types::geometric_orientation combined_orientation,
610 const RefinementCase<dim - 1> &ref_case)
611{
612 AssertIndexRange(face_no, reference_cell.n_faces());
613 AssertIndexRange(combined_orientation,
614 reference_cell.n_face_orientations(face_no));
615 AssertDimension(reference_cell.get_dimension(), dim);
616 AssertIndexRange(subface_no,
617 reference_cell.face_reference_cell(face_no)
618 .template n_children<dim - 1>(ref_case));
619
620 std::vector<Point<dim>> q_points;
621 std::vector<double> q_weights = quadrature.get_weights();
622 q_points.reserve(quadrature.size());
623
624 if constexpr (dim == 1)
625 {
626 AssertDimension(quadrature.size(), 1);
627 q_points.emplace_back(static_cast<double>(face_no));
628 }
629 else if constexpr (dim == 2)
630 {
631 if (reference_cell == ReferenceCells::Triangle)
632 // use linear polynomial to map the reference quadrature points
633 // correctly on faces, i.e., BarycentricPolynomials<1>(1)
634 for (unsigned int p = 0; p < quadrature.size(); ++p)
635 {
636 if (face_no == 0)
637 {
638 if (subface_no == 0)
639 q_points.emplace_back(quadrature.point(p)[0] / 2, 0);
640 else
641 q_points.emplace_back(0.5 + quadrature.point(p)[0] / 2, 0);
642 }
643 else if (face_no == 1)
644 {
645 if (subface_no == 0)
646 q_points.emplace_back(1 - quadrature.point(p)[0] / 2,
647 quadrature.point(p)[0] / 2);
648 else
649 q_points.emplace_back(0.5 - quadrature.point(p)[0] / 2,
650 0.5 + quadrature.point(p)[0] / 2);
651 }
652 else if (face_no == 2)
653 {
654 if (subface_no == 0)
655 q_points.emplace_back(0, 1 - quadrature.point(p)[0] / 2);
656 else
657 q_points.emplace_back(0, 0.5 - quadrature.point(p)[0] / 2);
658 }
659 else
661 }
662 else if (reference_cell == ReferenceCells::Quadrilateral)
663 for (unsigned int p = 0; p < quadrature.size(); ++p)
664 {
665 if (face_no == 0)
666 {
667 if (subface_no == 0)
668 q_points.emplace_back(0, quadrature.point(p)[0] / 2);
669 else
670 q_points.emplace_back(0, quadrature.point(p)[0] / 2 + 0.5);
671 }
672 else if (face_no == 1)
673 {
674 if (subface_no == 0)
675 q_points.emplace_back(1, quadrature.point(p)[0] / 2);
676 else
677 q_points.emplace_back(1, quadrature.point(p)[0] / 2 + 0.5);
678 }
679 else if (face_no == 2)
680 {
681 if (subface_no == 0)
682 q_points.emplace_back(quadrature.point(p)[0] / 2, 0);
683 else
684 q_points.emplace_back(quadrature.point(p)[0] / 2 + 0.5, 0);
685 }
686 else if (face_no == 3)
687 {
688 if (subface_no == 0)
689 q_points.emplace_back(quadrature.point(p)[0] / 2, 1);
690 else
691 q_points.emplace_back(quadrature.point(p)[0] / 2 + 0.5, 1);
692 }
693 else
695 }
696 else
698
699 if (combined_orientation == numbers::reverse_line_orientation)
700 {
701 std::reverse(q_points.begin(), q_points.end());
702 std::reverse(q_weights.begin(), q_weights.end());
703 }
704 for (auto &w : q_weights)
705 w *= reference_cell.face_measure(face_no);
706 }
707 else if constexpr (dim == 3)
708 {
710 internal::QProjector::project_to_hex_face_and_append(
711 quadrature.get_points(), face_no, q_points, ref_case, subface_no);
712 }
713 else
714 {
716 }
717
718 return Quadrature<dim>(std::move(q_points), std::move(q_weights));
719}
720
721
722
723template <int dim>
726 const ReferenceCell &reference_cell,
727 const hp::QCollection<dim - 1> &quadrature)
728{
729 std::vector<Point<dim>> points;
730 std::vector<double> weights;
731
732 for (const unsigned int face_no : reference_cell.face_indices())
733 {
734 const ReferenceCell face_reference_cell =
735 reference_cell.face_reference_cell(face_no);
736 std::vector<Point<dim>> face_vertices(face_reference_cell.n_vertices());
737 for (const unsigned int vertex_no : face_reference_cell.vertex_indices())
738 face_vertices[vertex_no] =
739 reference_cell.face_vertex_location<dim>(face_no, vertex_no);
740
741 for (types::geometric_orientation combined_orientation = 0;
742 combined_orientation < reference_cell.n_face_orientations(face_no);
743 ++combined_orientation)
744 internal::QProjector::append_subobject_rule(
745 face_reference_cell,
746 quadrature[quadrature.size() == 1 ? 0 : face_no],
747 face_vertices,
748 reference_cell.face_measure(face_no),
749 combined_orientation,
750 points,
751 weights);
752 }
753
754 return Quadrature<dim>(std::move(points), std::move(weights));
755}
756
757
758
759template <>
762 const Quadrature<0> &quadrature)
763{
764 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
765 (void)reference_cell;
766
767 const unsigned int dim = 1;
768
769 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
770 subfaces_per_face =
772
773 // first fix quadrature points
774 std::vector<Point<dim>> q_points;
775 q_points.reserve(n_points * n_faces * subfaces_per_face);
776 std::vector<Point<dim>> help(n_points);
777
778 // project to each face and copy
779 // results
780 for (unsigned int face = 0; face < n_faces; ++face)
781 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
782 {
783 project_to_subface(reference_cell, quadrature, face, subface, help);
784 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
785 }
786
787 // next copy over weights
788 std::vector<double> weights;
789 weights.reserve(n_points * n_faces * subfaces_per_face);
790 for (unsigned int face = 0; face < n_faces; ++face)
791 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
792 std::copy(quadrature.get_weights().begin(),
793 quadrature.get_weights().end(),
794 std::back_inserter(weights));
795
796 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
798 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
800
801 return Quadrature<dim>(std::move(q_points), std::move(weights));
802}
803
804
805
806template <>
809 const SubQuadrature &quadrature)
810{
811 Assert(reference_cell == ReferenceCells::Quadrilateral ||
812 reference_cell == ReferenceCells::Triangle,
814
815 const unsigned int dim = 2;
816
817 std::vector<Point<dim>> q_points;
818 std::vector<double> weights;
819
820 // project to each face and copy
821 // results
822 for (unsigned int face = 0; face < reference_cell.n_faces(); ++face)
823 for (types::geometric_orientation orientation = 0;
824 orientation < reference_cell.n_face_orientations(face);
825 ++orientation)
826 for (unsigned int subface = 0;
827 subface <
828 reference_cell.face_reference_cell(face).n_isotropic_children();
829 ++subface)
830 {
831 const unsigned int local_subface =
832 orientation == numbers::reverse_line_orientation ? 1 - subface :
833 subface;
834 const auto sub_quadrature =
835 project_to_subface(reference_cell,
836 quadrature,
837 face,
838 local_subface,
839 orientation,
841 q_points.insert(q_points.end(),
842 sub_quadrature.get_points().begin(),
843 sub_quadrature.get_points().end());
844 weights.insert(weights.end(),
845 sub_quadrature.get_weights().begin(),
846 sub_quadrature.get_weights().end());
847 }
848
849 return Quadrature<dim>(std::move(q_points), std::move(weights));
850}
851
852
853
854template <>
857 const SubQuadrature &quadrature)
858{
859 if (reference_cell == ReferenceCells::Triangle ||
860 reference_cell == ReferenceCells::Tetrahedron)
861 return Quadrature<3>(); // nothing to do
862
864
865 const unsigned int dim = 3;
866
867 const unsigned int n_points = quadrature.size(),
869 total_subfaces_per_face = 2 + 2 + 4;
870
871 // first fix quadrature points
872 std::vector<Point<dim>> q_points;
873 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
874
875 std::vector<double> weights;
876 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
877
878 // do the following for all possible mutations of a face (mutation==0
879 // corresponds to a face with standard orientation, no flip and no rotation)
880 for (unsigned char offset = 0; offset < 8; ++offset)
881 {
882 const auto mutation =
883 internal::QProjector::mutate_points_with_offset(quadrature.get_points(),
884 offset);
885
886 // project to each face and copy results
887 for (unsigned int face = 0; face < n_faces; ++face)
888 for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
889 ref_case >= RefinementCase<dim - 1>::cut_x;
890 --ref_case)
891 for (unsigned int subface = 0;
893 RefinementCase<dim - 1>(ref_case));
894 ++subface)
895 {
896 internal::QProjector::project_to_hex_face_and_append(
897 mutation,
898 face,
899 q_points,
900 RefinementCase<dim - 1>(ref_case),
901 subface);
902
903 // next copy over weights
904 std::copy(quadrature.get_weights().begin(),
905 quadrature.get_weights().end(),
906 std::back_inserter(weights));
907 }
908 }
909
910 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
912 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
914
915 return Quadrature<dim>(std::move(q_points), std::move(weights));
916}
917
918
919
920template <int dim>
923 const Quadrature<dim> &quadrature,
924 const unsigned int child_no)
925{
926 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
928 (void)reference_cell;
929
931
932 const unsigned int n_q_points = quadrature.size();
933
934 std::vector<Point<dim>> q_points(n_q_points);
935 for (unsigned int i = 0; i < n_q_points; ++i)
936 q_points[i] =
938 child_no);
939
940 // for the weights, things are
941 // equally simple: copy them and
942 // scale them
943 std::vector<double> weights = quadrature.get_weights();
944 for (unsigned int i = 0; i < n_q_points; ++i)
945 weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
946
947 return Quadrature<dim>(q_points, weights);
948}
949
950
951
952template <int dim>
955 const Quadrature<dim> &quadrature)
956{
957 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
959 (void)reference_cell;
960
961 const unsigned int n_points = quadrature.size(),
963
964 std::vector<Point<dim>> q_points(n_points * n_children);
965 std::vector<double> weights(n_points * n_children);
966
967 // project to each child and copy
968 // results
969 for (unsigned int child = 0; child < n_children; ++child)
970 {
971 Quadrature<dim> help =
972 project_to_child(reference_cell, quadrature, child);
973 for (unsigned int i = 0; i < n_points; ++i)
974 {
975 q_points[child * n_points + i] = help.point(i);
976 weights[child * n_points + i] = help.weight(i);
977 }
978 }
979 return Quadrature<dim>(q_points, weights);
980}
981
982
983
984template <int dim>
987 const Quadrature<1> &quadrature,
988 const Point<dim> &p1,
989 const Point<dim> &p2)
990{
991 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
993 (void)reference_cell;
994
995 const unsigned int n = quadrature.size();
996 std::vector<Point<dim>> points(n);
997 std::vector<double> weights(n);
998 const double length = p1.distance(p2);
999
1000 for (unsigned int k = 0; k < n; ++k)
1001 {
1002 const double alpha = quadrature.point(k)[0];
1003 points[k] = alpha * p2;
1004 points[k] += (1. - alpha) * p1;
1005 weights[k] = length * quadrature.weight(k);
1006 }
1007 return Quadrature<dim>(points, weights);
1008}
1009
1010
1011
1012template <int dim>
1015 const unsigned int face_no,
1016 const bool face_orientation,
1017 const bool face_flip,
1018 const bool face_rotation,
1019 const unsigned int n_quadrature_points)
1020{
1021 return face(reference_cell,
1022 face_no,
1023 internal::combined_face_orientation(face_orientation,
1024 face_rotation,
1025 face_flip),
1026 n_quadrature_points);
1027}
1028
1029
1030
1031template <int dim>
1034 const ReferenceCell &reference_cell,
1035 const unsigned int face_no,
1036 const types::geometric_orientation combined_orientation,
1037 const unsigned int n_quadrature_points)
1038{
1039 AssertIndexRange(face_no, reference_cell.n_faces());
1040 AssertIndexRange(combined_orientation,
1041 reference_cell.n_face_orientations(face_no));
1042 AssertDimension(reference_cell.get_dimension(), dim);
1043
1044
1045 return {(reference_cell.n_face_orientations(face_no) * face_no +
1046 combined_orientation) *
1047 n_quadrature_points};
1048}
1049
1050
1051
1052template <int dim>
1055 const ReferenceCell &reference_cell,
1056 const unsigned int face_no,
1057 const bool face_orientation,
1058 const bool face_flip,
1059 const bool face_rotation,
1060 const hp::QCollection<dim - 1> &quadrature)
1061{
1062 return face(reference_cell,
1063 face_no,
1064 internal::combined_face_orientation(face_orientation,
1065 face_rotation,
1066 face_flip),
1067 quadrature);
1068}
1069
1070
1071
1072template <int dim>
1075 const ReferenceCell &reference_cell,
1076 const unsigned int face_no,
1077 const types::geometric_orientation combined_orientation,
1078 const hp::QCollection<dim - 1> &quadrature)
1079{
1080 AssertIndexRange(face_no, reference_cell.n_faces());
1081 AssertIndexRange(combined_orientation,
1082 reference_cell.n_face_orientations(face_no));
1083 AssertDimension(reference_cell.get_dimension(), dim);
1084
1085 unsigned int offset = 0;
1086 if (quadrature.size() == 1)
1087 offset =
1088 reference_cell.n_face_orientations(0) * quadrature[0].size() * face_no;
1089 else
1090 for (unsigned int i = 0; i < face_no; ++i)
1091 offset += reference_cell.n_face_orientations(i) * quadrature[i].size();
1092
1093 return {offset + combined_orientation *
1094 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1095}
1096
1097
1098
1099template <int dim>
1102 const ReferenceCell &reference_cell,
1103 const unsigned int face_no,
1104 const unsigned int subface_no,
1105 const bool face_orientation,
1106 const bool face_flip,
1107 const bool face_rotation,
1108 const unsigned int n_quadrature_points,
1109 const internal::SubfaceCase<dim> ref_case)
1110{
1112 reference_cell,
1113 face_no,
1114 subface_no,
1115 internal::combined_face_orientation(face_orientation,
1116 face_rotation,
1117 face_flip),
1118 n_quadrature_points,
1119 ref_case);
1120}
1121
1122
1123
1124template <>
1127 const ReferenceCell &reference_cell,
1128 const unsigned int face_no,
1129 const unsigned int subface_no,
1130 const types::geometric_orientation /*combined_orientation*/,
1131 const unsigned int n_quadrature_points,
1133{
1134 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1135 (void)reference_cell;
1136
1140
1141 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1142 n_quadrature_points);
1143}
1144
1145
1146
1147template <>
1150 const ReferenceCell &reference_cell,
1151 const unsigned int face_no,
1152 const unsigned int subface_no,
1153 const types::geometric_orientation combined_orientation,
1154 const unsigned int n_quadrature_points,
1156{
1157 Assert(reference_cell == ReferenceCells::Quadrilateral ||
1158 reference_cell == ReferenceCells::Triangle,
1160
1161 const unsigned int n_faces = reference_cell.n_faces();
1162 const unsigned int n_subfaces =
1163 reference_cell.face_reference_cell(face_no).n_isotropic_children();
1164 const unsigned int n_orientations =
1165 reference_cell.n_face_orientations(face_no);
1166
1167 AssertIndexRange(face_no, n_faces);
1168 AssertIndexRange(subface_no, n_subfaces);
1169 AssertIndexRange(combined_orientation, n_orientations);
1170
1171 return ((face_no * n_orientations * n_subfaces +
1172 combined_orientation * n_subfaces + subface_no) *
1173 n_quadrature_points);
1174}
1175
1176
1177
1178template <>
1181 const ReferenceCell &reference_cell,
1182 const unsigned int face_no,
1183 const unsigned int subface_no,
1184 const types::geometric_orientation combined_orientation,
1185 const unsigned int n_quadrature_points,
1186 const internal::SubfaceCase<3> ref_case)
1187{
1188 const unsigned int dim = 3;
1189
1190 const auto [face_orientation, face_rotation, face_flip] =
1191 internal::split_face_orientation(combined_orientation);
1192
1194 (void)reference_cell;
1195
1199
1200 // in 3d, we have to account for faces that
1201 // have non-standard orientation. thus, we
1202 // have to store _eight_ data sets per face
1203 // or subface already for the isotropic
1204 // case. Additionally, we have three
1205 // different refinement cases, resulting in
1206 // <tt>4 + 2 + 2 = 8</tt> different subfaces
1207 // for each face.
1208 const unsigned int total_subfaces_per_face = 8;
1209
1210 // set up a table with the offsets for a
1211 // given refinement case respecting the
1212 // corresponding number of subfaces. the
1213 // index corresponds to (RefineCase::Type - 1)
1214
1215 // note, that normally we should use the
1216 // obvious offsets 0,2,6. However, prior to
1217 // the implementation of anisotropic
1218 // refinement, in many places of the library
1219 // the convention was used, that the first
1220 // dataset with offset 0 corresponds to a
1221 // standard (isotropic) face
1222 // refinement. therefore we use the offsets
1223 // 6,4,0 here to stick to that (implicit)
1224 // convention
1225 static const unsigned int ref_case_offset[3] = {
1226 6, // cut_x
1227 4, // cut_y
1228 0 // cut_xy
1229 };
1230
1231 const std::pair<unsigned int, RefinementCase<2>>
1232 final_subface_no_and_ref_case =
1233 internal::QProjector::select_subface_no_and_refinement_case(
1234 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1235
1236 return ((face_no * total_subfaces_per_face +
1237 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1238 final_subface_no_and_ref_case.first) +
1239 GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face *
1240 combined_orientation) *
1241 n_quadrature_points;
1242}
1243
1244
1245
1246template <int dim>
1249 const SubQuadrature &quadrature,
1250 const unsigned int face_no)
1251{
1252 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1254 (void)reference_cell;
1255
1256 std::vector<Point<dim>> points(quadrature.size());
1257 project_to_face(reference_cell, quadrature, face_no, points);
1258 return Quadrature<dim>(points, quadrature.get_weights());
1259}
1260
1261
1262
1263template <int dim>
1266 const SubQuadrature &quadrature,
1267 const unsigned int face_no,
1268 const unsigned int subface_no,
1269 const RefinementCase<dim - 1> &ref_case)
1270{
1271 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1273 (void)reference_cell;
1274
1275 std::vector<Point<dim>> points(quadrature.size());
1277 reference_cell, quadrature, face_no, subface_no, points, ref_case);
1278 return Quadrature<dim>(points, quadrature.get_weights());
1279}
1280
1281
1282// explicit instantiations; note: we need them all for all dimensions
1283template class QProjector<1>;
1284template class QProjector<2>;
1285template class QProjector<3>;
1286
auto make_const_array_view(const Container &container) -> decltype(make_array_view(container))
Definition point.h:113
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Class storing the offset index into a Quadrature rule created by project_to_all_faces() or project_to...
Definition qprojector.h:334
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)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:69
static Quadrature< dim > project_to_all_subfaces(const ReferenceCell &reference_cell, const SubQuadrature &quadrature)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &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_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
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
unsigned int n_vertices() const
double d_linear_shape_function(const Point< dim > &xi, const unsigned int i) 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
double volume() const
unsigned int size() const
Definition collection.h:316
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
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)
constexpr ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Tetrahedron
constexpr 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:238
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:365
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
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)