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