Reference documentation for deal.II version 9.4.0
qprojector.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2022 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
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
24
25namespace internal
26{
27 namespace QProjector
28 {
29 namespace
30 {
33 {
34 // Take the points and reflect them by the diagonal
35 std::vector<Point<2>> q_points(q.get_points());
36 for (Point<2> &p : q_points)
37 std::swap(p[0], p[1]);
38
40 }
41
42
44 rotate(const Quadrature<2> &q, const unsigned int n_times)
45 {
46 std::vector<Point<2>> q_points(q.size());
47 for (unsigned int i = 0; i < q.size(); ++i)
48 {
49 switch (n_times % 4)
50 {
51 case 0:
52 // 0 degree. the point remains as it is.
53 q_points[i] = q.point(i);
54 break;
55
56 case 1:
57 // 90 degree counterclockwise
58 q_points[i][0] = 1.0 - q.point(i)[1];
59 q_points[i][1] = q.point(i)[0];
60 break;
61 case 2:
62 // 180 degree counterclockwise
63 q_points[i][0] = 1.0 - q.point(i)[0];
64 q_points[i][1] = 1.0 - q.point(i)[1];
65 break;
66 case 3:
67 // 270 degree counterclockwise
68 q_points[i][0] = q.point(i)[1];
69 q_points[i][1] = 1.0 - q.point(i)[0];
70 break;
71 }
72 }
73
75 }
76 } // namespace
77 } // namespace QProjector
78} // namespace internal
79
80
81
82template <>
83void
85 const unsigned int face_no,
86 std::vector<Point<1>> &q_points)
87{
89}
90
91
92
93template <>
94void
97 const unsigned int face_no,
98 std::vector<Point<1>> &q_points)
99{
101 (void)reference_cell;
102
103 const unsigned int dim = 1;
105 AssertDimension(q_points.size(), 1);
106
107 q_points[0] = Point<dim>(static_cast<double>(face_no));
108}
109
110
111
112template <>
113void
115 const unsigned int face_no,
116 std::vector<Point<2>> &q_points)
117{
119}
120
121
122
123template <>
124void
127 const unsigned int face_no,
128 std::vector<Point<2>> &q_points)
129{
130 const unsigned int dim = 2;
134
136 {
137 // use linear polynomial to map the reference quadrature points correctly
138 // on faces, i.e., BarycentricPolynomials<1>(1)
139 for (unsigned int p = 0; p < quadrature.size(); ++p)
140 switch (face_no)
141 {
142 case 0:
144 break;
145 case 1:
146 q_points[p] =
148 break;
149 case 2:
150 q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0));
151 break;
152 default:
153 Assert(false, ExcInternalError());
154 }
155 }
157 {
158 for (unsigned int p = 0; p < quadrature.size(); ++p)
159 switch (face_no)
160 {
161 case 0:
163 break;
164 case 1:
166 break;
167 case 2:
169 break;
170 case 3:
172 break;
173 default:
174 Assert(false, ExcInternalError());
175 }
176 }
177 else
178 {
179 Assert(false, ExcInternalError());
180 }
181}
182
183
184
185template <>
186void
188 const unsigned int face_no,
189 std::vector<Point<3>> &q_points)
190{
192}
193
194
195
196template <>
197void
200 const unsigned int face_no,
201 std::vector<Point<3>> &q_points)
202{
204 (void)reference_cell;
205
206 const unsigned int dim = 3;
210
211 for (unsigned int p = 0; p < quadrature.size(); ++p)
212 switch (face_no)
213 {
214 case 0:
215 q_points[p] =
217 break;
218 case 1:
219 q_points[p] =
221 break;
222 case 2:
223 q_points[p] =
225 break;
226 case 3:
227 q_points[p] =
229 break;
230 case 4:
231 q_points[p] =
233 break;
234 case 5:
235 q_points[p] =
237 break;
238
239 default:
240 Assert(false, ExcInternalError());
241 }
242}
243
244
245
246template <>
247void
249 const unsigned int face_no,
250 const unsigned int subface_no,
251 std::vector<Point<1>> & q_points,
252 const RefinementCase<0> &ref_case)
253{
254 project_to_subface(
255 ReferenceCells::Line, quadrature, face_no, subface_no, q_points, ref_case);
256}
257
258
259
260template <>
261void
264 const unsigned int face_no,
265 const unsigned int,
266 std::vector<Point<1>> &q_points,
267 const RefinementCase<0> &)
268{
270 (void)reference_cell;
271
272 const unsigned int dim = 1;
274 AssertDimension(q_points.size(), 1);
275
276 q_points[0] = Point<dim>(static_cast<double>(face_no));
277}
278
279
280
281template <>
282void
284 const unsigned int face_no,
285 const unsigned int subface_no,
286 std::vector<Point<2>> & q_points,
287 const RefinementCase<1> &ref_case)
288{
291 face_no,
292 subface_no,
293 q_points,
294 ref_case);
295}
296
297
298
299template <>
300void
303 const unsigned int face_no,
304 const unsigned int subface_no,
305 std::vector<Point<2>> &q_points,
306 const RefinementCase<1> &)
307{
308 const unsigned int dim = 2;
311
314
316 {
317 // use linear polynomial to map the reference quadrature points correctly
318 // on faces, i.e., BarycentricPolynomials<1>(1)
319 for (unsigned int p = 0; p < quadrature.size(); ++p)
320 switch (face_no)
321 {
322 case 0:
323 switch (subface_no)
324 {
325 case 0:
326 q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
327 break;
328 case 1:
329 q_points[p] =
330 Point<dim>(0.5 + quadrature.point(p)(0) / 2, 0);
331 break;
332 default:
333 Assert(false, ExcInternalError());
334 }
335 break;
336 case 1:
337 switch (subface_no)
338 {
339 case 0:
340 q_points[p] = Point<dim>(1 - quadrature.point(p)(0) / 2,
342 break;
343 case 1:
344 q_points[p] = Point<dim>(0.5 - quadrature.point(p)(0) / 2,
345 0.5 + quadrature.point(p)(0) / 2);
346 break;
347 default:
348 Assert(false, ExcInternalError());
349 }
350 break;
351 case 2:
352 switch (subface_no)
353 {
354 case 0:
355 q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0) / 2);
356 break;
357 case 1:
358 q_points[p] =
359 Point<dim>(0, 0.5 - quadrature.point(p)(0) / 2);
360 break;
361 default:
362 Assert(false, ExcInternalError());
363 }
364 break;
365 default:
366 Assert(false, ExcInternalError());
367 }
368 }
370 {
371 for (unsigned int p = 0; p < quadrature.size(); ++p)
372 switch (face_no)
373 {
374 case 0:
375 switch (subface_no)
376 {
377 case 0:
378 q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
379 break;
380 case 1:
381 q_points[p] =
382 Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
383 break;
384 default:
385 Assert(false, ExcInternalError());
386 }
387 break;
388 case 1:
389 switch (subface_no)
390 {
391 case 0:
392 q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
393 break;
394 case 1:
395 q_points[p] =
396 Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
397 break;
398 default:
399 Assert(false, ExcInternalError());
400 }
401 break;
402 case 2:
403 switch (subface_no)
404 {
405 case 0:
406 q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
407 break;
408 case 1:
409 q_points[p] =
410 Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
411 break;
412 default:
413 Assert(false, ExcInternalError());
414 }
415 break;
416 case 3:
417 switch (subface_no)
418 {
419 case 0:
420 q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
421 break;
422 case 1:
423 q_points[p] =
424 Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
425 break;
426 default:
427 Assert(false, ExcInternalError());
428 }
429 break;
430
431 default:
432 Assert(false, ExcInternalError());
433 }
434 }
435 else
436 {
437 Assert(false, ExcInternalError());
438 }
439}
440
441
442
443template <>
444void
446 const unsigned int face_no,
447 const unsigned int subface_no,
448 std::vector<Point<3>> & q_points,
449 const RefinementCase<2> &ref_case)
450{
451 project_to_subface(ReferenceCells::Hexahedron,
453 face_no,
454 subface_no,
455 q_points,
456 ref_case);
457}
458
459
460
461template <>
462void
465 const unsigned int face_no,
466 const unsigned int subface_no,
467 std::vector<Point<3>> & q_points,
468 const RefinementCase<2> &ref_case)
469{
471 (void)reference_cell;
472
473 const unsigned int dim = 3;
478
479 // one coordinate is at a const value. for
480 // faces 0, 2 and 4 this value is 0.0, for
481 // faces 1, 3 and 5 it is 1.0
482 double const_value = face_no % 2;
483 // local 2d coordinates are xi and eta,
484 // global 3d coordinates are x, y and
485 // z. those have to be mapped. the following
486 // indices tell, which global coordinate
487 // (0->x, 1->y, 2->z) corresponds to which
488 // local one
489 unsigned int xi_index = numbers::invalid_unsigned_int,
491 const_index = face_no / 2;
492 // the xi and eta values have to be scaled
493 // (by factor 0.5 or factor 1.0) depending on
494 // the refinement case and translated (by 0.0
495 // or 0.5) depending on the refinement case
496 // and subface_no.
497 double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
498 eta_translation = 0.0;
499 // set the index mapping between local and
500 // global coordinates
501 switch (face_no / 2)
502 {
503 case 0:
504 xi_index = 1;
505 eta_index = 2;
506 break;
507 case 1:
508 xi_index = 2;
509 eta_index = 0;
510 break;
511 case 2:
512 xi_index = 0;
513 eta_index = 1;
514 break;
515 }
516 // set the scale and translation parameter
517 // for individual subfaces
518 switch (ref_case)
519 {
520 case RefinementCase<dim - 1>::cut_x:
521 xi_scale = 0.5;
522 xi_translation = subface_no % 2 * 0.5;
523 break;
524 case RefinementCase<dim - 1>::cut_y:
525 eta_scale = 0.5;
526 eta_translation = subface_no % 2 * 0.5;
527 break;
528 case RefinementCase<dim - 1>::cut_xy:
529 xi_scale = 0.5;
530 eta_scale = 0.5;
531 xi_translation = int(subface_no % 2) * 0.5;
532 eta_translation = int(subface_no / 2) * 0.5;
533 break;
534 default:
535 Assert(false, ExcInternalError());
536 break;
537 }
538 // finally, compute the scaled, translated,
540 for (unsigned int p = 0; p < quadrature.size(); ++p)
541 {
542 q_points[p][xi_index] =
543 xi_scale * quadrature.point(p)(0) + xi_translation;
544 q_points[p][eta_index] =
545 eta_scale * quadrature.point(p)(1) + eta_translation;
546 q_points[p][const_index] = const_value;
547 }
548}
549
550
551template <>
555{
558 (void)reference_cell;
559
560 const unsigned int dim = 1;
561
562 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;
563
564 // first fix quadrature points
565 std::vector<Point<dim>> q_points;
566 q_points.reserve(n_points * n_faces);
567 std::vector<Point<dim>> help(n_points);
568
569
570 // project to each face and append
571 // results
572 for (unsigned int face = 0; face < n_faces; ++face)
573 {
575 face,
576 help);
577 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
578 }
579
580 // next copy over weights
581 std::vector<double> weights;
582 weights.reserve(n_points * n_faces);
583 for (unsigned int face = 0; face < n_faces; ++face)
584 std::copy(
587 std::back_inserter(weights));
588
589 Assert(q_points.size() == n_points * n_faces, ExcInternalError());
590 Assert(weights.size() == n_points * n_faces, ExcInternalError());
591
593}
594
595
596
597template <>
601{
603 {
604 const auto support_points_line =
605 [](const auto &face, const auto &orientation) -> std::vector<Point<2>> {
606 std::array<Point<2>, 2> vertices;
607 std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
608 const auto temp =
610 orientation);
611 return std::vector<Point<2>>(temp.begin(),
612 temp.begin() + face.first.size());
613 };
614
615 // reference faces (defined by its support points and arc length)
616 const std::array<std::pair<std::array<Point<2>, 2>, double>, 3> faces = {
617 {{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0},
618 {{{Point<2>(1.0, 0.0), Point<2>(0.0, 1.0)}}, std::sqrt(2.0)},
619 {{{Point<2>(0.0, 1.0), Point<2>(0.0, 0.0)}}, 1.0}}};
620
621 // linear polynomial to map the reference quadrature points correctly
622 // on faces
624
625 // new (projected) quadrature points and weights
626 std::vector<Point<2>> points;
627 std::vector<double> weights;
628
629 // loop over all faces (lines) ...
630 for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
631 // ... and over all possible orientations
632 for (unsigned int orientation = 0; orientation < 2; ++orientation)
633 {
634 const auto &face = faces[face_no];
635
636 // determine support point of the current line with the correct
637 // orientation
638 std::vector<Point<2>> support_points =
639 support_points_line(face, orientation);
640
641 // the quadrature rule to be projected ...
646
647 // loop over all quadrature points
648 for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
649 {
650 Point<2> mapped_point;
651
652 // map reference quadrature point
653 for (unsigned int i = 0; i < 2; ++i)
654 mapped_point +=
655 support_points[i] *
657
658 points.emplace_back(mapped_point);
659
660 // scale weight by arc length
662 }
663 }
664
665 // construct new quadrature rule
666 return {points, weights};
667 }
668
670
671 const unsigned int dim = 2;
672
673 const unsigned int n_faces = GeometryInfo<dim>::faces_per_cell;
674
675 unsigned int n_points_total = 0;
676
678 n_points_total = quadrature[0].size() * GeometryInfo<dim>::faces_per_cell;
679 else
680 {
682 for (unsigned int i = 0; i < quadrature.size(); ++i)
684 }
685
686 // first fix quadrature points
687 std::vector<Point<dim>> q_points;
688 q_points.reserve(n_points_total);
689 std::vector<Point<dim>> help;
691
692 // project to each face and append
693 // results
694 for (unsigned int face = 0; face < n_faces; ++face)
695 {
698 face,
699 help);
700 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
701 }
702
703 // next copy over weights
704 std::vector<double> weights;
705 weights.reserve(n_points_total);
706 for (unsigned int face = 0; face < n_faces; ++face)
707 std::copy(
710 std::back_inserter(weights));
711
712 Assert(q_points.size() == n_points_total, ExcInternalError());
713 Assert(weights.size() == n_points_total, ExcInternalError());
714
716}
717
718
719
720template <>
724{
725 const auto support_points_tri =
726 [](const auto &face, const auto &orientation) -> std::vector<Point<3>> {
727 std::array<Point<3>, 3> vertices;
728 std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
729 const auto temp =
731 orientation);
732 return std::vector<Point<3>>(temp.begin(),
733 temp.begin() + face.first.size());
734 };
735
737 [](const auto &face, const auto &orientation) -> std::vector<Point<3>> {
738 std::array<Point<3>, 4> vertices;
739 std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
740 const auto temp =
742 orientation);
743 return std::vector<Point<3>>(temp.begin(),
744 temp.begin() + face.first.size());
745 };
746
747 const auto process = [&](const auto &faces) {
748 // new (projected) quadrature points and weights
749 std::vector<Point<3>> points;
750 std::vector<double> weights;
751
752 const auto poly_tri = BarycentricPolynomials<2>::get_fe_p_basis(1);
755 {Point<1>(0.0), Point<1>(1.0)}));
756
757 // loop over all faces (triangles) ...
758 for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
759 {
760 // linear polynomial to map the reference quadrature points correctly
761 // on faces
762 const unsigned int n_shape_functions = faces[face_no].first.size();
763
764 const auto &poly =
765 n_shape_functions == 3 ?
766 static_cast<const ScalarPolynomialsBase<2> &>(poly_tri) :
768
769 // ... and over all possible orientations
770 for (unsigned int orientation = 0;
771 orientation < (n_shape_functions * 2);
772 ++orientation)
773 {
774 const auto &face = faces[face_no];
775
776 const auto support_points =
777 n_shape_functions == 3 ? support_points_tri(face, orientation) :
779
780 // the quadrature rule to be projected ...
785
786 // loop over all quadrature points
787 for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
788 {
789 Point<3> mapped_point;
790
791 // map reference quadrature point
792 for (unsigned int i = 0; i < n_shape_functions; ++i)
793 mapped_point +=
794 support_points[i] *
796
797 points.push_back(mapped_point);
798
800 const double scaling = [&]() {
801 const auto & supp_pts = support_points;
802 const unsigned int dim_ = 2;
803 const unsigned int spacedim = 3;
804
805 double result[spacedim][dim_];
806
807 std::vector<Tensor<1, dim_>> shape_derivatives(
808 n_shape_functions);
809
810 for (unsigned int i = 0; i < n_shape_functions; ++i)
811 shape_derivatives[i] =
813
814 for (unsigned int i = 0; i < spacedim; ++i)
815 for (unsigned int j = 0; j < dim_; ++j)
816 result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
817 for (unsigned int k = 1; k < n_shape_functions; ++k)
818 for (unsigned int i = 0; i < spacedim; ++i)
819 for (unsigned int j = 0; j < dim_; ++j)
820 result[i][j] +=
821 shape_derivatives[k][j] * supp_pts[k][i];
822
824
825 for (unsigned int i = 0; i < spacedim; ++i)
826 for (unsigned int j = 0; j < dim_; ++j)
827 contravariant[i][j] = result[i][j];
828
829
830 Tensor<1, spacedim> DX_t[dim_];
831 for (unsigned int i = 0; i < spacedim; ++i)
832 for (unsigned int j = 0; j < dim_; ++j)
833 DX_t[j][i] = contravariant[i][j];
834
836 for (unsigned int i = 0; i < dim_; ++i)
837 for (unsigned int j = 0; j < dim_; ++j)
838 G[i][j] = DX_t[i] * DX_t[j];
839
840 return std::sqrt(determinant(G));
841 }();
842
844 }
845 }
846 }
847
848 // construct new quadrature rule
850 };
851
853 {
854 // reference faces (defined by its support points and its area)
855 // note: the area is later not used as a scaling factor but recomputed
856 const std::vector<std::pair<std::vector<Point<3>>, double>> faces = {
857 {{{{Point<3>(0.0, 0.0, 0.0),
858 Point<3>(1.0, 0.0, 0.0),
859 Point<3>(0.0, 1.0, 0.0)}},
860 0.5},
861 {{{Point<3>(1.0, 0.0, 0.0),
862 Point<3>(0.0, 0.0, 0.0),
863 Point<3>(0.0, 0.0, 1.0)}},
864 0.5},
865 {{{Point<3>(0.0, 0.0, 0.0),
866 Point<3>(0.0, 1.0, 0.0),
867 Point<3>(0.0, 0.0, 1.0)}},
868 0.5},
869 {{{Point<3>(0.0, 1.0, 0.0),
870 Point<3>(1.0, 0.0, 0.0),
871 Point<3>(0.0, 0.0, 1.0)}},
872 0.5 * sqrt(3.0) /*equilateral triangle*/}}};
873
874 return process(faces);
875 }
877 {
878 const std::vector<std::pair<std::vector<Point<3>>, double>> faces = {
879 {{{{Point<3>(1.0, 0.0, 0.0),
880 Point<3>(0.0, 0.0, 0.0),
881 Point<3>(0.0, 1.0, 0.0)}},
882 0.5},
883 {{{Point<3>(0.0, 0.0, 1.0),
884 Point<3>(1.0, 0.0, 1.0),
885 Point<3>(0.0, 1.0, 1.0)}},
886 0.5},
887 {{{Point<3>(0.0, 0.0, 0.0),
888 Point<3>(1.0, 0.0, 0.0),
889 Point<3>(0.0, 0.0, 1.0),
890 Point<3>(1.0, 0.0, 1.0)}},
891 1.0},
892 {{{Point<3>(1.0, 0.0, 0.0),
893 Point<3>(0.0, 1.0, 0.0),
894 Point<3>(1.0, 0.0, 1.0),
895 Point<3>(0.0, 1.0, 1.0)}},
896 std::sqrt(2.0)},
897 {{{Point<3>(0.0, 1.0, 0.0),
898 Point<3>(0.0, 0.0, 0.0),
899 Point<3>(0.0, 1.0, 1.0),
900 Point<3>(0.0, 0.0, 1.0)}},
901 1.0}}};
902
903 return process(faces);
904 }
906 {
907 const std::vector<std::pair<std::vector<Point<3>>, double>> faces = {
908 {{{{Point<3>(-1.0, -1.0, 0.0),
909 Point<3>(+1.0, -1.0, 0.0),
910 Point<3>(-1.0, +1.0, 0.0),
911 Point<3>(+1.0, +1.0, 0.0)}},
912 4.0},
913 {{{Point<3>(-1.0, -1.0, 0.0),
914 Point<3>(-1.0, +1.0, 0.0),
915 Point<3>(+0.0, +0.0, 1.0)}},
916 std::sqrt(2.0)},
917 {{{Point<3>(+1.0, +1.0, 0.0),
918 Point<3>(+1.0, -1.0, 0.0),
919 Point<3>(+0.0, +0.0, 1.0)}},
920 std::sqrt(2.0)},
921 {{{Point<3>(+1.0, -1.0, 0.0),
922 Point<3>(-1.0, -1.0, 0.0),
923 Point<3>(+0.0, +0.0, 1.0)}},
924 std::sqrt(2.0)},
925 {{{Point<3>(-1.0, +1.0, 0.0),
926 Point<3>(+1.0, +1.0, 0.0),
927 Point<3>(+0.0, +0.0, 1.0)}},
928 std::sqrt(2.0)}}};
929
930 return process(faces);
931 }
932
933
935
936 const unsigned int dim = 3;
937
938 unsigned int n_points_total = 0;
939
941 n_points_total = quadrature[0].size() * GeometryInfo<dim>::faces_per_cell;
942 else
943 {
945 for (unsigned int i = 0; i < quadrature.size(); ++i)
947 }
948
949 n_points_total *= 8;
950
951 // first fix quadrature points
952 std::vector<Point<dim>> q_points;
953 q_points.reserve(n_points_total);
954 std::vector<Point<dim>> help;
956
957 std::vector<double> weights;
958 weights.reserve(n_points_total);
959
960 // do the following for all possible
961 // mutations of a face (mutation==0
962 // corresponds to a face with standard
963 // orientation, no flip and no rotation)
964 for (unsigned int i = 0; i < 8; ++i)
965 {
966 // project to each face and append
967 // results
968 for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
969 ++face)
970 {
972
975 switch (i)
976 {
977 case 0:
979 break;
980 case 1:
982 break;
983 case 2:
985 break;
986 case 3:
988 break;
989 case 4:
991 break;
992 case 5:
995 break;
996 case 6:
999 break;
1000 case 7:
1003 break;
1004 default:
1005 Assert(false, ExcInternalError())
1006 }
1007
1009 project_to_face(mutation, face, help);
1010 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1011
1012 std::copy(mutation.get_weights().begin(),
1013 mutation.get_weights().end(),
1014 std::back_inserter(weights));
1015 }
1016 }
1017
1018
1019 Assert(q_points.size() == n_points_total, ExcInternalError());
1020 Assert(weights.size() == n_points_total, ExcInternalError());
1021
1023}
1024
1025
1026
1027template <>
1030{
1032}
1033
1034
1035
1036template <>
1040{
1042 (void)reference_cell;
1043
1044 const unsigned int dim = 1;
1045
1046 const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
1047 subfaces_per_face =
1049
1050 // first fix quadrature points
1051 std::vector<Point<dim>> q_points;
1052 q_points.reserve(n_points * n_faces * subfaces_per_face);
1053 std::vector<Point<dim>> help(n_points);
1054
1055 // project to each face and copy
1056 // results
1057 for (unsigned int face = 0; face < n_faces; ++face)
1058 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1059 {
1061 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1062 }
1063
1064 // next copy over weights
1065 std::vector<double> weights;
1066 weights.reserve(n_points * n_faces * subfaces_per_face);
1067 for (unsigned int face = 0; face < n_faces; ++face)
1068 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1071 std::back_inserter(weights));
1072
1073 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1075 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1077
1079}
1080
1081
1082
1083template <>
1087{
1090 return Quadrature<2>(); // nothing to do
1091
1093
1094 const unsigned int dim = 2;
1095
1096 const unsigned int n_points = quadrature.size(),
1098 subfaces_per_face =
1100
1101 // first fix quadrature points
1102 std::vector<Point<dim>> q_points;
1103 q_points.reserve(n_points * n_faces * subfaces_per_face);
1104 std::vector<Point<dim>> help(n_points);
1105
1106 // project to each face and copy
1107 // results
1108 for (unsigned int face = 0; face < n_faces; ++face)
1109 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1110 {
1112 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1113 }
1114
1115 // next copy over weights
1116 std::vector<double> weights;
1117 weights.reserve(n_points * n_faces * subfaces_per_face);
1118 for (unsigned int face = 0; face < n_faces; ++face)
1119 for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
1122 std::back_inserter(weights));
1123
1124 Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
1126 Assert(weights.size() == n_points * n_faces * subfaces_per_face,
1128
1130}
1131
1132
1133
1134template <>
1137{
1139}
1140
1141
1142
1143template <>
1147{
1150 return Quadrature<3>(); // nothing to do
1151
1153
1154 const unsigned int dim = 3;
1160 q_reflected,
1161 internal::QProjector::rotate(q_reflected, 3),
1162 internal::QProjector::rotate(q_reflected, 2),
1163 internal::QProjector::rotate(q_reflected, 1)};
1164
1165 const unsigned int n_points = quadrature.size(),
1167 total_subfaces_per_face = 2 + 2 + 4;
1168
1169 // first fix quadrature points
1170 std::vector<Point<dim>> q_points;
1171 q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1172 std::vector<Point<dim>> help(n_points);
1173
1174 std::vector<double> weights;
1175 weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1176
1177 // do the following for all possible
1178 // mutations of a face (mutation==0
1179 // corresponds to a face with standard
1180 // orientation, no flip and no rotation)
1181 for (const auto &mutation : q)
1182 {
1183 // project to each face and copy
1184 // results
1185 for (unsigned int face = 0; face < n_faces; ++face)
1186 for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1187 ref_case >= RefinementCase<dim - 1>::cut_x;
1188 --ref_case)
1189 for (unsigned int subface = 0;
1191 RefinementCase<dim - 1>(ref_case));
1192 ++subface)
1193 {
1194 project_to_subface(mutation,
1195 face,
1196 subface,
1197 help,
1198 RefinementCase<dim - 1>(ref_case));
1199 std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1200 }
1201
1202 // next copy over weights
1203 for (unsigned int face = 0; face < n_faces; ++face)
1204 for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1205 ref_case >= RefinementCase<dim - 1>::cut_x;
1206 --ref_case)
1207 for (unsigned int subface = 0;
1209 RefinementCase<dim - 1>(ref_case));
1210 ++subface)
1211 std::copy(mutation.get_weights().begin(),
1212 mutation.get_weights().end(),
1213 std::back_inserter(weights));
1214 }
1215
1216 Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1218 Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1220
1222}
1223
1224
1225
1226template <>
1229{
1231}
1232
1233
1234
1235// This function is not used in the library
1236template <int dim>
1239 const unsigned int child_no)
1240{
1241 return project_to_child(ReferenceCells::get_hypercube<dim>(),
1243 child_no);
1244}
1245
1246
1247
1248template <int dim>
1252 const unsigned int child_no)
1253{
1254 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1256 (void)reference_cell;
1257
1259
1260 const unsigned int n_q_points = quadrature.size();
1261
1262 std::vector<Point<dim>> q_points(n_q_points);
1263 for (unsigned int i = 0; i < n_q_points; ++i)
1264 q_points[i] =
1266 child_no);
1267
1268 // for the weights, things are
1269 // equally simple: copy them and
1270 // scale them
1272 for (unsigned int i = 0; i < n_q_points; ++i)
1273 weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1274
1276}
1277
1278
1279
1280template <int dim>
1283{
1284 return project_to_all_children(ReferenceCells::get_hypercube<dim>(),
1286}
1287
1288
1289
1290template <int dim>
1294{
1295 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1297 (void)reference_cell;
1298
1299 const unsigned int n_points = quadrature.size(),
1301
1302 std::vector<Point<dim>> q_points(n_points * n_children);
1303 std::vector<double> weights(n_points * n_children);
1304
1305 // project to each child and copy
1306 // results
1307 for (unsigned int child = 0; child < n_children; ++child)
1308 {
1310 for (unsigned int i = 0; i < n_points; ++i)
1311 {
1312 q_points[child * n_points + i] = help.point(i);
1313 weights[child * n_points + i] = help.weight(i);
1314 }
1315 }
1317}
1318
1319
1320
1321template <int dim>
1324 const Point<dim> & p1,
1325 const Point<dim> & p2)
1326{
1327 return project_to_line(ReferenceCells::get_hypercube<dim>(),
1329 p1,
1330 p2);
1331}
1332
1333
1334
1335template <int dim>
1339 const Point<dim> & p1,
1340 const Point<dim> & p2)
1341{
1342 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1344 (void)reference_cell;
1345
1346 const unsigned int n = quadrature.size();
1347 std::vector<Point<dim>> points(n);
1348 std::vector<double> weights(n);
1349 const double length = p1.distance(p2);
1350
1351 for (unsigned int k = 0; k < n; ++k)
1352 {
1353 const double alpha = quadrature.point(k)(0);
1354 points[k] = alpha * p2;
1355 points[k] += (1. - alpha) * p1;
1356 weights[k] = length * quadrature.weight(k);
1357 }
1359}
1360
1361
1362
1363template <int dim>
1366 const bool face_orientation,
1367 const bool face_flip,
1368 const bool face_rotation,
1370{
1371 return face(ReferenceCells::get_hypercube<dim>(),
1372 face_no,
1373 face_orientation,
1374 face_flip,
1375 face_rotation,
1377}
1378
1379
1380
1381template <int dim>
1384 const unsigned int face_no,
1385 const bool face_orientation,
1386 const bool face_flip,
1387 const bool face_rotation,
1389{
1392 {
1393 if (dim == 2)
1394 return {(2 * face_no + (face_orientation ? 1 : 0)) *
1396 else if (dim == 3)
1397 {
1398 const unsigned int orientation = (face_flip ? 4 : 0) +
1399 (face_rotation ? 2 : 0) +
1400 (face_orientation ? 1 : 0);
1401 return {(6 * face_no + orientation) * n_quadrature_points};
1402 }
1403 }
1404
1405 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1407
1409
1410 switch (dim)
1411 {
1412 case 1:
1413 case 2:
1415
1416
1417 case 3:
1418 {
1419 // in 3d, we have to account for faces that
1420 // have non-standard face orientation, flip
1421 // and rotation. thus, we have to store
1422 // _eight_ data sets per face or subface
1423
1424 // set up a table with the according offsets
1425 // for non-standard orientation, first index:
1426 // face_orientation (standard true=1), second
1427 // index: face_flip (standard false=0), third
1428 // index: face_rotation (standard false=0)
1429 //
1430 // note, that normally we should use the
1431 // obvious offsets 0,1,2,3,4,5,6,7. However,
1432 // prior to the changes enabling flipped and
1433 // rotated faces, in many places of the
1434 // library the convention was used, that the
1435 // first dataset with offset 0 corresponds to
1436 // a face in standard orientation. therefore
1437 // we use the offsets 4,5,6,7,0,1,2,3 here to
1438 // stick to that (implicit) convention
1439 static const unsigned int offset[2][2][2] = {
1441 5 * GeometryInfo<dim>::
1442 faces_per_cell}, // face_orientation=false; face_flip=false;
1443 // face_rotation=false and true
1445 7 * GeometryInfo<dim>::
1446 faces_per_cell}}, // face_orientation=false; face_flip=true;
1447 // face_rotation=false and true
1449 1 * GeometryInfo<dim>::
1450 faces_per_cell}, // face_orientation=true; face_flip=false;
1451 // face_rotation=false and true
1453 3 * GeometryInfo<dim>::
1454 faces_per_cell}}}; // face_orientation=true; face_flip=true;
1455 // face_rotation=false and true
1456
1457 return (
1458 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1460 }
1461
1462 default:
1463 Assert(false, ExcInternalError());
1464 }
1466}
1467
1468
1469
1470template <int dim>
1474 const unsigned int face_no,
1475 const bool face_orientation,
1476 const bool face_flip,
1477 const bool face_rotation,
1478 const hp::QCollection<dim - 1> &quadrature)
1479{
1484 {
1485 unsigned int offset = 0;
1486
1487 static const unsigned int X = numbers::invalid_unsigned_int;
1488 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1489 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1490 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1491 static const std::array<unsigned int, 5> scale_pyramid = {
1492 {8, 6, 6, 6, 6}};
1493
1494 const auto &scale =
1496 scale_tri :
1498 scale_tet :
1499 ((reference_cell == ReferenceCells::Wedge) ? scale_wedge :
1500 scale_pyramid));
1501
1503 offset = scale[0] * quadrature[0].size() * face_no;
1504 else
1505 for (unsigned int i = 0; i < face_no; ++i)
1506 offset += scale[i] * quadrature[i].size();
1507
1508 if (dim == 2)
1509 return {offset +
1510 face_orientation *
1512 else if (dim == 3)
1513 {
1514 const unsigned int orientation = (face_flip ? 4 : 0) +
1515 (face_rotation ? 2 : 0) +
1516 (face_orientation ? 1 : 0);
1517
1518 return {offset +
1519 orientation *
1521 }
1522 }
1523
1524 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1526
1528
1529 switch (dim)
1530 {
1531 case 1:
1532 case 2:
1533 {
1536 else
1537 {
1538 unsigned int result = 0;
1539 for (unsigned int i = 0; i < face_no; ++i)
1541 return result;
1542 }
1543 }
1544 case 3:
1545 {
1546 // in 3d, we have to account for faces that
1547 // have non-standard face orientation, flip
1548 // and rotation. thus, we have to store
1549 // _eight_ data sets per face or subface
1550
1551 // set up a table with the according offsets
1552 // for non-standard orientation, first index:
1553 // face_orientation (standard true=1), second
1554 // index: face_flip (standard false=0), third
1555 // index: face_rotation (standard false=0)
1556 //
1557 // note, that normally we should use the
1558 // obvious offsets 0,1,2,3,4,5,6,7. However,
1559 // prior to the changes enabling flipped and
1560 // rotated faces, in many places of the
1561 // library the convention was used, that the
1562 // first dataset with offset 0 corresponds to
1563 // a face in standard orientation. therefore
1564 // we use the offsets 4,5,6,7,0,1,2,3 here to
1565 // stick to that (implicit) convention
1566 static const unsigned int offset[2][2][2] = {
1567 {{4, 5}, // face_orientation=false; face_flip=false;
1568 // face_rotation=false and true
1569 {6, 7}}, // face_orientation=false; face_flip=true;
1570 // face_rotation=false and true
1571 {{0, 1}, // face_orientation=true; face_flip=false;
1572 // face_rotation=false and true
1573 {2, 3}}}; // face_orientation=true; face_flip=true;
1574 // face_rotation=false and true
1575
1576
1578 return (face_no +
1579 offset[face_orientation][face_flip][face_rotation] *
1582 else
1583 {
1584 unsigned int n_points_i = 0;
1585 for (unsigned int i = 0; i < face_no; ++i)
1587
1588 unsigned int n_points = 0;
1589 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1590 ++i)
1592
1593 return (n_points_i +
1594 offset[face_orientation][face_flip][face_rotation] *
1595 n_points);
1596 }
1597 }
1598
1599 default:
1600 Assert(false, ExcInternalError());
1601 }
1603}
1604
1605
1606
1607template <>
1611 const unsigned int face_no,
1612 const unsigned int subface_no,
1613 const bool,
1614 const bool,
1615 const bool,
1618{
1620 (void)reference_cell;
1621
1625
1626 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1628}
1629
1630
1631
1632template <>
1635 const unsigned int face_no,
1636 const unsigned int subface_no,
1637 const bool face_orientation,
1638 const bool face_flip,
1639 const bool face_rotation,
1641 const internal::SubfaceCase<1> ref_case)
1642{
1643 return subface(ReferenceCells::Line,
1644 face_no,
1645 subface_no,
1646 face_orientation,
1647 face_flip,
1648 face_rotation,
1650 ref_case);
1651}
1652
1653
1654
1655template <>
1659 const unsigned int face_no,
1660 const unsigned int subface_no,
1661 const bool,
1662 const bool,
1663 const bool,
1666{
1668 (void)reference_cell;
1669
1673
1674 return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1676}
1677
1678
1679
1680template <>
1683 const unsigned int face_no,
1684 const unsigned int subface_no,
1685 const bool face_orientation,
1686 const bool face_flip,
1687 const bool face_rotation,
1689 const internal::SubfaceCase<2> ref_case)
1690{
1692 face_no,
1693 subface_no,
1694 face_orientation,
1695 face_flip,
1696 face_rotation,
1698 ref_case);
1699}
1700
1701
1702template <>
1706 const unsigned int face_no,
1707 const unsigned int subface_no,
1708 const bool face_orientation,
1709 const bool face_flip,
1710 const bool face_rotation,
1712 const internal::SubfaceCase<3> ref_case)
1713{
1714 const unsigned int dim = 3;
1715
1717 (void)reference_cell;
1718
1722
1723 // As the quadrature points created by
1724 // QProjector are on subfaces in their
1725 // "standard location" we have to use a
1726 // permutation of the equivalent subface
1727 // number in order to respect face
1728 // orientation, flip and rotation. The
1729 // information we need here is exactly the
1730 // same as the
1731 // GeometryInfo<3>::child_cell_on_face info
1732 // for the bottom face (face 4) of a hex, as
1733 // on this the RefineCase of the cell matches
1734 // that of the face and the subfaces are
1735 // numbered in the same way as the child
1736 // cells.
1737
1738 // in 3d, we have to account for faces that
1739 // have non-standard face orientation, flip
1740 // and rotation. thus, we have to store
1741 // _eight_ data sets per face or subface
1742 // already for the isotropic
1743 // case. Additionally, we have three
1744 // different refinement cases, resulting in
1745 // <tt>4 + 2 + 2 = 8</tt> different subfaces
1746 // for each face.
1747 const unsigned int total_subfaces_per_face = 8;
1748
1749 // set up a table with the according offsets
1750 // for non-standard orientation, first index:
1751 // face_orientation (standard true=1), second
1752 // index: face_flip (standard false=0), third
1753 // index: face_rotation (standard false=0)
1754 //
1755 // note, that normally we should use the
1756 // obvious offsets 0,1,2,3,4,5,6,7. However,
1757 // prior to the changes enabling flipped and
1758 // rotated faces, in many places of the
1759 // library the convention was used, that the
1760 // first dataset with offset 0 corresponds to
1761 // a face in standard orientation. therefore
1762 // we use the offsets 4,5,6,7,0,1,2,3 here to
1763 // stick to that (implicit) convention
1764 static const unsigned int orientation_offset[2][2][2] = {
1765 {// face_orientation=false; face_flip=false; face_rotation=false and true
1766 {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1767 5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1768 // face_orientation=false; face_flip=true; face_rotation=false and true
1769 {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1770 7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
1771 {// face_orientation=true; face_flip=false; face_rotation=false and true
1772 {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1773 1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1774 // face_orientation=true; face_flip=true; face_rotation=false and true
1775 {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1776 3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};
1777
1778 // set up a table with the offsets for a
1779 // given refinement case respecting the
1780 // corresponding number of subfaces. the
1781 // index corresponds to (RefineCase::Type - 1)
1782
1783 // note, that normally we should use the
1784 // obvious offsets 0,2,6. However, prior to
1785 // the implementation of anisotropic
1786 // refinement, in many places of the library
1787 // the convention was used, that the first
1788 // dataset with offset 0 corresponds to a
1789 // standard (isotropic) face
1790 // refinement. therefore we use the offsets
1791 // 6,4,0 here to stick to that (implicit)
1792 // convention
1793 static const unsigned int ref_case_offset[3] = {
1794 6, // cut_x
1795 4, // cut_y
1796 0 // cut_xy
1797 };
1798
1799
1800 // for each subface of a given FaceRefineCase
1801 // there is a corresponding equivalent
1802 // subface number of one of the "standard"
1803 // RefineCases (cut_x, cut_y, cut_xy). Map
1804 // the given values to those equivalent
1805 // ones.
1806
1807 // first, define an invalid number
1808 static const unsigned int e = numbers::invalid_unsigned_int;
1809
1810 static const RefinementCase<dim - 1>
1811 equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
1813 // case_none. there should be only
1814 // invalid values here. However, as
1815 // this function is also called (in
1816 // tests) for cells which have no
1817 // refined faces, use isotropic
1819 {RefinementCase<dim - 1>::cut_xy,
1820 RefinementCase<dim - 1>::cut_xy,
1821 RefinementCase<dim - 1>::cut_xy,
1822 RefinementCase<dim - 1>::cut_xy},
1823 // case_x
1824 {RefinementCase<dim - 1>::cut_x,
1825 RefinementCase<dim - 1>::cut_x,
1826 RefinementCase<dim - 1>::no_refinement,
1827 RefinementCase<dim - 1>::no_refinement},
1828 // case_x1y
1829 {RefinementCase<dim - 1>::cut_xy,
1830 RefinementCase<dim - 1>::cut_xy,
1831 RefinementCase<dim - 1>::cut_x,
1832 RefinementCase<dim - 1>::no_refinement},
1833 // case_x2y
1834 {RefinementCase<dim - 1>::cut_x,
1835 RefinementCase<dim - 1>::cut_xy,
1836 RefinementCase<dim - 1>::cut_xy,
1837 RefinementCase<dim - 1>::no_refinement},
1838 // case_x1y2y
1839 {RefinementCase<dim - 1>::cut_xy,
1840 RefinementCase<dim - 1>::cut_xy,
1841 RefinementCase<dim - 1>::cut_xy,
1842 RefinementCase<dim - 1>::cut_xy},
1843 // case_y
1844 {RefinementCase<dim - 1>::cut_y,
1845 RefinementCase<dim - 1>::cut_y,
1846 RefinementCase<dim - 1>::no_refinement,
1847 RefinementCase<dim - 1>::no_refinement},
1848 // case_y1x
1849 {RefinementCase<dim - 1>::cut_xy,
1850 RefinementCase<dim - 1>::cut_xy,
1851 RefinementCase<dim - 1>::cut_y,
1852 RefinementCase<dim - 1>::no_refinement},
1853 // case_y2x
1854 {RefinementCase<dim - 1>::cut_y,
1855 RefinementCase<dim - 1>::cut_xy,
1856 RefinementCase<dim - 1>::cut_xy,
1857 RefinementCase<dim - 1>::no_refinement},
1858 // case_y1x2x
1859 {RefinementCase<dim - 1>::cut_xy,
1860 RefinementCase<dim - 1>::cut_xy,
1861 RefinementCase<dim - 1>::cut_xy,
1862 RefinementCase<dim - 1>::cut_xy},
1863 // case_xy (case_isotropic)
1864 {RefinementCase<dim - 1>::cut_xy,
1865 RefinementCase<dim - 1>::cut_xy,
1866 RefinementCase<dim - 1>::cut_xy,
1867 RefinementCase<dim - 1>::cut_xy}};
1868
1869 static const unsigned int
1870 equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
1872 // case_none, see above
1873 {0, 1, 2, 3},
1874 // case_x
1875 {0, 1, e, e},
1876 // case_x1y
1877 {0, 2, 1, e},
1878 // case_x2y
1879 {0, 1, 3, e},
1880 // case_x1y2y
1881 {0, 2, 1, 3},
1882 // case_y
1883 {0, 1, e, e},
1884 // case_y1x
1885 {0, 1, 1, e},
1886 // case_y2x
1887 {0, 2, 3, e},
1888 // case_y1x2x
1889 {0, 1, 2, 3},
1890 // case_xy (case_isotropic)
1891 {0, 1, 2, 3}};
1892
1893 // If face-orientation or face_rotation are
1894 // non-standard, cut_x and cut_y have to be
1895 // exchanged.
1896 static const RefinementCase<dim - 1> ref_case_permutation[4] = {
1897 RefinementCase<dim - 1>::no_refinement,
1898 RefinementCase<dim - 1>::cut_y,
1899 RefinementCase<dim - 1>::cut_x,
1900 RefinementCase<dim - 1>::cut_xy};
1901
1902 // set a corresponding (equivalent)
1903 // RefineCase and subface number
1904 const RefinementCase<dim - 1> equ_ref_case =
1905 equivalent_refine_case[ref_case][subface_no];
1906 const unsigned int equ_subface_no =
1907 equivalent_subface_number[ref_case][subface_no];
1908 // make sure, that we got a valid subface and RefineCase
1911 Assert(equ_subface_no != e, ExcInternalError());
1912 // now, finally respect non-standard faces
1913 const RefinementCase<dim - 1> final_ref_case =
1914 (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1915 equ_ref_case);
1916
1917 // what we have now is the number of
1918 // the subface in the natural
1919 // orientation of the *face*. what we
1920 // need to know is the number of the
1921 // subface concerning the standard face
1922 // orientation as seen from the *cell*.
1923
1924 // this mapping is not trivial, but we
1925 // have done exactly this stuff in the
1926 // child_cell_on_face function. in
1927 // order to reduce the amount of code
1928 // as well as to make maintaining the
1929 // functionality easier we want to
1930 // reuse that information. So we note
1931 // that on the bottom face (face 4) of
1932 // a hex cell the local x and y
1933 // coordinates of the face and the cell
1934 // coincide, thus also the refinement
1935 // case of the face corresponds to the
1936 // refinement case of the cell
1937 // (ignoring cell refinement along the
1938 // z direction). Using this knowledge
1939 // we can (ab)use the
1940 // child_cell_on_face function to do
1941 // exactly the transformation we are in
1942 // need of now
1943 const unsigned int final_subface_no =
1945 4,
1946 equ_subface_no,
1947 face_orientation,
1948 face_flip,
1949 face_rotation,
1950 equ_ref_case);
1951
1952 return (((face_no * total_subfaces_per_face +
1953 ref_case_offset[final_ref_case - 1] + final_subface_no) +
1954 orientation_offset[face_orientation][face_flip][face_rotation]) *
1956}
1957
1958
1959template <>
1962 const unsigned int face_no,
1963 const unsigned int subface_no,
1964 const bool face_orientation,
1965 const bool face_flip,
1966 const bool face_rotation,
1968 const internal::SubfaceCase<3> ref_case)
1969{
1970 return subface(ReferenceCells::Hexahedron,
1971 face_no,
1972 subface_no,
1973 face_orientation,
1974 face_flip,
1975 face_rotation,
1977 ref_case);
1978}
1979
1980
1981
1982template <int dim>
1985 const unsigned int face_no)
1986{
1987 return project_to_face(ReferenceCells::get_hypercube<dim>(),
1989 face_no);
1990}
1991
1992
1993
1994template <int dim>
1998 const unsigned int face_no)
1999{
2000 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
2002 (void)reference_cell;
2003
2007}
2008
2009
2010
2011template <int dim>
2014 const unsigned int face_no,
2015 const unsigned int subface_no,
2016 const RefinementCase<dim - 1> &ref_case)
2017{
2018 return project_to_subface(ReferenceCells::get_hypercube<dim>(),
2020 face_no,
2021 subface_no,
2022 ref_case);
2023}
2024
2025
2026
2027template <int dim>
2031 const unsigned int face_no,
2032 const unsigned int subface_no,
2033 const RefinementCase<dim - 1> &ref_case)
2034{
2035 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
2037 (void)reference_cell;
2038
2040 project_to_subface(quadrature, face_no, subface_no, points, ref_case);
2042}
2043
2044
2045// explicit instantiations; note: we need them all for all dimensions
2046template class QProjector<1>;
2047template class QProjector<2>;
2048template class QProjector<3>;
2049
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
static DataSetDescriptor subface(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 unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: qprojector.cc:1365
static void project_to_subface(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)
Definition: qprojector.h:579
static Quadrature< dim > project_to_line(const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
Definition: qprojector.cc:1323
Definition: qprojector.cc:1282
Definition: qprojector.cc:1238
static void project_to_face(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::array< T, N > permute_according_orientation(const std::array< T, N > &vertices, const unsigned int orientation) const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
CollectionIterator< T > begin() const
Definition: collection.h:283
unsigned int size() const
Definition: collection.h:264
CollectionIterator< T > end() const
Definition: collection.h:292
Definition: q_collection.h:174
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2084
void rotate(const double angle, Triangulation< dim > &triangulation)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition: polynomial.cc:702
constexpr const ReferenceCell Tetrahedron