Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3087-ga35b476a78 2025-04-19 20:30:01+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
190 mutate_quadrature(const Quadrature<2> &quadrature,
191 const types::geometric_orientation combined_orientation)
192 {
193 return Quadrature<2>(mutate_points_with_offset(quadrature.get_points(),
194 combined_orientation),
195 quadrature.get_weights());
196 }
197
198 std::pair<unsigned int, RefinementCase<2>>
199 select_subface_no_and_refinement_case(
200 const unsigned int subface_no,
201 const bool face_orientation,
202 const bool face_flip,
203 const bool face_rotation,
204 const internal::SubfaceCase<3> ref_case)
205 {
206 constexpr int dim = 3;
207 // for each subface of a given FaceRefineCase
208 // there is a corresponding equivalent
209 // subface number of one of the "standard"
210 // RefineCases (cut_x, cut_y, cut_xy). Map
211 // the given values to those equivalent
212 // ones.
213
214 // first, define an invalid number
215 static const unsigned int e = numbers::invalid_unsigned_int;
216
217 static const RefinementCase<dim - 1>
218 equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
220 // case_none. there should be only
221 // invalid values here. However, as
222 // this function is also called (in
223 // tests) for cells which have no
224 // refined faces, use isotropic
225 // refinement instead
226 {RefinementCase<dim - 1>::cut_xy,
227 RefinementCase<dim - 1>::cut_xy,
228 RefinementCase<dim - 1>::cut_xy,
229 RefinementCase<dim - 1>::cut_xy},
230 // case_x
231 {RefinementCase<dim - 1>::cut_x,
232 RefinementCase<dim - 1>::cut_x,
233 RefinementCase<dim - 1>::no_refinement,
234 RefinementCase<dim - 1>::no_refinement},
235 // case_x1y
236 {RefinementCase<dim - 1>::cut_xy,
237 RefinementCase<dim - 1>::cut_xy,
238 RefinementCase<dim - 1>::cut_x,
239 RefinementCase<dim - 1>::no_refinement},
240 // case_x2y
241 {RefinementCase<dim - 1>::cut_x,
242 RefinementCase<dim - 1>::cut_xy,
243 RefinementCase<dim - 1>::cut_xy,
244 RefinementCase<dim - 1>::no_refinement},
245 // case_x1y2y
246 {RefinementCase<dim - 1>::cut_xy,
247 RefinementCase<dim - 1>::cut_xy,
248 RefinementCase<dim - 1>::cut_xy,
249 RefinementCase<dim - 1>::cut_xy},
250 // case_y
251 {RefinementCase<dim - 1>::cut_y,
252 RefinementCase<dim - 1>::cut_y,
253 RefinementCase<dim - 1>::no_refinement,
254 RefinementCase<dim - 1>::no_refinement},
255 // case_y1x
256 {RefinementCase<dim - 1>::cut_xy,
257 RefinementCase<dim - 1>::cut_xy,
258 RefinementCase<dim - 1>::cut_y,
259 RefinementCase<dim - 1>::no_refinement},
260 // case_y2x
261 {RefinementCase<dim - 1>::cut_y,
262 RefinementCase<dim - 1>::cut_xy,
263 RefinementCase<dim - 1>::cut_xy,
264 RefinementCase<dim - 1>::no_refinement},
265 // case_y1x2x
266 {RefinementCase<dim - 1>::cut_xy,
267 RefinementCase<dim - 1>::cut_xy,
268 RefinementCase<dim - 1>::cut_xy,
269 RefinementCase<dim - 1>::cut_xy},
270 // case_xy (case_isotropic)
271 {RefinementCase<dim - 1>::cut_xy,
272 RefinementCase<dim - 1>::cut_xy,
273 RefinementCase<dim - 1>::cut_xy,
274 RefinementCase<dim - 1>::cut_xy}};
275
276 static const unsigned int
277 equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic +
279 {// case_none, see above
280 {0, 1, 2, 3},
281 // case_x
282 {0, 1, e, e},
283 // case_x1y
284 {0, 2, 1, e},
285 // case_x2y
286 {0, 1, 3, e},
287 // case_x1y2y
288 {0, 2, 1, 3},
289 // case_y
290 {0, 1, e, e},
291 // case_y1x
292 {0, 1, 1, e},
293 // case_y2x
294 {0, 2, 3, e},
295 // case_y1x2x
296 {0, 1, 2, 3},
297 // case_xy (case_isotropic)
298 {0, 1, 2, 3}};
299
300 // If face-orientation or face_rotation are
301 // non-standard, cut_x and cut_y have to be
302 // exchanged.
303 static const RefinementCase<dim - 1> ref_case_permutation[4] = {
304 RefinementCase<dim - 1>::no_refinement,
305 RefinementCase<dim - 1>::cut_y,
306 RefinementCase<dim - 1>::cut_x,
307 RefinementCase<dim - 1>::cut_xy};
308
309 // set a corresponding (equivalent)
310 // RefineCase and subface number
311 const RefinementCase<dim - 1> equ_ref_case =
312 equivalent_refine_case[ref_case][subface_no];
313 const unsigned int equ_subface_no =
314 equivalent_subface_number[ref_case][subface_no];
315 // make sure, that we got a valid subface and RefineCase
318 Assert(equ_subface_no != e, ExcInternalError());
319 // now, finally respect non-standard faces
320 const RefinementCase<dim - 1> final_ref_case =
321 (face_orientation == face_rotation ?
322 ref_case_permutation[equ_ref_case] :
323 equ_ref_case);
324
325 const unsigned int final_subface_no =
327 final_ref_case),
328 4,
329 equ_subface_no,
330 face_orientation,
331 face_flip,
332 face_rotation,
333 equ_ref_case);
334
335 return std::make_pair(final_subface_no, final_ref_case);
336 }
337 } // namespace
338 } // namespace QProjector
339} // namespace internal
340
341
342
343template <>
344void
346 const Quadrature<0> &quadrature,
347 const unsigned int face_no,
348 std::vector<Point<1>> &q_points)
349{
350 AssertDimension(quadrature.size(), q_points.size());
351 const auto face_quadrature =
352 QProjector<1>::project_to_face(reference_cell,
353 quadrature,
354 face_no,
356 q_points = face_quadrature.get_points();
357}
358
359
360
361template <>
362void
364 const Quadrature<1> &quadrature,
365 const unsigned int face_no,
366 std::vector<Point<2>> &q_points)
367{
368 AssertDimension(quadrature.size(), q_points.size());
369 const auto face_quadrature =
370 QProjector<2>::project_to_face(reference_cell,
371 quadrature,
372 face_no,
374 q_points = face_quadrature.get_points();
375}
376
377
378
379template <>
380void
382 const Quadrature<2> &quadrature,
383 const unsigned int face_no,
384 std::vector<Point<3>> &q_points)
385{
387 (void)reference_cell;
388
390 Assert(q_points.size() == quadrature.size(),
391 ExcDimensionMismatch(q_points.size(), quadrature.size()));
392 q_points.clear();
393 internal::QProjector::project_to_hex_face_and_append(quadrature.get_points(),
394 face_no,
395 q_points);
396}
397
398
399
400template <int dim>
403 const Quadrature<dim - 1> &quadrature,
404 const unsigned int face_no,
405 const bool face_orientation,
406 const bool face_flip,
407 const bool face_rotation)
408{
410 reference_cell,
411 quadrature,
412 face_no,
414 face_rotation,
415 face_flip));
416}
417
418
419
420template <int dim>
423 const Quadrature<dim - 1> &quadrature,
424 const unsigned int face_no,
425 const types::geometric_orientation orientation)
426{
427 AssertIndexRange(face_no, reference_cell.n_faces());
428 AssertIndexRange(orientation, reference_cell.n_face_orientations(face_no));
429 AssertDimension(reference_cell.get_dimension(), dim);
430
431 std::vector<Point<dim>> q_points;
432 std::vector<double> q_weights = quadrature.get_weights();
433 q_points.reserve(quadrature.size());
434 if constexpr (dim == 1)
435 {
436 AssertDimension(quadrature.size(), 1);
437 q_points.emplace_back(static_cast<double>(face_no));
438 }
439 else if constexpr (dim == 2)
440 {
441 if (reference_cell == ReferenceCells::Triangle)
442 // use linear polynomial to map the reference quadrature points
443 // correctly on faces, i.e., BarycentricPolynomials<1>(1)
444 for (unsigned int p = 0; p < quadrature.size(); ++p)
445 {
446 if (face_no == 0)
447 q_points.emplace_back(quadrature.point(p)[0], 0);
448 else if (face_no == 1)
449 q_points.emplace_back(1.0 - quadrature.point(p)[0],
450 quadrature.point(p)[0]);
451 else if (face_no == 2)
452 q_points.emplace_back(0, 1.0 - quadrature.point(p)[0]);
453 else
455 }
456 else if (reference_cell == ReferenceCells::Quadrilateral)
457 for (unsigned int p = 0; p < quadrature.size(); ++p)
458 {
459 if (face_no == 0)
460 q_points.emplace_back(0, quadrature.point(p)[0]);
461 else if (face_no == 1)
462 q_points.emplace_back(1, quadrature.point(p)[0]);
463 else if (face_no == 2)
464 q_points.emplace_back(quadrature.point(p)[0], 0);
465 else if (face_no == 3)
466 q_points.emplace_back(quadrature.point(p)[0], 1);
467 else
469 }
470 else
472
473 if (orientation == numbers::reverse_line_orientation)
474 {
475 std::reverse(q_points.begin(), q_points.end());
476 std::reverse(q_weights.begin(), q_weights.end());
477 }
478 }
479 else if constexpr (dim == 3)
480 {
482 const Quadrature<2> mutation =
483 internal::QProjector::mutate_quadrature(quadrature, orientation);
484 return QProjector<3>::project_to_face(reference_cell, mutation, face_no);
485 }
486 else
487 {
489 }
490
491 return Quadrature<dim>(std::move(q_points), std::move(q_weights));
492}
493
494
495
496template <>
497void
499 const Quadrature<0> &,
500 const unsigned int face_no,
501 const unsigned int,
502 std::vector<Point<1>> &q_points,
503 const RefinementCase<0> &)
504{
505 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
506 (void)reference_cell;
507
508 const unsigned int dim = 1;
510 AssertDimension(q_points.size(), 1);
511
512 q_points[0] = Point<dim>(static_cast<double>(face_no));
513}
514
515
516
517template <>
518void
520 const Quadrature<1> &quadrature,
521 const unsigned int face_no,
522 const unsigned int subface_no,
523 std::vector<Point<2>> &q_points,
524 const RefinementCase<1> &ref_case)
525{
526 AssertDimension(quadrature.size(), q_points.size());
527 const auto face_quadrature =
528 project_to_subface(reference_cell,
529 quadrature,
530 face_no,
531 subface_no,
533 ref_case);
534 q_points = face_quadrature.get_points();
535}
536
537
538
539template <>
540void
542 const Quadrature<2> &quadrature,
543 const unsigned int face_no,
544 const unsigned int subface_no,
545 std::vector<Point<3>> &q_points,
546 const RefinementCase<2> &ref_case)
547{
549 (void)reference_cell;
550 AssertDimension(quadrature.size(), q_points.size());
551 const auto face_quadrature =
552 project_to_subface(reference_cell,
553 quadrature,
554 face_no,
555 subface_no,
557 ref_case);
558 q_points = face_quadrature.get_points();
559}
560
561
562
563template <int dim>
566 const ReferenceCell &reference_cell,
567 const Quadrature<dim - 1> &quadrature,
568 const unsigned int face_no,
569 const unsigned int subface_no,
570 const bool,
571 const bool,
572 const bool,
574{
576 reference_cell,
577 quadrature,
578 face_no,
579 subface_no,
581}
582
583
584
585template <int dim>
588 const ReferenceCell &reference_cell,
589 const SubQuadrature &quadrature,
590 const unsigned int face_no,
591 const unsigned int subface_no,
592 const types::geometric_orientation orientation,
593 const RefinementCase<dim - 1> &ref_case)
594{
595 AssertIndexRange(face_no, reference_cell.n_faces());
596 AssertIndexRange(orientation, reference_cell.n_face_orientations(face_no));
597 AssertDimension(reference_cell.get_dimension(), dim);
598 AssertIndexRange(subface_no,
599 reference_cell.face_reference_cell(face_no)
600 .template n_children<dim - 1>(ref_case));
601
602 std::vector<Point<dim>> q_points;
603 std::vector<double> q_weights = quadrature.get_weights();
604 q_points.reserve(quadrature.size());
605
606 if constexpr (dim == 1)
607 {
608 AssertDimension(quadrature.size(), 1);
609 q_points.emplace_back(static_cast<double>(face_no));
610 }
611 else if constexpr (dim == 2)
612 {
613 if (reference_cell == ReferenceCells::Triangle)
614 // use linear polynomial to map the reference quadrature points
615 // correctly on faces, i.e., BarycentricPolynomials<1>(1)
616 for (unsigned int p = 0; p < quadrature.size(); ++p)
617 {
618 if (face_no == 0)
619 {
620 if (subface_no == 0)
621 q_points.emplace_back(quadrature.point(p)[0] / 2, 0);
622 else
623 q_points.emplace_back(0.5 + quadrature.point(p)[0] / 2, 0);
624 }
625 else if (face_no == 1)
626 {
627 if (subface_no == 0)
628 q_points.emplace_back(1 - quadrature.point(p)[0] / 2,
629 quadrature.point(p)[0] / 2);
630 else
631 q_points.emplace_back(0.5 - quadrature.point(p)[0] / 2,
632 0.5 + quadrature.point(p)[0] / 2);
633 }
634 else if (face_no == 2)
635 {
636 if (subface_no == 0)
637 q_points.emplace_back(0, 1 - quadrature.point(p)[0] / 2);
638 else
639 q_points.emplace_back(0, 0.5 - quadrature.point(p)[0] / 2);
640 }
641 else
643 }
644 else if (reference_cell == ReferenceCells::Quadrilateral)
645 for (unsigned int p = 0; p < quadrature.size(); ++p)
646 {
647 if (face_no == 0)
648 {
649 if (subface_no == 0)
650 q_points.emplace_back(0, quadrature.point(p)[0] / 2);
651 else
652 q_points.emplace_back(0, quadrature.point(p)[0] / 2 + 0.5);
653 }
654 else if (face_no == 1)
655 {
656 if (subface_no == 0)
657 q_points.emplace_back(1, quadrature.point(p)[0] / 2);
658 else
659 q_points.emplace_back(1, quadrature.point(p)[0] / 2 + 0.5);
660 }
661 else if (face_no == 2)
662 {
663 if (subface_no == 0)
664 q_points.emplace_back(quadrature.point(p)[0] / 2, 0);
665 else
666 q_points.emplace_back(quadrature.point(p)[0] / 2 + 0.5, 0);
667 }
668 else if (face_no == 3)
669 {
670 if (subface_no == 0)
671 q_points.emplace_back(quadrature.point(p)[0] / 2, 1);
672 else
673 q_points.emplace_back(quadrature.point(p)[0] / 2 + 0.5, 1);
674 }
675 else
677 }
678 else
680
681 if (orientation == numbers::reverse_line_orientation)
682 {
683 std::reverse(q_points.begin(), q_points.end());
684 std::reverse(q_weights.begin(), q_weights.end());
685 }
686 }
687 else if constexpr (dim == 3)
688 {
690 internal::QProjector::project_to_hex_face_and_append(
691 quadrature.get_points(), face_no, q_points, ref_case, subface_no);
692 }
693 else
694 {
696 }
697
698 return Quadrature<dim>(std::move(q_points), std::move(q_weights));
699}
700
701
702
703template <int dim>
706 const ReferenceCell &reference_cell,
707 const hp::QCollection<dim - 1> &quadrature)
708{
709 const auto process = [&](const std::vector<std::vector<Point<dim>>> &faces) {
710 // new (projected) quadrature points and weights
711 std::vector<Point<dim>> points;
712 std::vector<double> weights;
713
714 // loop over all faces (triangles) ...
715 for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
716 {
717 const ReferenceCell face_reference_cell =
718 reference_cell.face_reference_cell(face_no);
719
720 // ... and over all possible orientations
721 for (types::geometric_orientation orientation = 0;
722 orientation < reference_cell.n_face_orientations(face_no);
723 ++orientation)
724 {
725 const auto &face = faces[face_no];
726
727 // The goal of this function is to compute identical sets of
728 // quadrature points on the common face of two abutting cells. Our
729 // orientation convention is that, given such a pair of abutting
730 // cells:
731 //
732 // 1. The shared face, from the perspective of the first cell, is
733 // in the default orientation.
734 // 2. The shared face, from the perspective of the second cell, has
735 // its orientation computed relative to the first cell: i.e.,
736 // 'orientation' is the vertex permutation applied to the first
737 // cell's face to get the second cell's face.
738 //
739 // The first case is trivial since points do not need to be
740 // oriented. However, in the second case, we need to use the
741 // *reverse* of the stored orientation (i.e., the permutation
742 // applied to the second cell's face which yields the first cell's
743 // face) so that we get identical quadrature points.
744 //
745 // For more information see connectivity.h.
746 const boost::container::small_vector<Point<dim>, 8> support_points =
747 face_reference_cell.permute_by_combined_orientation<Point<dim>>(
748 face,
749 face_reference_cell.get_inverse_combined_orientation(
750 orientation));
751
752 // the quadrature rule to be projected ...
753 const auto &sub_quadrature_points =
754 quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
755 const auto &sub_quadrature_weights =
756 quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
757
758 // loop over all quadrature points
759 for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
760 {
761 Point<dim> mapped_point;
762
763 // map reference quadrature point
764 for (const unsigned int i :
765 face_reference_cell.vertex_indices())
766 mapped_point += support_points[i] *
767 face_reference_cell.d_linear_shape_function(
768 sub_quadrature_points[j], i);
769
770 points.push_back(mapped_point);
771
772 // rescale quadrature weights so that the sum of the weights on
773 // each face equals the measure of that face.
774 const double scaling = reference_cell.face_measure(face_no) /
775 face_reference_cell.volume();
776 weights.push_back(sub_quadrature_weights[j] * scaling);
777 }
778 }
779 }
780
781 // construct new quadrature rule
782 return Quadrature<dim>(std::move(points), std::move(weights));
783 };
784
785 std::vector<std::vector<Point<dim>>> face_vertex_locations(
786 reference_cell.n_faces());
787 for (const unsigned int f : reference_cell.face_indices())
788 {
789 face_vertex_locations[f].resize(
790 reference_cell.face_reference_cell(f).n_vertices());
791 for (const unsigned int v :
792 reference_cell.face_reference_cell(f).vertex_indices())
793 face_vertex_locations[f][v] =
794 reference_cell.face_vertex_location<dim>(f, v);
795 }
796
797 return process(face_vertex_locations);
798}
799
800
801
802template <>
805 const Quadrature<0> &quadrature)
806{
807 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
808 (void)reference_cell;
809
810 const unsigned int dim = 1;
811
812 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
813 subfaces_per_face =
815
816 // first fix quadrature points
817 std::vector<Point<dim>> q_points;
818 q_points.reserve(n_points * n_faces * subfaces_per_face);
819 std::vector<Point<dim>> help(n_points);
820
821 // project to each face and copy
822 // results
823 for (unsigned int face = 0; face < n_faces; ++face)
824 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
825 {
826 project_to_subface(reference_cell, quadrature, face, subface, help);
827 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
828 }
829
830 // next copy over weights
831 std::vector<double> weights;
832 weights.reserve(n_points * n_faces * subfaces_per_face);
833 for (unsigned int face = 0; face < n_faces; ++face)
834 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
835 std::copy(quadrature.get_weights().begin(),
836 quadrature.get_weights().end(),
837 std::back_inserter(weights));
838
839 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
841 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
843
844 return Quadrature<dim>(std::move(q_points), std::move(weights));
845}
846
847
848
849template <>
852 const SubQuadrature &quadrature)
853{
854 if (reference_cell == ReferenceCells::Triangle ||
855 reference_cell == ReferenceCells::Tetrahedron)
856 return Quadrature<2>(); // nothing to do
857
859
860 const unsigned int dim = 2;
861
862 const unsigned int n_points = quadrature.size(),
864 subfaces_per_face =
866
867 // first fix quadrature points
868 std::vector<Point<dim>> q_points;
869 q_points.reserve(n_points * n_faces * subfaces_per_face);
870 std::vector<Point<dim>> help(n_points);
871
872 // project to each face and copy
873 // results
874 for (unsigned int face = 0; face < n_faces; ++face)
875 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
876 {
877 project_to_subface(reference_cell, quadrature, face, subface, help);
878 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
879 }
880
881 // next copy over weights
882 std::vector<double> weights;
883 weights.reserve(n_points * n_faces * subfaces_per_face);
884 for (unsigned int face = 0; face < n_faces; ++face)
885 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
886 std::copy(quadrature.get_weights().begin(),
887 quadrature.get_weights().end(),
888 std::back_inserter(weights));
889
890 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
892 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
894
895 return Quadrature<dim>(std::move(q_points), std::move(weights));
896}
897
898
899
900template <>
903 const SubQuadrature &quadrature)
904{
905 if (reference_cell == ReferenceCells::Triangle ||
906 reference_cell == ReferenceCells::Tetrahedron)
907 return Quadrature<3>(); // nothing to do
908
910
911 const unsigned int dim = 3;
912
913 const unsigned int n_points = quadrature.size(),
915 total_subfaces_per_face = 2 + 2 + 4;
916
917 // first fix quadrature points
918 std::vector<Point<dim>> q_points;
919 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
920
921 std::vector<double> weights;
922 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
923
924 // do the following for all possible mutations of a face (mutation==0
925 // corresponds to a face with standard orientation, no flip and no rotation)
926 for (unsigned char offset = 0; offset < 8; ++offset)
927 {
928 const auto mutation =
929 internal::QProjector::mutate_points_with_offset(quadrature.get_points(),
930 offset);
931
932 // project to each face and copy results
933 for (unsigned int face = 0; face < n_faces; ++face)
934 for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
935 ref_case >= RefinementCase<dim - 1>::cut_x;
936 --ref_case)
937 for (unsigned int subface = 0;
939 RefinementCase<dim - 1>(ref_case));
940 ++subface)
941 {
942 internal::QProjector::project_to_hex_face_and_append(
943 mutation,
944 face,
945 q_points,
946 RefinementCase<dim - 1>(ref_case),
947 subface);
948
949 // next copy over weights
950 std::copy(quadrature.get_weights().begin(),
951 quadrature.get_weights().end(),
952 std::back_inserter(weights));
953 }
954 }
955
956 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
958 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
960
961 return Quadrature<dim>(std::move(q_points), std::move(weights));
962}
963
964
965
966template <int dim>
969 const Quadrature<dim> &quadrature,
970 const unsigned int child_no)
971{
972 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
974 (void)reference_cell;
975
977
978 const unsigned int n_q_points = quadrature.size();
979
980 std::vector<Point<dim>> q_points(n_q_points);
981 for (unsigned int i = 0; i < n_q_points; ++i)
982 q_points[i] =
984 child_no);
985
986 // for the weights, things are
987 // equally simple: copy them and
988 // scale them
989 std::vector<double> weights = quadrature.get_weights();
990 for (unsigned int i = 0; i < n_q_points; ++i)
991 weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
992
993 return Quadrature<dim>(q_points, weights);
994}
995
996
997
998template <int dim>
1001 const Quadrature<dim> &quadrature)
1002{
1003 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1005 (void)reference_cell;
1006
1007 const unsigned int n_points = quadrature.size(),
1009
1010 std::vector<Point<dim>> q_points(n_points * n_children);
1011 std::vector<double> weights(n_points * n_children);
1012
1013 // project to each child and copy
1014 // results
1015 for (unsigned int child = 0; child < n_children; ++child)
1016 {
1017 Quadrature<dim> help =
1018 project_to_child(reference_cell, quadrature, child);
1019 for (unsigned int i = 0; i < n_points; ++i)
1020 {
1021 q_points[child * n_points + i] = help.point(i);
1022 weights[child * n_points + i] = help.weight(i);
1023 }
1024 }
1025 return Quadrature<dim>(q_points, weights);
1026}
1027
1028
1029
1030template <int dim>
1033 const Quadrature<1> &quadrature,
1034 const Point<dim> &p1,
1035 const Point<dim> &p2)
1036{
1037 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1039 (void)reference_cell;
1040
1041 const unsigned int n = quadrature.size();
1042 std::vector<Point<dim>> points(n);
1043 std::vector<double> weights(n);
1044 const double length = p1.distance(p2);
1045
1046 for (unsigned int k = 0; k < n; ++k)
1047 {
1048 const double alpha = quadrature.point(k)[0];
1049 points[k] = alpha * p2;
1050 points[k] += (1. - alpha) * p1;
1051 weights[k] = length * quadrature.weight(k);
1052 }
1053 return Quadrature<dim>(points, weights);
1054}
1055
1056
1057
1058template <int dim>
1061 const unsigned int face_no,
1062 const bool face_orientation,
1063 const bool face_flip,
1064 const bool face_rotation,
1065 const unsigned int n_quadrature_points)
1066{
1067 return face(reference_cell,
1068 face_no,
1069 internal::combined_face_orientation(face_orientation,
1070 face_rotation,
1071 face_flip),
1072 n_quadrature_points);
1073}
1074
1075
1076
1077template <int dim>
1080 const ReferenceCell &reference_cell,
1081 const unsigned int face_no,
1082 const types::geometric_orientation combined_orientation,
1083 const unsigned int n_quadrature_points)
1084{
1085 AssertIndexRange(face_no, reference_cell.n_faces());
1086 AssertIndexRange(combined_orientation,
1087 reference_cell.n_face_orientations(face_no));
1088 AssertDimension(reference_cell.get_dimension(), dim);
1089
1090
1091 return {(reference_cell.n_face_orientations(face_no) * face_no +
1092 combined_orientation) *
1093 n_quadrature_points};
1094}
1095
1096
1097
1098template <int dim>
1101 const ReferenceCell &reference_cell,
1102 const unsigned int face_no,
1103 const bool face_orientation,
1104 const bool face_flip,
1105 const bool face_rotation,
1106 const hp::QCollection<dim - 1> &quadrature)
1107{
1108 return face(reference_cell,
1109 face_no,
1110 internal::combined_face_orientation(face_orientation,
1111 face_rotation,
1112 face_flip),
1113 quadrature);
1114}
1115
1116
1117
1118template <int dim>
1121 const ReferenceCell &reference_cell,
1122 const unsigned int face_no,
1123 const types::geometric_orientation combined_orientation,
1124 const hp::QCollection<dim - 1> &quadrature)
1125{
1126 AssertIndexRange(face_no, reference_cell.n_faces());
1127 AssertIndexRange(combined_orientation,
1128 reference_cell.n_face_orientations(face_no));
1129 AssertDimension(reference_cell.get_dimension(), dim);
1130
1131 unsigned int offset = 0;
1132 if (quadrature.size() == 1)
1133 offset =
1134 reference_cell.n_face_orientations(0) * quadrature[0].size() * face_no;
1135 else
1136 for (unsigned int i = 0; i < face_no; ++i)
1137 offset += reference_cell.n_face_orientations(i) * quadrature[i].size();
1138
1139 return {offset + combined_orientation *
1140 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1141}
1142
1143
1144
1145template <int dim>
1148 const ReferenceCell &reference_cell,
1149 const unsigned int face_no,
1150 const unsigned int subface_no,
1151 const bool face_orientation,
1152 const bool face_flip,
1153 const bool face_rotation,
1154 const unsigned int n_quadrature_points,
1155 const internal::SubfaceCase<dim> ref_case)
1156{
1158 reference_cell,
1159 face_no,
1160 subface_no,
1161 internal::combined_face_orientation(face_orientation,
1162 face_rotation,
1163 face_flip),
1164 n_quadrature_points,
1165 ref_case);
1166}
1167
1168
1169
1170template <>
1173 const ReferenceCell &reference_cell,
1174 const unsigned int face_no,
1175 const unsigned int subface_no,
1176 const types::geometric_orientation /*combined_orientation*/,
1177 const unsigned int n_quadrature_points,
1179{
1180 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1181 (void)reference_cell;
1182
1186
1187 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1188 n_quadrature_points);
1189}
1190
1191
1192
1193template <>
1196 const ReferenceCell &reference_cell,
1197 const unsigned int face_no,
1198 const unsigned int subface_no,
1199 const types::geometric_orientation /*combined_orientation*/,
1200 const unsigned int n_quadrature_points,
1202{
1204 (void)reference_cell;
1205
1209
1210 return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1211 n_quadrature_points);
1212}
1213
1214
1215
1216template <>
1219 const ReferenceCell &reference_cell,
1220 const unsigned int face_no,
1221 const unsigned int subface_no,
1222 const types::geometric_orientation combined_orientation,
1223 const unsigned int n_quadrature_points,
1224 const internal::SubfaceCase<3> ref_case)
1225{
1226 const unsigned int dim = 3;
1227
1228 const auto [face_orientation, face_rotation, face_flip] =
1229 internal::split_face_orientation(combined_orientation);
1230
1232 (void)reference_cell;
1233
1237
1238 // in 3d, we have to account for faces that
1239 // have non-standard orientation. thus, we
1240 // have to store _eight_ data sets per face
1241 // or subface already for the isotropic
1242 // case. Additionally, we have three
1243 // different refinement cases, resulting in
1244 // <tt>4 + 2 + 2 = 8</tt> different subfaces
1245 // for each face.
1246 const unsigned int total_subfaces_per_face = 8;
1247
1248 // set up a table with the offsets for a
1249 // given refinement case respecting the
1250 // corresponding number of subfaces. the
1251 // index corresponds to (RefineCase::Type - 1)
1252
1253 // note, that normally we should use the
1254 // obvious offsets 0,2,6. However, prior to
1255 // the implementation of anisotropic
1256 // refinement, in many places of the library
1257 // the convention was used, that the first
1258 // dataset with offset 0 corresponds to a
1259 // standard (isotropic) face
1260 // refinement. therefore we use the offsets
1261 // 6,4,0 here to stick to that (implicit)
1262 // convention
1263 static const unsigned int ref_case_offset[3] = {
1264 6, // cut_x
1265 4, // cut_y
1266 0 // cut_xy
1267 };
1268
1269 const std::pair<unsigned int, RefinementCase<2>>
1270 final_subface_no_and_ref_case =
1271 internal::QProjector::select_subface_no_and_refinement_case(
1272 subface_no, face_orientation, face_flip, face_rotation, ref_case);
1273
1274 return ((face_no * total_subfaces_per_face +
1275 ref_case_offset[final_subface_no_and_ref_case.second - 1] +
1276 final_subface_no_and_ref_case.first) +
1277 GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face *
1278 combined_orientation) *
1279 n_quadrature_points;
1280}
1281
1282
1283
1284template <int dim>
1287 const SubQuadrature &quadrature,
1288 const unsigned int face_no)
1289{
1290 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1292 (void)reference_cell;
1293
1294 std::vector<Point<dim>> points(quadrature.size());
1295 project_to_face(reference_cell, quadrature, face_no, points);
1296 return Quadrature<dim>(points, quadrature.get_weights());
1297}
1298
1299
1300
1301template <int dim>
1304 const SubQuadrature &quadrature,
1305 const unsigned int face_no,
1306 const unsigned int subface_no,
1307 const RefinementCase<dim - 1> &ref_case)
1308{
1309 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1311 (void)reference_cell;
1312
1313 std::vector<Point<dim>> points(quadrature.size());
1315 reference_cell, quadrature, face_no, subface_no, points, ref_case);
1316 return Quadrature<dim>(points, quadrature.get_weights());
1317}
1318
1319
1320// explicit instantiations; note: we need them all for all dimensions
1321template class QProjector<1>;
1322template class QProjector<2>;
1323template class QProjector<3>;
1324
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:340
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
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)