Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.4.1
\(\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 - 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
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{
100 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
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
135 if (reference_cell == ReferenceCells::Triangle)
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 }
156 else if (reference_cell == ReferenceCells::Quadrilateral)
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{
269 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
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
315 if (reference_cell == ReferenceCells::Triangle)
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 }
369 else if (reference_cell == ReferenceCells::Quadrilateral)
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);
557 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
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{
602 if (reference_cell == ReferenceCells::Triangle)
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
852 if (reference_cell == ReferenceCells::Tetrahedron)
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 }
876 else if (reference_cell == ReferenceCells::Wedge)
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 }
905 else if (reference_cell == ReferenceCells::Pyramid)
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:
993 mutation = internal::QProjector::rotate(
994 internal::QProjector::reflect(quadrature_f), 3);
995 break;
996 case 6:
997 mutation = internal::QProjector::rotate(
998 internal::QProjector::reflect(quadrature_f), 2);
999 break;
1000 case 7:
1001 mutation = internal::QProjector::rotate(
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{
1041 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
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{
1088 if (reference_cell == ReferenceCells::Triangle ||
1089 reference_cell == ReferenceCells::Tetrahedron)
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{
1148 if (reference_cell == ReferenceCells::Triangle ||
1149 reference_cell == ReferenceCells::Tetrahedron)
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{
1390 if (reference_cell == ReferenceCells::Triangle ||
1391 reference_cell == ReferenceCells::Tetrahedron)
1392 {
1393 if (dim == 2)
1394 return {(2 * face_no + (face_orientation ? 1 : 0)) *
1395 n_quadrature_points};
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:
1414 return face_no * n_quadrature_points;
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]) *
1459 n_quadrature_points);
1460 }
1461
1462 default:
1463 Assert(false, ExcInternalError());
1464 }
1466}
1467
1468
1469
1470template <int dim>
1473 const ReferenceCell reference_cell,
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{
1480 if (reference_cell == ReferenceCells::Triangle ||
1481 reference_cell == ReferenceCells::Tetrahedron ||
1482 reference_cell == ReferenceCells::Wedge ||
1483 reference_cell == ReferenceCells::Pyramid)
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 =
1495 (reference_cell == ReferenceCells::Triangle) ?
1496 scale_tri :
1497 ((reference_cell == ReferenceCells::Tetrahedron) ?
1498 scale_tet :
1499 ((reference_cell == ReferenceCells::Wedge) ? scale_wedge :
1500 scale_pyramid));
1501
1502 if (quadrature.size() == 1)
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 *
1511 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
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 *
1520 quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1521 }
1522 }
1523
1524 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1526
1528
1529 switch (dim)
1530 {
1531 case 1:
1532 case 2:
1533 {
1534 if (quadrature.size() == 1)
1535 return quadrature[0].size() * face_no;
1536 else
1537 {
1538 unsigned int result = 0;
1539 for (unsigned int i = 0; i < face_no; ++i)
1540 result += quadrature[i].size();
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
1577 if (quadrature.size() == 1)
1578 return (face_no +
1579 offset[face_orientation][face_flip][face_rotation] *
1581 quadrature[0].size();
1582 else
1583 {
1584 unsigned int n_points_i = 0;
1585 for (unsigned int i = 0; i < face_no; ++i)
1586 n_points_i += quadrature[i].size();
1587
1588 unsigned int n_points = 0;
1589 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1590 ++i)
1591 n_points += quadrature[i].size();
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 <>
1610 const ReferenceCell reference_cell,
1611 const unsigned int face_no,
1612 const unsigned int subface_no,
1613 const bool,
1614 const bool,
1615 const bool,
1616 const unsigned int n_quadrature_points,
1618{
1619 Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1620 (void)reference_cell;
1621
1625
1626 return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1627 n_quadrature_points);
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,
1640 const unsigned int n_quadrature_points,
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,
1649 n_quadrature_points,
1650 ref_case);
1651}
1652
1653
1654
1655template <>
1658 const ReferenceCell reference_cell,
1659 const unsigned int face_no,
1660 const unsigned int subface_no,
1661 const bool,
1662 const bool,
1663 const bool,
1664 const unsigned int n_quadrature_points,
1666{
1668 (void)reference_cell;
1669
1673
1674 return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1675 n_quadrature_points);
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,
1688 const unsigned int n_quadrature_points,
1689 const internal::SubfaceCase<2> ref_case)
1690{
1691 return subface(ReferenceCells::Quadrilateral,
1692 face_no,
1693 subface_no,
1694 face_orientation,
1695 face_flip,
1696 face_rotation,
1697 n_quadrature_points,
1698 ref_case);
1699}
1700
1701
1702template <>
1705 const ReferenceCell reference_cell,
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,
1711 const unsigned int n_quadrature_points,
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
1818 // refinement instead
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]) *
1955 n_quadrature_points);
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,
1967 const unsigned int n_quadrature_points,
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,
1976 n_quadrature_points,
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>(),
1988 quadrature,
1989 face_no);
1990}
1991
1992
1993
1994template <int dim>
1997 const SubQuadrature &quadrature,
1998 const unsigned int face_no)
1999{
2000 Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
2002 (void)reference_cell;
2003
2004 std::vector<Point<dim>> points(quadrature.size());
2005 project_to_face(quadrature, face_no, points);
2006 return Quadrature<dim>(points, quadrature.get_weights());
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>(),
2019 quadrature,
2020 face_no,
2021 subface_no,
2022 ref_case);
2023}
2024
2025
2026
2027template <int dim>
2030 const SubQuadrature &quadrature,
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
2039 std::vector<Point<dim>> points(quadrature.size());
2040 project_to_subface(quadrature, face_no, subface_no, points, ref_case);
2041 return Quadrature<dim>(points, quadrature.get_weights());
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)
Definition: point.h:111
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
Definition: tensor.h:503
CollectionIterator< T > begin() const
Definition: collection.h:283
unsigned int size() const
Definition: collection.h:264
CollectionIterator< T > end() const
Definition: collection.h:292
unsigned int max_n_quadrature_points() const
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)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition: polynomial.cc:702
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Line
static const unsigned int invalid_unsigned_int
Definition: types.h:201
::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)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)