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