Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
qprojector.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 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
24
25namespace internal
26{
27 namespace QProjector
28 {
29 namespace
30 {
32 reflect(const Quadrature<2> &q)
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
39 return Quadrature<2>(q_points, q.get_weights());
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
74 return Quadrature<2>(q_points, q.get_weights());
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{
88 project_to_face(ReferenceCells::Line, quadrature, face_no, q_points);
89}
90
91
92
93template <>
94void
96 const Quadrature<0> &,
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{
118 project_to_face(ReferenceCells::Quadrilateral, quadrature, face_no, q_points);
119}
120
121
122
123template <>
124void
126 const Quadrature<1> & quadrature,
127 const unsigned int face_no,
128 std::vector<Point<2>> &q_points)
129{
130 const unsigned int dim = 2;
132 Assert(q_points.size() == quadrature.size(),
133 ExcDimensionMismatch(q_points.size(), quadrature.size()));
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:
143 q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
144 break;
145 case 1:
146 q_points[p] =
147 Point<dim>(1 - quadrature.point(p)(0), quadrature.point(p)(0));
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:
162 q_points[p] = Point<dim>(0, quadrature.point(p)(0));
163 break;
164 case 1:
165 q_points[p] = Point<dim>(1, quadrature.point(p)(0));
166 break;
167 case 2:
168 q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
169 break;
170 case 3:
171 q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
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{
191 project_to_face(ReferenceCells::Hexahedron, quadrature, face_no, q_points);
192}
193
194
195
196template <>
197void
199 const Quadrature<2> & quadrature,
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;
208 Assert(q_points.size() == quadrature.size(),
209 ExcDimensionMismatch(q_points.size(), quadrature.size()));
210
211 for (unsigned int p = 0; p < quadrature.size(); ++p)
212 switch (face_no)
213 {
214 case 0:
215 q_points[p] =
216 Point<dim>(0, quadrature.point(p)(0), quadrature.point(p)(1));
217 break;
218 case 1:
219 q_points[p] =
220 Point<dim>(1, quadrature.point(p)(0), quadrature.point(p)(1));
221 break;
222 case 2:
223 q_points[p] =
224 Point<dim>(quadrature.point(p)(1), 0, quadrature.point(p)(0));
225 break;
226 case 3:
227 q_points[p] =
228 Point<dim>(quadrature.point(p)(1), 1, quadrature.point(p)(0));
229 break;
230 case 4:
231 q_points[p] =
232 Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 0);
233 break;
234 case 5:
235 q_points[p] =
236 Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 1);
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
263 const Quadrature<0> &,
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{
289 project_to_subface(ReferenceCells::Quadrilateral,
290 quadrature,
291 face_no,
292 subface_no,
293 q_points,
294 ref_case);
295}
296
297
298
299template <>
300void
302 const Quadrature<1> & quadrature,
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
312 Assert(q_points.size() == quadrature.size(),
313 ExcDimensionMismatch(q_points.size(), quadrature.size()));
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,
341 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,
452 quadrature,
453 face_no,
454 subface_no,
455 q_points,
456 ref_case);
457}
458
459
460
461template <>
462void
464 const Quadrature<2> & quadrature,
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;
476 Assert(q_points.size() == quadrature.size(),
477 ExcDimensionMismatch(q_points.size(), quadrature.size()));
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,
539 // projected quadrature points
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 <>
554 const hp::QCollection<0> &quadrature)
555{
556 AssertDimension(quadrature.size(), 1);
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 {
574 project_to_face(quadrature[quadrature.size() == 1 ? 0 : face],
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(
585 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().begin(),
586 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().end(),
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
592 return Quadrature<dim>(q_points, weights);
593}
594
595
596
597template <>
600 const hp::QCollection<1> &quadrature)
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 ...
642 const auto &sub_quadrature_points =
643 quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
644 const auto &sub_quadrature_weights =
645 quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
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] *
656 poly.compute_value(i, sub_quadrature_points[j]);
657
658 points.emplace_back(mapped_point);
659
660 // scale weight by arc length
661 weights.emplace_back(sub_quadrature_weights[j] * face.second);
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
677 if (quadrature.size() == 1)
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)
683 n_points_total += quadrature[i].size();
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;
690 help.reserve(quadrature.max_n_quadrature_points());
691
692 // project to each face and append
693 // results
694 for (unsigned int face = 0; face < n_faces; ++face)
695 {
696 help.resize(quadrature[quadrature.size() == 1 ? 0 : face].size());
697 project_to_face(quadrature[quadrature.size() == 1 ? 0 : face],
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(
708 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().begin(),
709 quadrature[quadrature.size() == 1 ? 0 : face].get_weights().end(),
710 std::back_inserter(weights));
711
712 Assert(q_points.size() == n_points_total, ExcInternalError());
713 Assert(weights.size() == n_points_total, ExcInternalError());
714
715 return Quadrature<dim>(q_points, weights);
716}
717
718
719
720template <>
723 const hp::QCollection<2> &quadrature)
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
736 const auto support_points_quad =
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);
753 const TensorProductPolynomials<2> poly_quad(
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) :
767 static_cast<const ScalarPolynomialsBase<2> &>(poly_quad);
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) :
778 support_points_quad(face, orientation);
779
780 // the quadrature rule to be projected ...
781 const auto &sub_quadrature_points =
782 quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
783 const auto &sub_quadrature_weights =
784 quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
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] *
795 poly.compute_value(i, sub_quadrature_points[j]);
796
797 points.push_back(mapped_point);
798
799 // scale quadrature weight
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] =
812 poly.compute_1st_derivative(i, sub_quadrature_points[j]);
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
843 weights.push_back(sub_quadrature_weights[j] * scaling);
844 }
845 }
846 }
847
848 // construct new quadrature rule
849 return Quadrature<3>(points, weights);
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
940 if (quadrature.size() == 1)
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)
946 n_points_total += quadrature[i].size();
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;
955 help.reserve(quadrature.max_n_quadrature_points());
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 {
971 SubQuadrature mutation;
972
973 const auto quadrature_f =
974 quadrature[quadrature.size() == 1 ? 0 : face];
975 switch (i)
976 {
977 case 0:
978 mutation = quadrature_f;
979 break;
980 case 1:
981 mutation = internal::QProjector::rotate(quadrature_f, 1);
982 break;
983 case 2:
984 mutation = internal::QProjector::rotate(quadrature_f, 2);
985 break;
986 case 3:
987 mutation = internal::QProjector::rotate(quadrature_f, 3);
988 break;
989 case 4:
990 mutation = internal::QProjector::reflect(quadrature_f);
991 break;
992 case 5:
994 internal::QProjector::reflect(quadrature_f), 3);
995 break;
996 case 6:
998 internal::QProjector::reflect(quadrature_f), 2);
999 break;
1000 case 7:
1002 internal::QProjector::reflect(quadrature_f), 1);
1003 break;
1004 default:
1005 Assert(false, ExcInternalError())
1006 }
1007
1008 help.resize(quadrature[quadrature.size() == 1 ? 0 : face].size());
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
1022 return Quadrature<dim>(q_points, weights);
1023}
1024
1025
1026
1027template <>
1030{
1031 return project_to_all_subfaces(ReferenceCells::Line, quadrature);
1032}
1033
1034
1035
1036template <>
1039 const Quadrature<0> &quadrature)
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 {
1060 project_to_subface(quadrature, face, subface, help);
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)
1069 std::copy(quadrature.get_weights().begin(),
1070 quadrature.get_weights().end(),
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
1078 return Quadrature<dim>(q_points, weights);
1079}
1080
1081
1082
1083template <>
1086 const SubQuadrature &quadrature)
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 {
1111 project_to_subface(quadrature, face, subface, help);
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)
1120 std::copy(quadrature.get_weights().begin(),
1121 quadrature.get_weights().end(),
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
1129 return Quadrature<dim>(q_points, weights);
1130}
1131
1132
1133
1134template <>
1137{
1138 return project_to_all_subfaces(ReferenceCells::Quadrilateral, quadrature);
1139}
1140
1141
1142
1143template <>
1146 const SubQuadrature &quadrature)
1147{
1150 return Quadrature<3>(); // nothing to do
1151
1153
1154 const unsigned int dim = 3;
1155 SubQuadrature q_reflected = internal::QProjector::reflect(quadrature);
1156 SubQuadrature q[8] = {quadrature,
1157 internal::QProjector::rotate(quadrature, 1),
1158 internal::QProjector::rotate(quadrature, 2),
1159 internal::QProjector::rotate(quadrature, 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
1221 return Quadrature<dim>(q_points, weights);
1222}
1223
1224
1225
1226template <>
1229{
1230 return project_to_all_subfaces(ReferenceCells::Hexahedron, quadrature);
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>(),
1242 quadrature,
1243 child_no);
1244}
1245
1246
1247
1248template <int dim>
1251 const Quadrature<dim> &quadrature,
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
1271 std::vector<double> weights = quadrature.get_weights();
1272 for (unsigned int i = 0; i < n_q_points; ++i)
1273 weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1274
1275 return Quadrature<dim>(q_points, weights);
1276}
1277
1278
1279
1280template <int dim>
1283{
1284 return project_to_all_children(ReferenceCells::get_hypercube<dim>(),
1285 quadrature);
1286}
1287
1288
1289
1290template <int dim>
1293 const Quadrature<dim> &quadrature)
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 {
1309 Quadrature<dim> help = project_to_child(quadrature, child);
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 }
1316 return Quadrature<dim>(q_points, weights);
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>(),
1328 quadrature,
1329 p1,
1330 p2);
1331}
1332
1333
1334
1335template <int dim>
1338 const Quadrature<1> &quadrature,
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 }
1358 return Quadrature<dim>(points, weights);
1359}
1360
1361
1362
1363template <int dim>
1366 const bool face_orientation,
1367 const bool face_flip,
1368 const bool face_rotation,
1369 const unsigned int n_quadrature_points)
1370{
1371 return face(ReferenceCells::get_hypercube<dim>(),
1372 face_no,
1373 face_orientation,
1374 face_flip,
1375 face_rotation,
1376 n_quadrature_points);
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,
1388 const unsigned int n_quadrature_points)
1389{
1392 {
1393 if (dim == 2)
1394 return {(2 * face_no + face_orientation) * n_quadrature_points};
1395 else if (dim == 3)
1396 {
1397 const unsigned int orientation =
1398 (face_flip * 2 + face_rotation) * 2 + face_orientation;
1399 return {(6 * face_no + orientation) * n_quadrature_points};
1400 }
1401 }
1402
1403 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1405
1407
1408 switch (dim)
1409 {
1410 case 1:
1411 case 2:
1412 return face_no * n_quadrature_points;
1413
1414
1415 case 3:
1416 {
1417 // in 3d, we have to account for faces that
1418 // have non-standard face orientation, flip
1419 // and rotation. thus, we have to store
1420 // _eight_ data sets per face or subface
1421
1422 // set up a table with the according offsets
1423 // for non-standard orientation, first index:
1424 // face_orientation (standard true=1), second
1425 // index: face_flip (standard false=0), third
1426 // index: face_rotation (standard false=0)
1427 //
1428 // note, that normally we should use the
1429 // obvious offsets 0,1,2,3,4,5,6,7. However,
1430 // prior to the changes enabling flipped and
1431 // rotated faces, in many places of the
1432 // library the convention was used, that the
1433 // first dataset with offset 0 corresponds to
1434 // a face in standard orientation. therefore
1435 // we use the offsets 4,5,6,7,0,1,2,3 here to
1436 // stick to that (implicit) convention
1437 static const unsigned int offset[2][2][2] = {
1439 5 * GeometryInfo<dim>::
1440 faces_per_cell}, // face_orientation=false; face_flip=false;
1441 // face_rotation=false and true
1443 7 * GeometryInfo<dim>::
1444 faces_per_cell}}, // face_orientation=false; face_flip=true;
1445 // face_rotation=false and true
1447 1 * GeometryInfo<dim>::
1448 faces_per_cell}, // face_orientation=true; face_flip=false;
1449 // face_rotation=false and true
1451 3 * GeometryInfo<dim>::
1452 faces_per_cell}}}; // face_orientation=true; face_flip=true;
1453 // face_rotation=false and true
1454
1455 return (
1456 (face_no + offset[face_orientation][face_flip][face_rotation]) *
1457 n_quadrature_points);
1458 }
1459
1460 default:
1461 Assert(false, ExcInternalError());
1462 }
1464}
1465
1466
1467
1468template <int dim>
1472 const unsigned int face_no,
1473 const bool face_orientation,
1474 const bool face_flip,
1475 const bool face_rotation,
1476 const hp::QCollection<dim - 1> &quadrature)
1477{
1482 {
1483 unsigned int offset = 0;
1484
1485 static const unsigned int X = numbers::invalid_unsigned_int;
1486 static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1487 static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1488 static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1489 static const std::array<unsigned int, 5> scale_pyramid = {
1490 {8, 6, 6, 6, 6}};
1491
1492 const auto &scale =
1494 scale_tri :
1496 scale_tet :
1497 ((reference_cell == ReferenceCells::Wedge) ? scale_wedge :
1498 scale_pyramid));
1499
1500 if (quadrature.size() == 1)
1501 offset = scale[0] * quadrature[0].size() * face_no;
1502 else
1503 for (unsigned int i = 0; i < face_no; ++i)
1504 offset += scale[i] * quadrature[i].size();
1505
1506 if (dim == 2)
1507 return {offset +
1508 face_orientation *
1509 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1510 else if (dim == 3)
1511 {
1512 const unsigned int orientation =
1513 (face_flip * 2 + face_rotation) * 2 + face_orientation;
1514
1515 return {offset +
1516 orientation *
1517 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1518 }
1519 }
1520
1521 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1523
1525
1526 switch (dim)
1527 {
1528 case 1:
1529 case 2:
1530 {
1531 if (quadrature.size() == 1)
1532 return quadrature[0].size() * face_no;
1533 else
1534 {
1535 unsigned int result = 0;
1536 for (unsigned int i = 0; i < face_no; ++i)
1537 result += quadrature[i].size();
1538 return result;
1539 }
1540 }
1541 case 3:
1542 {
1543 // in 3d, we have to account for faces that
1544 // have non-standard face orientation, flip
1545 // and rotation. thus, we have to store
1546 // _eight_ data sets per face or subface
1547
1548 // set up a table with the according offsets
1549 // for non-standard orientation, first index:
1550 // face_orientation (standard true=1), second
1551 // index: face_flip (standard false=0), third
1552 // index: face_rotation (standard false=0)
1553 //
1554 // note, that normally we should use the
1555 // obvious offsets 0,1,2,3,4,5,6,7. However,
1556 // prior to the changes enabling flipped and
1557 // rotated faces, in many places of the
1558 // library the convention was used, that the
1559 // first dataset with offset 0 corresponds to
1560 // a face in standard orientation. therefore
1561 // we use the offsets 4,5,6,7,0,1,2,3 here to
1562 // stick to that (implicit) convention
1563 static const unsigned int offset[2][2][2] = {
1564 {{4, 5}, // face_orientation=false; face_flip=false;
1565 // face_rotation=false and true
1566 {6, 7}}, // face_orientation=false; face_flip=true;
1567 // face_rotation=false and true
1568 {{0, 1}, // face_orientation=true; face_flip=false;
1569 // face_rotation=false and true
1570 {2, 3}}}; // face_orientation=true; face_flip=true;
1571 // face_rotation=false and true
1572
1573
1574 if (quadrature.size() == 1)
1575 return (face_no +
1576 offset[face_orientation][face_flip][face_rotation] *
1578 quadrature[0].size();
1579 else
1580 {
1581 unsigned int n_points_i = 0;
1582 for (unsigned int i = 0; i < face_no; ++i)
1583 n_points_i += quadrature[i].size();
1584
1585 unsigned int n_points = 0;
1586 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1587 ++i)
1588 n_points += quadrature[i].size();
1589
1590 return (n_points_i +
1591 offset[face_orientation][face_flip][face_rotation] *
1592 n_points);
1593 }
1594 }
1595
1596 default:
1597 Assert(false, ExcInternalError());
1598 }
1600}
1601
1602
1603
1604template <>
1608 const unsigned int face_no,
1609 const unsigned int subface_no,
1610 const bool,
1611 const bool,
1612 const bool,
1613 const unsigned int n_quadrature_points,
1615{
1617 (void)reference_cell;
1618
1622
1623 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1624 n_quadrature_points);
1625}
1626
1627
1628
1629template <>
1632 const unsigned int face_no,
1633 const unsigned int subface_no,
1634 const bool face_orientation,
1635 const bool face_flip,
1636 const bool face_rotation,
1637 const unsigned int n_quadrature_points,
1638 const internal::SubfaceCase<1> ref_case)
1639{
1640 return subface(ReferenceCells::Line,
1641 face_no,
1642 subface_no,
1643 face_orientation,
1644 face_flip,
1645 face_rotation,
1646 n_quadrature_points,
1647 ref_case);
1648}
1649
1650
1651
1652template <>
1656 const unsigned int face_no,
1657 const unsigned int subface_no,
1658 const bool,
1659 const bool,
1660 const bool,
1661 const unsigned int n_quadrature_points,
1663{
1665 (void)reference_cell;
1666
1670
1671 return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1672 n_quadrature_points);
1673}
1674
1675
1676
1677template <>
1680 const unsigned int face_no,
1681 const unsigned int subface_no,
1682 const bool face_orientation,
1683 const bool face_flip,
1684 const bool face_rotation,
1685 const unsigned int n_quadrature_points,
1686 const internal::SubfaceCase<2> ref_case)
1687{
1688 return subface(ReferenceCells::Quadrilateral,
1689 face_no,
1690 subface_no,
1691 face_orientation,
1692 face_flip,
1693 face_rotation,
1694 n_quadrature_points,
1695 ref_case);
1696}
1697
1698
1699template <>
1703 const unsigned int face_no,
1704 const unsigned int subface_no,
1705 const bool face_orientation,
1706 const bool face_flip,
1707 const bool face_rotation,
1708 const unsigned int n_quadrature_points,
1709 const internal::SubfaceCase<3> ref_case)
1710{
1711 const unsigned int dim = 3;
1712
1714 (void)reference_cell;
1715
1719
1720 // As the quadrature points created by
1721 // QProjector are on subfaces in their
1722 // "standard location" we have to use a
1723 // permutation of the equivalent subface
1724 // number in order to respect face
1725 // orientation, flip and rotation. The
1726 // information we need here is exactly the
1727 // same as the
1728 // GeometryInfo<3>::child_cell_on_face info
1729 // for the bottom face (face 4) of a hex, as
1730 // on this the RefineCase of the cell matches
1731 // that of the face and the subfaces are
1732 // numbered in the same way as the child
1733 // cells.
1734
1735 // in 3d, we have to account for faces that
1736 // have non-standard face orientation, flip
1737 // and rotation. thus, we have to store
1738 // _eight_ data sets per face or subface
1739 // already for the isotropic
1740 // case. Additionally, we have three
1741 // different refinement cases, resulting in
1742 // <tt>4 + 2 + 2 = 8</tt> different subfaces
1743 // for each face.
1744 const unsigned int total_subfaces_per_face = 8;
1745
1746 // set up a table with the according offsets
1747 // for non-standard orientation, first index:
1748 // face_orientation (standard true=1), second
1749 // index: face_flip (standard false=0), third
1750 // index: face_rotation (standard false=0)
1751 //
1752 // note, that normally we should use the
1753 // obvious offsets 0,1,2,3,4,5,6,7. However,
1754 // prior to the changes enabling flipped and
1755 // rotated faces, in many places of the
1756 // library the convention was used, that the
1757 // first dataset with offset 0 corresponds to
1758 // a face in standard orientation. therefore
1759 // we use the offsets 4,5,6,7,0,1,2,3 here to
1760 // stick to that (implicit) convention
1761 static const unsigned int orientation_offset[2][2][2] = {
1762 {// face_orientation=false; face_flip=false; face_rotation=false and true
1763 {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1764 5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1765 // face_orientation=false; face_flip=true; face_rotation=false and true
1766 {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1767 7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
1768 {// face_orientation=true; face_flip=false; face_rotation=false and true
1769 {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1770 1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1771 // face_orientation=true; face_flip=true; face_rotation=false and true
1772 {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1773 3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};
1774
1775 // set up a table with the offsets for a
1776 // given refinement case respecting the
1777 // corresponding number of subfaces. the
1778 // index corresponds to (RefineCase::Type - 1)
1779
1780 // note, that normally we should use the
1781 // obvious offsets 0,2,6. However, prior to
1782 // the implementation of anisotropic
1783 // refinement, in many places of the library
1784 // the convention was used, that the first
1785 // dataset with offset 0 corresponds to a
1786 // standard (isotropic) face
1787 // refinement. therefore we use the offsets
1788 // 6,4,0 here to stick to that (implicit)
1789 // convention
1790 static const unsigned int ref_case_offset[3] = {
1791 6, // cut_x
1792 4, // cut_y
1793 0 // cut_xy
1794 };
1795
1796
1797 // for each subface of a given FaceRefineCase
1798 // there is a corresponding equivalent
1799 // subface number of one of the "standard"
1800 // RefineCases (cut_x, cut_y, cut_xy). Map
1801 // the given values to those equivalent
1802 // ones.
1803
1804 // first, define an invalid number
1805 static const unsigned int e = numbers::invalid_unsigned_int;
1806
1807 static const RefinementCase<dim - 1>
1808 equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
1810 // case_none. there should be only
1811 // invalid values here. However, as
1812 // this function is also called (in
1813 // tests) for cells which have no
1814 // refined faces, use isotropic
1815 // refinement instead
1816 {RefinementCase<dim - 1>::cut_xy,
1817 RefinementCase<dim - 1>::cut_xy,
1818 RefinementCase<dim - 1>::cut_xy,
1819 RefinementCase<dim - 1>::cut_xy},
1820 // case_x
1821 {RefinementCase<dim - 1>::cut_x,
1822 RefinementCase<dim - 1>::cut_x,
1823 RefinementCase<dim - 1>::no_refinement,
1824 RefinementCase<dim - 1>::no_refinement},
1825 // case_x1y
1826 {RefinementCase<dim - 1>::cut_xy,
1827 RefinementCase<dim - 1>::cut_xy,
1828 RefinementCase<dim - 1>::cut_x,
1829 RefinementCase<dim - 1>::no_refinement},
1830 // case_x2y
1831 {RefinementCase<dim - 1>::cut_x,
1832 RefinementCase<dim - 1>::cut_xy,
1833 RefinementCase<dim - 1>::cut_xy,
1834 RefinementCase<dim - 1>::no_refinement},
1835 // case_x1y2y
1836 {RefinementCase<dim - 1>::cut_xy,
1837 RefinementCase<dim - 1>::cut_xy,
1838 RefinementCase<dim - 1>::cut_xy,
1839 RefinementCase<dim - 1>::cut_xy},
1840 // case_y
1841 {RefinementCase<dim - 1>::cut_y,
1842 RefinementCase<dim - 1>::cut_y,
1843 RefinementCase<dim - 1>::no_refinement,
1844 RefinementCase<dim - 1>::no_refinement},
1845 // case_y1x
1846 {RefinementCase<dim - 1>::cut_xy,
1847 RefinementCase<dim - 1>::cut_xy,
1848 RefinementCase<dim - 1>::cut_y,
1849 RefinementCase<dim - 1>::no_refinement},
1850 // case_y2x
1851 {RefinementCase<dim - 1>::cut_y,
1852 RefinementCase<dim - 1>::cut_xy,
1853 RefinementCase<dim - 1>::cut_xy,
1854 RefinementCase<dim - 1>::no_refinement},
1855 // case_y1x2x
1856 {RefinementCase<dim - 1>::cut_xy,
1857 RefinementCase<dim - 1>::cut_xy,
1858 RefinementCase<dim - 1>::cut_xy,
1859 RefinementCase<dim - 1>::cut_xy},
1860 // case_xy (case_isotropic)
1861 {RefinementCase<dim - 1>::cut_xy,
1862 RefinementCase<dim - 1>::cut_xy,
1863 RefinementCase<dim - 1>::cut_xy,
1864 RefinementCase<dim - 1>::cut_xy}};
1865
1866 static const unsigned int
1867 equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
1869 // case_none, see above
1870 {0, 1, 2, 3},
1871 // case_x
1872 {0, 1, e, e},
1873 // case_x1y
1874 {0, 2, 1, e},
1875 // case_x2y
1876 {0, 1, 3, e},
1877 // case_x1y2y
1878 {0, 2, 1, 3},
1879 // case_y
1880 {0, 1, e, e},
1881 // case_y1x
1882 {0, 1, 1, e},
1883 // case_y2x
1884 {0, 2, 3, e},
1885 // case_y1x2x
1886 {0, 1, 2, 3},
1887 // case_xy (case_isotropic)
1888 {0, 1, 2, 3}};
1889
1890 // If face-orientation or face_rotation are
1891 // non-standard, cut_x and cut_y have to be
1892 // exchanged.
1893 static const RefinementCase<dim - 1> ref_case_permutation[4] = {
1894 RefinementCase<dim - 1>::no_refinement,
1895 RefinementCase<dim - 1>::cut_y,
1896 RefinementCase<dim - 1>::cut_x,
1897 RefinementCase<dim - 1>::cut_xy};
1898
1899 // set a corresponding (equivalent)
1900 // RefineCase and subface number
1901 const RefinementCase<dim - 1> equ_ref_case =
1902 equivalent_refine_case[ref_case][subface_no];
1903 const unsigned int equ_subface_no =
1904 equivalent_subface_number[ref_case][subface_no];
1905 // make sure, that we got a valid subface and RefineCase
1908 Assert(equ_subface_no != e, ExcInternalError());
1909 // now, finally respect non-standard faces
1910 const RefinementCase<dim - 1> final_ref_case =
1911 (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1912 equ_ref_case);
1913
1914 // what we have now is the number of
1915 // the subface in the natural
1916 // orientation of the *face*. what we
1917 // need to know is the number of the
1918 // subface concerning the standard face
1919 // orientation as seen from the *cell*.
1920
1921 // this mapping is not trivial, but we
1922 // have done exactly this stuff in the
1923 // child_cell_on_face function. in
1924 // order to reduce the amount of code
1925 // as well as to make maintaining the
1926 // functionality easier we want to
1927 // reuse that information. So we note
1928 // that on the bottom face (face 4) of
1929 // a hex cell the local x and y
1930 // coordinates of the face and the cell
1931 // coincide, thus also the refinement
1932 // case of the face corresponds to the
1933 // refinement case of the cell
1934 // (ignoring cell refinement along the
1935 // z direction). Using this knowledge
1936 // we can (ab)use the
1937 // child_cell_on_face function to do
1938 // exactly the transformation we are in
1939 // need of now
1940 const unsigned int final_subface_no =
1942 4,
1943 equ_subface_no,
1944 face_orientation,
1945 face_flip,
1946 face_rotation,
1947 equ_ref_case);
1948
1949 return (((face_no * total_subfaces_per_face +
1950 ref_case_offset[final_ref_case - 1] + final_subface_no) +
1951 orientation_offset[face_orientation][face_flip][face_rotation]) *
1952 n_quadrature_points);
1953}
1954
1955
1956template <>
1959 const unsigned int face_no,
1960 const unsigned int subface_no,
1961 const bool face_orientation,
1962 const bool face_flip,
1963 const bool face_rotation,
1964 const unsigned int n_quadrature_points,
1965 const internal::SubfaceCase<3> ref_case)
1966{
1967 return subface(ReferenceCells::Hexahedron,
1968 face_no,
1969 subface_no,
1970 face_orientation,
1971 face_flip,
1972 face_rotation,
1973 n_quadrature_points,
1974 ref_case);
1975}
1976
1977
1978
1979template <int dim>
1982 const unsigned int face_no)
1983{
1984 return project_to_face(ReferenceCells::get_hypercube<dim>(),
1985 quadrature,
1986 face_no);
1987}
1988
1989
1990
1991template <int dim>
1994 const SubQuadrature &quadrature,
1995 const unsigned int face_no)
1996{
1997 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1999 (void)reference_cell;
2000
2001 std::vector<Point<dim>> points(quadrature.size());
2002 project_to_face(quadrature, face_no, points);
2003 return Quadrature<dim>(points, quadrature.get_weights());
2004}
2005
2006
2007
2008template <int dim>
2011 const unsigned int face_no,
2012 const unsigned int subface_no,
2013 const RefinementCase<dim - 1> &ref_case)
2014{
2015 return project_to_subface(ReferenceCells::get_hypercube<dim>(),
2016 quadrature,
2017 face_no,
2018 subface_no,
2019 ref_case);
2020}
2021
2022
2023
2024template <int dim>
2027 const SubQuadrature &quadrature,
2028 const unsigned int face_no,
2029 const unsigned int subface_no,
2030 const RefinementCase<dim - 1> &ref_case)
2031{
2032 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
2034 (void)reference_cell;
2035
2036 std::vector<Point<dim>> points(quadrature.size());
2037 project_to_subface(quadrature, face_no, subface_no, points, ref_case);
2038 return Quadrature<dim>(points, quadrature.get_weights());
2039}
2040
2041
2042// explicit instantiations; note: we need them all for all dimensions
2043template class QProjector<1>;
2044template class QProjector<2>;
2045template class QProjector<3>;
2046
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)
static Quadrature< dim > project_to_all_faces(const Quadrature< dim - 1 > &quadrature)
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
static Quadrature< dim > project_to_all_children(const Quadrature< dim > &quadrature)
Definition: qprojector.cc:1282
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
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)
unsigned int size() const
Definition: collection.h:109
unsigned int max_n_quadrature_points() const
Definition: q_collection.h:181
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:2042
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:701
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
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
::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)