Reference documentation for deal.II version Git 932f7faded 2020-11-28 20:02:43 +0100
\(\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 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 
19 
21 
23 
24 
25 template <int dim>
28 {
29  std::vector<Point<2>> q_points(q.size());
30  std::vector<double> weights(q.size());
31  for (unsigned int i = 0; i < q.size(); ++i)
32  {
33  q_points[i][0] = q.point(i)[1];
34  q_points[i][1] = q.point(i)[0];
35 
36  weights[i] = q.weight(i);
37  }
38 
39  return Quadrature<2>(q_points, weights);
40 }
41 
42 
43 template <int dim>
45 QProjector<dim>::rotate(const Quadrature<2> &q, const unsigned int n_times)
46 {
47  std::vector<Point<2>> q_points(q.size());
48  std::vector<double> weights(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
55  q_points[i][0] = q.point(i)[0];
56  q_points[i][1] = q.point(i)[1];
57  break;
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  weights[i] = q.weight(i);
76  }
77 
78  return Quadrature<2>(q_points, weights);
79 }
80 
81 
82 template <>
83 void
85  const unsigned int face_no,
86  std::vector<Point<1>> &q_points)
87 {
88  project_to_face(ReferenceCell::Type::Line, quadrature, face_no, q_points);
89 }
90 
91 
92 
93 template <>
94 void
96  const Quadrature<0> &,
97  const unsigned int face_no,
98  std::vector<Point<1>> &q_points)
99 {
100  Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
101  (void)reference_cell_type;
102 
103  const unsigned int dim = 1;
105  AssertDimension(q_points.size(), 1);
106 
107  q_points[0] = Point<dim>(static_cast<double>(face_no));
108 }
109 
110 
111 
112 template <>
113 void
115  const unsigned int face_no,
116  std::vector<Point<2>> &q_points)
117 {
118  project_to_face(ReferenceCell::Type::Quad, quadrature, face_no, q_points);
119 }
120 
121 
122 
123 template <>
124 void
126  const Quadrature<1> & quadrature,
127  const unsigned int face_no,
128  std::vector<Point<2>> & q_points)
129 {
130  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
131  (void)reference_cell_type;
132 
133  const unsigned int dim = 2;
135  Assert(q_points.size() == quadrature.size(),
136  ExcDimensionMismatch(q_points.size(), quadrature.size()));
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 
158 
159 
160 template <>
161 void
163  const unsigned int face_no,
164  std::vector<Point<3>> &q_points)
165 {
166  project_to_face(ReferenceCell::Type::Hex, quadrature, face_no, q_points);
167 }
168 
169 
170 
171 template <>
172 void
174  const Quadrature<2> & quadrature,
175  const unsigned int face_no,
176  std::vector<Point<3>> & q_points)
177 {
178  Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
179  (void)reference_cell_type;
180 
181  const unsigned int dim = 3;
183  Assert(q_points.size() == quadrature.size(),
184  ExcDimensionMismatch(q_points.size(), quadrature.size()));
185 
186  for (unsigned int p = 0; p < quadrature.size(); ++p)
187  switch (face_no)
188  {
189  case 0:
190  q_points[p] =
191  Point<dim>(0, quadrature.point(p)(0), quadrature.point(p)(1));
192  break;
193  case 1:
194  q_points[p] =
195  Point<dim>(1, quadrature.point(p)(0), quadrature.point(p)(1));
196  break;
197  case 2:
198  q_points[p] =
199  Point<dim>(quadrature.point(p)(1), 0, quadrature.point(p)(0));
200  break;
201  case 3:
202  q_points[p] =
203  Point<dim>(quadrature.point(p)(1), 1, quadrature.point(p)(0));
204  break;
205  case 4:
206  q_points[p] =
207  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 0);
208  break;
209  case 5:
210  q_points[p] =
211  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 1);
212  break;
213 
214  default:
215  Assert(false, ExcInternalError());
216  }
217 }
218 
219 
220 
221 template <>
222 void
224  const unsigned int face_no,
225  const unsigned int subface_no,
226  std::vector<Point<1>> & q_points,
227  const RefinementCase<0> &ref_case)
228 {
229  project_to_subface(ReferenceCell::Type::Line,
230  quadrature,
231  face_no,
232  subface_no,
233  q_points,
234  ref_case);
235 }
236 
237 
238 
239 template <>
240 void
242  const Quadrature<0> &,
243  const unsigned int face_no,
244  const unsigned int,
245  std::vector<Point<1>> &q_points,
246  const RefinementCase<0> &)
247 {
248  Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
249  (void)reference_cell_type;
250 
251  const unsigned int dim = 1;
253  AssertDimension(q_points.size(), 1);
254 
255  q_points[0] = Point<dim>(static_cast<double>(face_no));
256 }
257 
258 
259 
260 template <>
261 void
263  const unsigned int face_no,
264  const unsigned int subface_no,
265  std::vector<Point<2>> & q_points,
266  const RefinementCase<1> &ref_case)
267 {
268  project_to_subface(ReferenceCell::Type::Quad,
269  quadrature,
270  face_no,
271  subface_no,
272  q_points,
273  ref_case);
274 }
275 
276 
277 
278 template <>
279 void
281  const Quadrature<1> & quadrature,
282  const unsigned int face_no,
283  const unsigned int subface_no,
284  std::vector<Point<2>> & q_points,
285  const RefinementCase<1> &)
286 {
287  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
288  (void)reference_cell_type;
289 
290  const unsigned int dim = 2;
293 
294  Assert(q_points.size() == quadrature.size(),
295  ExcDimensionMismatch(q_points.size(), quadrature.size()));
296 
297  for (unsigned int p = 0; p < quadrature.size(); ++p)
298  switch (face_no)
299  {
300  case 0:
301  switch (subface_no)
302  {
303  case 0:
304  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
305  break;
306  case 1:
307  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
308  break;
309  default:
310  Assert(false, ExcInternalError());
311  }
312  break;
313  case 1:
314  switch (subface_no)
315  {
316  case 0:
317  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
318  break;
319  case 1:
320  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
321  break;
322  default:
323  Assert(false, ExcInternalError());
324  }
325  break;
326  case 2:
327  switch (subface_no)
328  {
329  case 0:
330  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
331  break;
332  case 1:
333  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
334  break;
335  default:
336  Assert(false, ExcInternalError());
337  }
338  break;
339  case 3:
340  switch (subface_no)
341  {
342  case 0:
343  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
344  break;
345  case 1:
346  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
347  break;
348  default:
349  Assert(false, ExcInternalError());
350  }
351  break;
352 
353  default:
354  Assert(false, ExcInternalError());
355  }
356 }
357 
358 
359 
360 template <>
361 void
363  const unsigned int face_no,
364  const unsigned int subface_no,
365  std::vector<Point<3>> & q_points,
366  const RefinementCase<2> &ref_case)
367 {
368  project_to_subface(ReferenceCell::Type::Hex,
369  quadrature,
370  face_no,
371  subface_no,
372  q_points,
373  ref_case);
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 {
387  Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
388  (void)reference_cell_type;
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 {
472  return project_to_all_faces(ReferenceCell::Type::Line, quadrature);
473 }
474 
475 
476 template <>
479  const ReferenceCell::Type reference_cell_type,
480  const Quadrature<0> & quadrature)
481 {
482  Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
483  (void)reference_cell_type;
484 
485  const unsigned int dim = 1;
486 
487  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;
488 
489  // first fix quadrature points
490  std::vector<Point<dim>> q_points;
491  q_points.reserve(n_points * n_faces);
492  std::vector<Point<dim>> help(n_points);
493 
494 
495  // project to each face and append
496  // results
497  for (unsigned int face = 0; face < n_faces; ++face)
498  {
499  project_to_face(quadrature, face, help);
500  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
501  }
502 
503  // next copy over weights
504  std::vector<double> weights;
505  weights.reserve(n_points * n_faces);
506  for (unsigned int face = 0; face < n_faces; ++face)
507  std::copy(quadrature.get_weights().begin(),
508  quadrature.get_weights().end(),
509  std::back_inserter(weights));
510 
511  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
512  Assert(weights.size() == n_points * n_faces, ExcInternalError());
513 
514  return Quadrature<dim>(q_points, weights);
515 }
516 
517 
518 
519 template <>
522  const ReferenceCell::Type reference_cell_type,
523  const SubQuadrature & quadrature)
524 {
525  if (reference_cell_type == ReferenceCell::Type::Tri)
526  {
527  // the quadrature rule to be projected ...
528  const auto &sub_quadrature_points = quadrature.get_points();
529  const auto &sub_quadrature_weights = quadrature.get_weights();
530 
531  // ... on to the faces (defined by its support points and arc length)
532  const std::array<std::pair<std::array<Point<2>, 2>, double>, 3> faces = {
533  {{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0},
534  {{{Point<2>(1.0, 0.0), Point<2>(0.0, 1.0)}}, std::sqrt(2.0)},
535  {{{Point<2>(0.0, 1.0), Point<2>(0.0, 0.0)}}, 1.0}}};
536 
537  // linear polynomial to map the reference quadrature points correctly
538  // on faces
539  const Simplex::ScalarPolynomial<1> poly(1);
540 
541  // new (projected) quadrature points and weights
542  std::vector<Point<2>> points;
543  std::vector<double> weights;
544 
545  // loop over all faces (lines) ...
546  for (const auto &face : faces)
547  // ... and over all possible orientations
548  for (unsigned int orientation = 0; orientation < 2; ++orientation)
549  {
550  std::array<Point<2>, 2> support_points;
551 
552  // determine support point of the current line with the correct
553  // orientation
554  switch (orientation)
555  {
556  case 0:
557  support_points = {{face.first[1], face.first[0]}};
558  break;
559  case 1:
560  support_points = {{face.first[0], face.first[1]}};
561  break;
562  default:
563  Assert(false, ExcNotImplemented());
564  }
565 
566  // loop over all quadrature points
567  for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
568  {
569  Point<2> mapped_point;
570 
571  // map reference quadrature point
572  for (unsigned int i = 0; i < 2; ++i)
573  mapped_point +=
574  support_points[i] *
575  poly.compute_value(i, sub_quadrature_points[j]);
576 
577  points.emplace_back(mapped_point);
578 
579  // scale weight by arc length
580  weights.emplace_back(sub_quadrature_weights[j] * face.second);
581  }
582  }
583 
584  // construct new quadrature rule
585  return {points, weights};
586  }
587 
588  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
589 
590  const unsigned int dim = 2;
591 
592  const unsigned int n_points = quadrature.size(),
594 
595  // first fix quadrature points
596  std::vector<Point<dim>> q_points;
597  q_points.reserve(n_points * n_faces);
598  std::vector<Point<dim>> help(n_points);
599 
600  // project to each face and append
601  // results
602  for (unsigned int face = 0; face < n_faces; ++face)
603  {
604  project_to_face(quadrature, face, help);
605  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
606  }
607 
608  // next copy over weights
609  std::vector<double> weights;
610  weights.reserve(n_points * n_faces);
611  for (unsigned int face = 0; face < n_faces; ++face)
612  std::copy(quadrature.get_weights().begin(),
613  quadrature.get_weights().end(),
614  std::back_inserter(weights));
615 
616  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
617  Assert(weights.size() == n_points * n_faces, ExcInternalError());
618 
619  return Quadrature<dim>(q_points, weights);
620 }
621 
622 
623 
624 template <>
627 {
628  return project_to_all_faces(ReferenceCell::Type::Quad, quadrature);
629 }
630 
631 
632 
633 template <>
636  const ReferenceCell::Type reference_cell_type,
637  const SubQuadrature & quadrature)
638 {
639  if (reference_cell_type == ReferenceCell::Type::Tet)
640  {
641  // the quadrature rule to be projected ...
642  const auto &sub_quadrature_points = quadrature.get_points();
643  const auto &sub_quadrature_weights = quadrature.get_weights();
644 
645  // ... on to the faces (defined by its support points and its area)
646  // note: the area is later not used as a scaling factor but recomputed
647  const std::array<std::pair<std::array<Point<3>, 3>, double>, 4> faces = {
648  {{{{Point<3>(0.0, 0.0, 0.0),
649  Point<3>(1.0, 0.0, 0.0),
650  Point<3>(0.0, 1.0, 0.0)}},
651  0.5},
652  {{{Point<3>(1.0, 0.0, 0.0),
653  Point<3>(0.0, 0.0, 0.0),
654  Point<3>(0.0, 0.0, 1.0)}},
655  0.5},
656  {{{Point<3>(0.0, 0.0, 0.0),
657  Point<3>(0.0, 1.0, 0.0),
658  Point<3>(0.0, 0.0, 1.0)}},
659  0.5},
660  {{{Point<3>(0.0, 1.0, 0.0),
661  Point<3>(1.0, 0.0, 0.0),
662  Point<3>(0.0, 0.0, 1.0)}},
663  0.5 * sqrt(3.0) /*equilateral triangle*/}}};
664 
665  // linear polynomial to map the reference quadrature points correctly
666  // on faces
667  const Simplex::ScalarPolynomial<2> poly(1);
668 
669  // new (projected) quadrature points and weights
670  std::vector<Point<3>> points;
671  std::vector<double> weights;
672 
673  // loop over all faces (triangles) ...
674  for (const auto &face : faces)
675  // ... and over all possible orientations
676  for (unsigned int orientation = 0; orientation < 6; ++orientation)
677  {
678  std::array<Point<3>, 3> support_points;
679 
680  // determine support point of the current line with the correct
681  // orientation
682  switch (orientation)
683  {
684  case 1:
685  support_points = {
686  {face.first[0], face.first[1], face.first[2]}};
687  break;
688  case 3:
689  support_points = {
690  {face.first[1], face.first[0], face.first[2]}};
691  break;
692  case 5:
693  support_points = {
694  {face.first[2], face.first[0], face.first[1]}};
695  break;
696  case 0:
697  support_points = {
698  {face.first[0], face.first[2], face.first[1]}};
699  break;
700  case 2:
701  support_points = {
702  {face.first[1], face.first[2], face.first[0]}};
703  break;
704  case 4:
705  support_points = {
706  {face.first[2], face.first[1], face.first[0]}};
707  break;
708  default:
709  Assert(false, ExcNotImplemented());
710  }
711 
712  // loop over all quadrature points
713  for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
714  {
715  Point<3> mapped_point;
716 
717  // map reference quadrature point
718  for (unsigned int i = 0; i < 3; ++i)
719  mapped_point +=
720  support_points[i] *
721  poly.compute_value(i, sub_quadrature_points[j]);
722 
723  points.push_back(mapped_point);
724 
725  // scale quadrature weight
726  const double scaling = [&]() {
727  const auto &supp_pts = support_points;
728 
729  const unsigned int n_shape_functions = 3;
730  const unsigned int dim_ = 2;
731  const unsigned int spacedim = 3;
732 
733  double result[spacedim][dim_];
734 
735  std::vector<Tensor<1, dim_>> shape_derivatives(
736  n_shape_functions);
737 
738  for (unsigned int i = 0; i < 3; ++i)
739  shape_derivatives[i] =
740  poly.compute_1st_derivative(i, sub_quadrature_points[j]);
741 
742  for (unsigned int i = 0; i < spacedim; ++i)
743  for (unsigned int j = 0; j < dim_; ++j)
744  result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
745  for (unsigned int k = 1; k < n_shape_functions; ++k)
746  for (unsigned int i = 0; i < spacedim; ++i)
747  for (unsigned int j = 0; j < dim_; ++j)
748  result[i][j] +=
749  shape_derivatives[k][j] * supp_pts[k][i];
750 
751  DerivativeForm<1, dim_, spacedim> contravariant;
752 
753  for (unsigned int i = 0; i < spacedim; ++i)
754  for (unsigned int j = 0; j < dim_; ++j)
755  contravariant[i][j] = result[i][j];
756 
757 
758  Tensor<1, spacedim> DX_t[dim_];
759  for (unsigned int i = 0; i < spacedim; ++i)
760  for (unsigned int j = 0; j < dim_; ++j)
761  DX_t[j][i] = contravariant[i][j];
762 
763  Tensor<2, dim_> G;
764  for (unsigned int i = 0; i < dim_; ++i)
765  for (unsigned int j = 0; j < dim_; ++j)
766  G[i][j] = DX_t[i] * DX_t[j];
767 
768  return std::sqrt(determinant(G));
769  }();
770 
771  weights.push_back(sub_quadrature_weights[j] * scaling);
772  }
773  }
774 
775  // construct new quadrature rule
776  return {points, weights};
777  }
778 
779 
780  Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
781 
782  const unsigned int dim = 3;
783 
784  SubQuadrature q_reflected = reflect(quadrature);
785  SubQuadrature q[8] = {quadrature,
786  rotate(quadrature, 1),
787  rotate(quadrature, 2),
788  rotate(quadrature, 3),
789  q_reflected,
790  rotate(q_reflected, 3),
791  rotate(q_reflected, 2),
792  rotate(q_reflected, 1)};
793 
794 
795 
796  const unsigned int n_points = quadrature.size(),
798 
799  // first fix quadrature points
800  std::vector<Point<dim>> q_points;
801  q_points.reserve(n_points * n_faces * 8);
802  std::vector<Point<dim>> help(n_points);
803 
804  std::vector<double> weights;
805  weights.reserve(n_points * n_faces * 8);
806 
807  // do the following for all possible
808  // mutations of a face (mutation==0
809  // corresponds to a face with standard
810  // orientation, no flip and no rotation)
811  for (const auto &mutation : q)
812  {
813  // project to each face and append
814  // results
815  for (unsigned int face = 0; face < n_faces; ++face)
816  {
817  project_to_face(mutation, face, help);
818  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
819  }
820 
821  // next copy over weights
822  for (unsigned int face = 0; face < n_faces; ++face)
823  std::copy(mutation.get_weights().begin(),
824  mutation.get_weights().end(),
825  std::back_inserter(weights));
826  }
827 
828 
829  Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError());
830  Assert(weights.size() == n_points * n_faces * 8, ExcInternalError());
831 
832  return Quadrature<dim>(q_points, weights);
833 }
834 
835 
836 
837 template <>
840 {
841  return project_to_all_faces(ReferenceCell::Type::Hex, quadrature);
842 }
843 
844 
845 
846 template <>
849 {
850  return project_to_all_subfaces(ReferenceCell::Type::Line, quadrature);
851 }
852 
853 
854 
855 template <>
858  const ReferenceCell::Type reference_cell_type,
859  const Quadrature<0> & quadrature)
860 {
861  Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
862  (void)reference_cell_type;
863 
864  const unsigned int dim = 1;
865 
866  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
867  subfaces_per_face =
869 
870  // first fix quadrature points
871  std::vector<Point<dim>> q_points;
872  q_points.reserve(n_points * n_faces * subfaces_per_face);
873  std::vector<Point<dim>> help(n_points);
874 
875  // project to each face and copy
876  // results
877  for (unsigned int face = 0; face < n_faces; ++face)
878  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
879  {
880  project_to_subface(quadrature, face, subface, help);
881  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
882  }
883 
884  // next copy over weights
885  std::vector<double> weights;
886  weights.reserve(n_points * n_faces * subfaces_per_face);
887  for (unsigned int face = 0; face < n_faces; ++face)
888  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
889  std::copy(quadrature.get_weights().begin(),
890  quadrature.get_weights().end(),
891  std::back_inserter(weights));
892 
893  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
894  ExcInternalError());
895  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
896  ExcInternalError());
897 
898  return Quadrature<dim>(q_points, weights);
899 }
900 
901 
902 
903 template <>
906  const ReferenceCell::Type reference_cell_type,
907  const SubQuadrature & quadrature)
908 {
909  if (reference_cell_type == ReferenceCell::Type::Tri ||
910  reference_cell_type == ReferenceCell::Type::Tet)
911  return Quadrature<2>(); // nothing to do
912 
913  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
914 
915  const unsigned int dim = 2;
916 
917  const unsigned int n_points = quadrature.size(),
919  subfaces_per_face =
921 
922  // first fix quadrature points
923  std::vector<Point<dim>> q_points;
924  q_points.reserve(n_points * n_faces * subfaces_per_face);
925  std::vector<Point<dim>> help(n_points);
926 
927  // project to each face and copy
928  // results
929  for (unsigned int face = 0; face < n_faces; ++face)
930  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
931  {
932  project_to_subface(quadrature, face, subface, help);
933  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
934  }
935 
936  // next copy over weights
937  std::vector<double> weights;
938  weights.reserve(n_points * n_faces * subfaces_per_face);
939  for (unsigned int face = 0; face < n_faces; ++face)
940  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
941  std::copy(quadrature.get_weights().begin(),
942  quadrature.get_weights().end(),
943  std::back_inserter(weights));
944 
945  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
946  ExcInternalError());
947  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
948  ExcInternalError());
949 
950  return Quadrature<dim>(q_points, weights);
951 }
952 
953 
954 
955 template <>
958 {
959  return project_to_all_subfaces(ReferenceCell::Type::Quad, quadrature);
960 }
961 
962 
963 
964 template <>
967  const ReferenceCell::Type reference_cell_type,
968  const SubQuadrature & quadrature)
969 {
970  if (reference_cell_type == ReferenceCell::Type::Tri ||
971  reference_cell_type == ReferenceCell::Type::Tet)
972  return Quadrature<3>(); // nothing to do
973 
974  Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
975 
976  const unsigned int dim = 3;
977  SubQuadrature q_reflected = reflect(quadrature);
978  SubQuadrature q[8] = {quadrature,
979  rotate(quadrature, 1),
980  rotate(quadrature, 2),
981  rotate(quadrature, 3),
982  q_reflected,
983  rotate(q_reflected, 3),
984  rotate(q_reflected, 2),
985  rotate(q_reflected, 1)};
986 
987  const unsigned int n_points = quadrature.size(),
989  total_subfaces_per_face = 2 + 2 + 4;
990 
991  // first fix quadrature points
992  std::vector<Point<dim>> q_points;
993  q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
994  std::vector<Point<dim>> help(n_points);
995 
996  std::vector<double> weights;
997  weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
998 
999  // do the following for all possible
1000  // mutations of a face (mutation==0
1001  // corresponds to a face with standard
1002  // orientation, no flip and no rotation)
1003  for (const auto &mutation : q)
1004  {
1005  // project to each face and copy
1006  // results
1007  for (unsigned int face = 0; face < n_faces; ++face)
1008  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1009  ref_case >= RefinementCase<dim - 1>::cut_x;
1010  --ref_case)
1011  for (unsigned int subface = 0;
1013  RefinementCase<dim - 1>(ref_case));
1014  ++subface)
1015  {
1016  project_to_subface(mutation,
1017  face,
1018  subface,
1019  help,
1020  RefinementCase<dim - 1>(ref_case));
1021  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1022  }
1023 
1024  // next copy over weights
1025  for (unsigned int face = 0; face < n_faces; ++face)
1026  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1027  ref_case >= RefinementCase<dim - 1>::cut_x;
1028  --ref_case)
1029  for (unsigned int subface = 0;
1031  RefinementCase<dim - 1>(ref_case));
1032  ++subface)
1033  std::copy(mutation.get_weights().begin(),
1034  mutation.get_weights().end(),
1035  std::back_inserter(weights));
1036  }
1037 
1038  Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1039  ExcInternalError());
1040  Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1041  ExcInternalError());
1042 
1043  return Quadrature<dim>(q_points, weights);
1044 }
1045 
1046 
1047 
1048 template <>
1051 {
1052  return project_to_all_subfaces(ReferenceCell::Type::Hex, quadrature);
1053 }
1054 
1055 
1056 
1057 // This function is not used in the library
1058 template <int dim>
1061  const unsigned int child_no)
1062 {
1063  return project_to_child(ReferenceCell::get_hypercube(dim),
1064  quadrature,
1065  child_no);
1066 }
1067 
1068 
1069 
1070 template <int dim>
1073  const Quadrature<dim> & quadrature,
1074  const unsigned int child_no)
1075 {
1076  Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
1077  ExcNotImplemented());
1078  (void)reference_cell_type;
1079 
1081 
1082  const unsigned int n_q_points = quadrature.size();
1083 
1084  std::vector<Point<dim>> q_points(n_q_points);
1085  for (unsigned int i = 0; i < n_q_points; ++i)
1086  q_points[i] =
1088  child_no);
1089 
1090  // for the weights, things are
1091  // equally simple: copy them and
1092  // scale them
1093  std::vector<double> weights = quadrature.get_weights();
1094  for (unsigned int i = 0; i < n_q_points; ++i)
1095  weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1096 
1097  return Quadrature<dim>(q_points, weights);
1098 }
1099 
1100 
1101 
1102 template <int dim>
1105 {
1106  return project_to_all_children(ReferenceCell::get_hypercube(dim), quadrature);
1107 }
1108 
1109 
1110 
1111 template <int dim>
1114  const ReferenceCell::Type reference_cell_type,
1115  const Quadrature<dim> & quadrature)
1116 {
1117  Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
1118  ExcNotImplemented());
1119  (void)reference_cell_type;
1120 
1121  const unsigned int n_points = quadrature.size(),
1123 
1124  std::vector<Point<dim>> q_points(n_points * n_children);
1125  std::vector<double> weights(n_points * n_children);
1126 
1127  // project to each child and copy
1128  // results
1129  for (unsigned int child = 0; child < n_children; ++child)
1130  {
1131  Quadrature<dim> help = project_to_child(quadrature, child);
1132  for (unsigned int i = 0; i < n_points; ++i)
1133  {
1134  q_points[child * n_points + i] = help.point(i);
1135  weights[child * n_points + i] = help.weight(i);
1136  }
1137  }
1138  return Quadrature<dim>(q_points, weights);
1139 }
1140 
1141 
1142 
1143 template <int dim>
1146  const Point<dim> & p1,
1147  const Point<dim> & p2)
1148 {
1149  return project_to_line(ReferenceCell::get_hypercube(dim), quadrature, p1, p2);
1150 }
1151 
1152 
1153 
1154 template <int dim>
1157  const Quadrature<1> & quadrature,
1158  const Point<dim> & p1,
1159  const Point<dim> & p2)
1160 {
1161  Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
1162  ExcNotImplemented());
1163  (void)reference_cell_type;
1164 
1165  const unsigned int n = quadrature.size();
1166  std::vector<Point<dim>> points(n);
1167  std::vector<double> weights(n);
1168  const double length = p1.distance(p2);
1169 
1170  for (unsigned int k = 0; k < n; ++k)
1171  {
1172  const double alpha = quadrature.point(k)(0);
1173  points[k] = alpha * p2;
1174  points[k] += (1. - alpha) * p1;
1175  weights[k] = length * quadrature.weight(k);
1176  }
1177  return Quadrature<dim>(points, weights);
1178 }
1179 
1180 
1181 
1182 template <int dim>
1184 QProjector<dim>::DataSetDescriptor::face(const unsigned int face_no,
1185  const bool face_orientation,
1186  const bool face_flip,
1187  const bool face_rotation,
1188  const unsigned int n_quadrature_points)
1189 {
1190  return face(ReferenceCell::get_hypercube(dim),
1191  face_no,
1192  face_orientation,
1193  face_flip,
1194  face_rotation,
1195  n_quadrature_points);
1196 }
1197 
1198 
1199 
1200 template <int dim>
1203  const ReferenceCell::Type reference_cell_type,
1204  const unsigned int face_no,
1205  const bool face_orientation,
1206  const bool face_flip,
1207  const bool face_rotation,
1208  const unsigned int n_quadrature_points)
1209 {
1210  if (reference_cell_type == ReferenceCell::Type::Tri ||
1211  reference_cell_type == ReferenceCell::Type::Tet)
1212  {
1213  if (dim == 2)
1214  return {(2 * face_no + face_orientation) * n_quadrature_points};
1215  else if (dim == 3)
1216  {
1217  const unsigned int orientation =
1218  (face_flip * 2 + face_rotation) * 2 + face_orientation;
1219  return {(6 * face_no + orientation) * n_quadrature_points};
1220  }
1221  }
1222 
1223  Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
1224  ExcNotImplemented());
1225 
1227 
1228  switch (dim)
1229  {
1230  case 1:
1231  case 2:
1232  return face_no * n_quadrature_points;
1233 
1234 
1235  case 3:
1236  {
1237  // in 3d, we have to account for faces that
1238  // have non-standard face orientation, flip
1239  // and rotation. thus, we have to store
1240  // _eight_ data sets per face or subface
1241 
1242  // set up a table with the according offsets
1243  // for non-standard orientation, first index:
1244  // face_orientation (standard true=1), second
1245  // index: face_flip (standard false=0), third
1246  // index: face_rotation (standard false=0)
1247  //
1248  // note, that normally we should use the
1249  // obvious offsets 0,1,2,3,4,5,6,7. However,
1250  // prior to the changes enabling flipped and
1251  // rotated faces, in many places of the
1252  // library the convention was used, that the
1253  // first dataset with offset 0 corresponds to
1254  // a face in standard orientation. therefore
1255  // we use the offsets 4,5,6,7,0,1,2,3 here to
1256  // stick to that (implicit) convention
1257  static const unsigned int offset[2][2][2] = {
1259  5 * GeometryInfo<dim>::
1260  faces_per_cell}, // face_orientation=false; face_flip=false;
1261  // face_rotation=false and true
1263  7 * GeometryInfo<dim>::
1264  faces_per_cell}}, // face_orientation=false; face_flip=true;
1265  // face_rotation=false and true
1267  1 * GeometryInfo<dim>::
1268  faces_per_cell}, // face_orientation=true; face_flip=false;
1269  // face_rotation=false and true
1271  3 * GeometryInfo<dim>::
1272  faces_per_cell}}}; // face_orientation=true; face_flip=true;
1273  // face_rotation=false and true
1274 
1275  return (
1276  (face_no + offset[face_orientation][face_flip][face_rotation]) *
1277  n_quadrature_points);
1278  }
1279 
1280  default:
1281  Assert(false, ExcInternalError());
1282  }
1284 }
1285 
1286 
1287 
1288 template <>
1291  const ReferenceCell::Type reference_cell_type,
1292  const unsigned int face_no,
1293  const unsigned int subface_no,
1294  const bool,
1295  const bool,
1296  const bool,
1297  const unsigned int n_quadrature_points,
1299 {
1300  Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
1301  (void)reference_cell_type;
1302 
1305  ExcInternalError());
1306 
1307  return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1308  n_quadrature_points);
1309 }
1310 
1311 
1312 
1313 template <>
1316  const unsigned int face_no,
1317  const unsigned int subface_no,
1318  const bool face_orientation,
1319  const bool face_flip,
1320  const bool face_rotation,
1321  const unsigned int n_quadrature_points,
1322  const internal::SubfaceCase<1> ref_case)
1323 {
1325  face_no,
1326  subface_no,
1327  face_orientation,
1328  face_flip,
1329  face_rotation,
1330  n_quadrature_points,
1331  ref_case);
1332 }
1333 
1334 
1335 
1336 template <>
1339  const ReferenceCell::Type reference_cell_type,
1340  const unsigned int face_no,
1341  const unsigned int subface_no,
1342  const bool,
1343  const bool,
1344  const bool,
1345  const unsigned int n_quadrature_points,
1347 {
1348  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
1349  (void)reference_cell_type;
1350 
1353  ExcInternalError());
1354 
1355  return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1356  n_quadrature_points);
1357 }
1358 
1359 
1360 
1361 template <>
1364  const unsigned int face_no,
1365  const unsigned int subface_no,
1366  const bool face_orientation,
1367  const bool face_flip,
1368  const bool face_rotation,
1369  const unsigned int n_quadrature_points,
1370  const internal::SubfaceCase<2> ref_case)
1371 {
1373  face_no,
1374  subface_no,
1375  face_orientation,
1376  face_flip,
1377  face_rotation,
1378  n_quadrature_points,
1379  ref_case);
1380 }
1381 
1382 
1383 template <>
1386  const ReferenceCell::Type reference_cell_type,
1387  const unsigned int face_no,
1388  const unsigned int subface_no,
1389  const bool face_orientation,
1390  const bool face_flip,
1391  const bool face_rotation,
1392  const unsigned int n_quadrature_points,
1393  const internal::SubfaceCase<3> ref_case)
1394 {
1395  const unsigned int dim = 3;
1396 
1397  Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
1398  (void)reference_cell_type;
1399 
1402  ExcInternalError());
1403 
1404  // As the quadrature points created by
1405  // QProjector are on subfaces in their
1406  // "standard location" we have to use a
1407  // permutation of the equivalent subface
1408  // number in order to respect face
1409  // orientation, flip and rotation. The
1410  // information we need here is exactly the
1411  // same as the
1412  // GeometryInfo<3>::child_cell_on_face info
1413  // for the bottom face (face 4) of a hex, as
1414  // on this the RefineCase of the cell matches
1415  // that of the face and the subfaces are
1416  // numbered in the same way as the child
1417  // cells.
1418 
1419  // in 3d, we have to account for faces that
1420  // have non-standard face orientation, flip
1421  // and rotation. thus, we have to store
1422  // _eight_ data sets per face or subface
1423  // already for the isotropic
1424  // case. Additionally, we have three
1425  // different refinement cases, resulting in
1426  // <tt>4 + 2 + 2 = 8</tt> different subfaces
1427  // for each face.
1428  const unsigned int total_subfaces_per_face = 8;
1429 
1430  // set up a table with the according offsets
1431  // for non-standard orientation, first index:
1432  // face_orientation (standard true=1), second
1433  // index: face_flip (standard false=0), third
1434  // index: face_rotation (standard false=0)
1435  //
1436  // note, that normally we should use the
1437  // obvious offsets 0,1,2,3,4,5,6,7. However,
1438  // prior to the changes enabling flipped and
1439  // rotated faces, in many places of the
1440  // library the convention was used, that the
1441  // first dataset with offset 0 corresponds to
1442  // a face in standard orientation. therefore
1443  // we use the offsets 4,5,6,7,0,1,2,3 here to
1444  // stick to that (implicit) convention
1445  static const unsigned int orientation_offset[2][2][2] = {
1446  {// face_orientation=false; face_flip=false; face_rotation=false and true
1447  {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1448  5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1449  // face_orientation=false; face_flip=true; face_rotation=false and true
1450  {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1451  7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
1452  {// face_orientation=true; face_flip=false; face_rotation=false and true
1453  {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1454  1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1455  // face_orientation=true; face_flip=true; face_rotation=false and true
1456  {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1457  3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};
1458 
1459  // set up a table with the offsets for a
1460  // given refinement case respecting the
1461  // corresponding number of subfaces. the
1462  // index corresponds to (RefineCase::Type - 1)
1463 
1464  // note, that normally we should use the
1465  // obvious offsets 0,2,6. However, prior to
1466  // the implementation of anisotropic
1467  // refinement, in many places of the library
1468  // the convention was used, that the first
1469  // dataset with offset 0 corresponds to a
1470  // standard (isotropic) face
1471  // refinement. therefore we use the offsets
1472  // 6,4,0 here to stick to that (implicit)
1473  // convention
1474  static const unsigned int ref_case_offset[3] = {
1475  6, // cut_x
1476  4, // cut_y
1477  0 // cut_xy
1478  };
1479 
1480 
1481  // for each subface of a given FaceRefineCase
1482  // there is a corresponding equivalent
1483  // subface number of one of the "standard"
1484  // RefineCases (cut_x, cut_y, cut_xy). Map
1485  // the given values to those equivalent
1486  // ones.
1487 
1488  // first, define an invalid number
1489  static const unsigned int e = numbers::invalid_unsigned_int;
1490 
1491  static const RefinementCase<dim - 1>
1492  equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
1494  // case_none. there should be only
1495  // invalid values here. However, as
1496  // this function is also called (in
1497  // tests) for cells which have no
1498  // refined faces, use isotropic
1499  // refinement instead
1500  {RefinementCase<dim - 1>::cut_xy,
1501  RefinementCase<dim - 1>::cut_xy,
1502  RefinementCase<dim - 1>::cut_xy,
1503  RefinementCase<dim - 1>::cut_xy},
1504  // case_x
1505  {RefinementCase<dim - 1>::cut_x,
1506  RefinementCase<dim - 1>::cut_x,
1507  RefinementCase<dim - 1>::no_refinement,
1508  RefinementCase<dim - 1>::no_refinement},
1509  // case_x1y
1510  {RefinementCase<dim - 1>::cut_xy,
1511  RefinementCase<dim - 1>::cut_xy,
1512  RefinementCase<dim - 1>::cut_x,
1513  RefinementCase<dim - 1>::no_refinement},
1514  // case_x2y
1515  {RefinementCase<dim - 1>::cut_x,
1516  RefinementCase<dim - 1>::cut_xy,
1517  RefinementCase<dim - 1>::cut_xy,
1518  RefinementCase<dim - 1>::no_refinement},
1519  // case_x1y2y
1520  {RefinementCase<dim - 1>::cut_xy,
1521  RefinementCase<dim - 1>::cut_xy,
1522  RefinementCase<dim - 1>::cut_xy,
1523  RefinementCase<dim - 1>::cut_xy},
1524  // case_y
1525  {RefinementCase<dim - 1>::cut_y,
1526  RefinementCase<dim - 1>::cut_y,
1527  RefinementCase<dim - 1>::no_refinement,
1528  RefinementCase<dim - 1>::no_refinement},
1529  // case_y1x
1530  {RefinementCase<dim - 1>::cut_xy,
1531  RefinementCase<dim - 1>::cut_xy,
1532  RefinementCase<dim - 1>::cut_y,
1533  RefinementCase<dim - 1>::no_refinement},
1534  // case_y2x
1535  {RefinementCase<dim - 1>::cut_y,
1536  RefinementCase<dim - 1>::cut_xy,
1537  RefinementCase<dim - 1>::cut_xy,
1538  RefinementCase<dim - 1>::no_refinement},
1539  // case_y1x2x
1540  {RefinementCase<dim - 1>::cut_xy,
1541  RefinementCase<dim - 1>::cut_xy,
1542  RefinementCase<dim - 1>::cut_xy,
1543  RefinementCase<dim - 1>::cut_xy},
1544  // case_xy (case_isotropic)
1545  {RefinementCase<dim - 1>::cut_xy,
1546  RefinementCase<dim - 1>::cut_xy,
1547  RefinementCase<dim - 1>::cut_xy,
1548  RefinementCase<dim - 1>::cut_xy}};
1549 
1550  static const unsigned int
1551  equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
1553  // case_none, see above
1554  {0, 1, 2, 3},
1555  // case_x
1556  {0, 1, e, e},
1557  // case_x1y
1558  {0, 2, 1, e},
1559  // case_x2y
1560  {0, 1, 3, e},
1561  // case_x1y2y
1562  {0, 2, 1, 3},
1563  // case_y
1564  {0, 1, e, e},
1565  // case_y1x
1566  {0, 1, 1, e},
1567  // case_y2x
1568  {0, 2, 3, e},
1569  // case_y1x2x
1570  {0, 1, 2, 3},
1571  // case_xy (case_isotropic)
1572  {0, 1, 2, 3}};
1573 
1574  // If face-orientation or face_rotation are
1575  // non-standard, cut_x and cut_y have to be
1576  // exchanged.
1577  static const RefinementCase<dim - 1> ref_case_permutation[4] = {
1578  RefinementCase<dim - 1>::no_refinement,
1579  RefinementCase<dim - 1>::cut_y,
1580  RefinementCase<dim - 1>::cut_x,
1581  RefinementCase<dim - 1>::cut_xy};
1582 
1583  // set a corresponding (equivalent)
1584  // RefineCase and subface number
1585  const RefinementCase<dim - 1> equ_ref_case =
1586  equivalent_refine_case[ref_case][subface_no];
1587  const unsigned int equ_subface_no =
1588  equivalent_subface_number[ref_case][subface_no];
1589  // make sure, that we got a valid subface and RefineCase
1591  ExcInternalError());
1592  Assert(equ_subface_no != e, ExcInternalError());
1593  // now, finally respect non-standard faces
1594  const RefinementCase<dim - 1> final_ref_case =
1595  (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1596  equ_ref_case);
1597 
1598  // what we have now is the number of
1599  // the subface in the natural
1600  // orientation of the *face*. what we
1601  // need to know is the number of the
1602  // subface concerning the standard face
1603  // orientation as seen from the *cell*.
1604 
1605  // this mapping is not trivial, but we
1606  // have done exactly this stuff in the
1607  // child_cell_on_face function. in
1608  // order to reduce the amount of code
1609  // as well as to make maintaining the
1610  // functionality easier we want to
1611  // reuse that information. So we note
1612  // that on the bottom face (face 4) of
1613  // a hex cell the local x and y
1614  // coordinates of the face and the cell
1615  // coincide, thus also the refinement
1616  // case of the face corresponds to the
1617  // refinement case of the cell
1618  // (ignoring cell refinement along the
1619  // z direction). Using this knowledge
1620  // we can (ab)use the
1621  // child_cell_on_face function to do
1622  // exactly the transformation we are in
1623  // need of now
1624  const unsigned int final_subface_no =
1626  4,
1627  equ_subface_no,
1628  face_orientation,
1629  face_flip,
1630  face_rotation,
1631  equ_ref_case);
1632 
1633  return (((face_no * total_subfaces_per_face +
1634  ref_case_offset[final_ref_case - 1] + final_subface_no) +
1635  orientation_offset[face_orientation][face_flip][face_rotation]) *
1636  n_quadrature_points);
1637 }
1638 
1639 
1640 template <>
1643  const unsigned int face_no,
1644  const unsigned int subface_no,
1645  const bool face_orientation,
1646  const bool face_flip,
1647  const bool face_rotation,
1648  const unsigned int n_quadrature_points,
1649  const internal::SubfaceCase<3> ref_case)
1650 {
1652  face_no,
1653  subface_no,
1654  face_orientation,
1655  face_flip,
1656  face_rotation,
1657  n_quadrature_points,
1658  ref_case);
1659 }
1660 
1661 
1662 
1663 template <int dim>
1666  const unsigned int face_no)
1667 {
1669  quadrature,
1670  face_no);
1671 }
1672 
1673 
1674 
1675 template <int dim>
1678  const SubQuadrature & quadrature,
1679  const unsigned int face_no)
1680 {
1681  Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
1682  ExcNotImplemented());
1683  (void)reference_cell_type;
1684 
1685  std::vector<Point<dim>> points(quadrature.size());
1686  project_to_face(quadrature, face_no, points);
1687  return Quadrature<dim>(points, quadrature.get_weights());
1688 }
1689 
1690 
1691 
1692 template <int dim>
1695  const unsigned int face_no,
1696  const unsigned int subface_no,
1697  const RefinementCase<dim - 1> &ref_case)
1698 {
1700  quadrature,
1701  face_no,
1702  subface_no,
1703  ref_case);
1704 }
1705 
1706 
1707 
1708 template <int dim>
1711  const ReferenceCell::Type reference_cell_type,
1712  const SubQuadrature & quadrature,
1713  const unsigned int face_no,
1714  const unsigned int subface_no,
1715  const RefinementCase<dim - 1> &ref_case)
1716 {
1717  Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
1718  ExcNotImplemented());
1719  (void)reference_cell_type;
1720 
1721  std::vector<Point<dim>> points(quadrature.size());
1722  project_to_subface(quadrature, face_no, subface_no, points, ref_case);
1723  return Quadrature<dim>(points, quadrature.get_weights());
1724 }
1725 
1726 
1727 // explicit instantiations; note: we need them all for all dimensions
1728 template class QProjector<1>;
1729 template class QProjector<2>;
1730 template class QProjector<3>;
1731 
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: qprojector.cc:1060
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
const std::vector< Point< dim > > & get_points() const
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim >> &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
const std::vector< double > & get_weights() const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Type get_hypercube(const unsigned int dim)
static Quadrature< 2 > rotate(const Quadrature< 2 > &q, const unsigned int n_times)
Definition: qprojector.cc:45
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
const Point< dim > & point(const unsigned int i) const
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
Definition: polynomials.cc:419
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
static Quadrature< dim > project_to_line(const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
Definition: qprojector.cc:1145
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
void rotate(const double angle, Triangulation< dim > &triangulation)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
unsigned int size() const
static Quadrature< 2 > reflect(const Quadrature< 2 > &q)
Definition: qprojector.cc:27
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
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)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static ::ExceptionBase & ExcNotImplemented()
double compute_value(const unsigned int i, const Point< dim > &p) const override
Definition: polynomials.cc:67
void copy(const T *begin, const T *end, U *dest)
static Quadrature< dim > project_to_all_children(const Quadrature< dim > &quadrature)
Definition: qprojector.cc:1104
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: qprojector.cc:1184