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