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