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