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