Reference documentation for deal.II version 9.3.2
\(\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 - 2021 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) * n_quadrature_points};
1395  else if (dim == 3)
1396  {
1397  const unsigned int orientation =
1398  (face_flip * 2 + face_rotation) * 2 + face_orientation;
1399  return {(6 * face_no + orientation) * n_quadrature_points};
1400  }
1401  }
1402 
1403  Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1404  ExcNotImplemented());
1405 
1406  Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
1407 
1408  switch (dim)
1409  {
1410  case 1:
1411  case 2:
1412  return face_no * n_quadrature_points;
1413 
1414 
1415  case 3:
1416  {
1417  // in 3d, we have to account for faces that
1418  // have non-standard face orientation, flip
1419  // and rotation. thus, we have to store
1420  // _eight_ data sets per face or subface
1421 
1422  // set up a table with the according offsets
1423  // for non-standard orientation, first index:
1424  // face_orientation (standard true=1), second
1425  // index: face_flip (standard false=0), third
1426  // index: face_rotation (standard false=0)
1427  //
1428  // note, that normally we should use the
1429  // obvious offsets 0,1,2,3,4,5,6,7. However,
1430  // prior to the changes enabling flipped and
1431  // rotated faces, in many places of the
1432  // library the convention was used, that the
1433  // first dataset with offset 0 corresponds to
1434  // a face in standard orientation. therefore
1435  // we use the offsets 4,5,6,7,0,1,2,3 here to
1436  // stick to that (implicit) convention
1437  static const unsigned int offset[2][2][2] = {
1438  {{4 * GeometryInfo<dim>::faces_per_cell,
1439  5 * GeometryInfo<dim>::
1440  faces_per_cell}, // face_orientation=false; face_flip=false;
1441  // face_rotation=false and true
1442  {6 * GeometryInfo<dim>::faces_per_cell,
1443  7 * GeometryInfo<dim>::
1444  faces_per_cell}}, // face_orientation=false; face_flip=true;
1445  // face_rotation=false and true
1446  {{0 * GeometryInfo<dim>::faces_per_cell,
1447  1 * GeometryInfo<dim>::
1448  faces_per_cell}, // face_orientation=true; face_flip=false;
1449  // face_rotation=false and true
1450  {2 * GeometryInfo<dim>::faces_per_cell,
1451  3 * GeometryInfo<dim>::
1452  faces_per_cell}}}; // face_orientation=true; face_flip=true;
1453  // face_rotation=false and true
1454 
1455  return (
1456  (face_no + offset[face_orientation][face_flip][face_rotation]) *
1457  n_quadrature_points);
1458  }
1459 
1460  default:
1461  Assert(false, ExcInternalError());
1462  }
1463  return numbers::invalid_unsigned_int;
1464 }
1465 
1466 
1467 
1468 template <int dim>
1469 typename QProjector<dim>::DataSetDescriptor
1470 QProjector<dim>::DataSetDescriptor::face(
1471  const ReferenceCell reference_cell,
1472  const unsigned int face_no,
1473  const bool face_orientation,
1474  const bool face_flip,
1475  const bool face_rotation,
1476  const hp::QCollection<dim - 1> &quadrature)
1477 {
1478  if (reference_cell == ReferenceCells::Triangle ||
1479  reference_cell == ReferenceCells::Tetrahedron ||
1480  reference_cell == ReferenceCells::Wedge ||
1481  reference_cell == ReferenceCells::Pyramid)
1482  {
1483  unsigned int offset = 0;
1484 
1485  static const unsigned int X = numbers::invalid_unsigned_int;
1486  static const std::array<unsigned int, 5> scale_tri = {{2, 2, 2, X, X}};
1487  static const std::array<unsigned int, 5> scale_tet = {{6, 6, 6, 6, X}};
1488  static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
1489  static const std::array<unsigned int, 5> scale_pyramid = {
1490  {8, 6, 6, 6, 6}};
1491 
1492  const auto &scale =
1493  (reference_cell == ReferenceCells::Triangle) ?
1494  scale_tri :
1495  ((reference_cell == ReferenceCells::Tetrahedron) ?
1496  scale_tet :
1497  ((reference_cell == ReferenceCells::Wedge) ? scale_wedge :
1498  scale_pyramid));
1499 
1500  if (quadrature.size() == 1)
1501  offset = scale[0] * quadrature[0].size() * face_no;
1502  else
1503  for (unsigned int i = 0; i < face_no; ++i)
1504  offset += scale[i] * quadrature[i].size();
1505 
1506  if (dim == 2)
1507  return {offset +
1508  face_orientation *
1509  quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1510  else if (dim == 3)
1511  {
1512  const unsigned int orientation =
1513  (face_flip * 2 + face_rotation) * 2 + face_orientation;
1514 
1515  return {offset +
1516  orientation *
1517  quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
1518  }
1519  }
1520 
1521  Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1522  ExcNotImplemented());
1523 
1524  Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
1525 
1526  switch (dim)
1527  {
1528  case 1:
1529  case 2:
1530  {
1531  if (quadrature.size() == 1)
1532  return quadrature[0].size() * face_no;
1533  else
1534  {
1535  unsigned int result = 0;
1536  for (unsigned int i = 0; i < face_no; ++i)
1537  result += quadrature[i].size();
1538  return result;
1539  }
1540  }
1541  case 3:
1542  {
1543  // in 3d, we have to account for faces that
1544  // have non-standard face orientation, flip
1545  // and rotation. thus, we have to store
1546  // _eight_ data sets per face or subface
1547 
1548  // set up a table with the according offsets
1549  // for non-standard orientation, first index:
1550  // face_orientation (standard true=1), second
1551  // index: face_flip (standard false=0), third
1552  // index: face_rotation (standard false=0)
1553  //
1554  // note, that normally we should use the
1555  // obvious offsets 0,1,2,3,4,5,6,7. However,
1556  // prior to the changes enabling flipped and
1557  // rotated faces, in many places of the
1558  // library the convention was used, that the
1559  // first dataset with offset 0 corresponds to
1560  // a face in standard orientation. therefore
1561  // we use the offsets 4,5,6,7,0,1,2,3 here to
1562  // stick to that (implicit) convention
1563  static const unsigned int offset[2][2][2] = {
1564  {{4, 5}, // face_orientation=false; face_flip=false;
1565  // face_rotation=false and true
1566  {6, 7}}, // face_orientation=false; face_flip=true;
1567  // face_rotation=false and true
1568  {{0, 1}, // face_orientation=true; face_flip=false;
1569  // face_rotation=false and true
1570  {2, 3}}}; // face_orientation=true; face_flip=true;
1571  // face_rotation=false and true
1572 
1573 
1574  if (quadrature.size() == 1)
1575  return (face_no +
1576  offset[face_orientation][face_flip][face_rotation] *
1577  GeometryInfo<dim>::faces_per_cell) *
1578  quadrature[0].size();
1579  else
1580  {
1581  unsigned int n_points_i = 0;
1582  for (unsigned int i = 0; i < face_no; ++i)
1583  n_points_i += quadrature[i].size();
1584 
1585  unsigned int n_points = 0;
1586  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1587  ++i)
1588  n_points += quadrature[i].size();
1589 
1590  return (n_points_i +
1591  offset[face_orientation][face_flip][face_rotation] *
1592  n_points);
1593  }
1594  }
1595 
1596  default:
1597  Assert(false, ExcInternalError());
1598  }
1599  return numbers::invalid_unsigned_int;
1600 }
1601 
1602 
1603 
1604 template <>
1605 QProjector<1>::DataSetDescriptor
1606 QProjector<1>::DataSetDescriptor::subface(
1607  const ReferenceCell reference_cell,
1608  const unsigned int face_no,
1609  const unsigned int subface_no,
1610  const bool,
1611  const bool,
1612  const bool,
1613  const unsigned int n_quadrature_points,
1614  const internal::SubfaceCase<1>)
1615 {
1616  Assert(reference_cell == ReferenceCells::Line, ExcNotImplemented());
1617  (void)reference_cell;
1618 
1619  Assert(face_no < GeometryInfo<1>::faces_per_cell, ExcInternalError());
1620  Assert(subface_no < GeometryInfo<1>::max_children_per_face,
1621  ExcInternalError());
1622 
1623  return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1624  n_quadrature_points);
1625 }
1626 
1627 
1628 
1629 template <>
1630 QProjector<1>::DataSetDescriptor
1631 QProjector<1>::DataSetDescriptor::subface(
1632  const unsigned int face_no,
1633  const unsigned int subface_no,
1634  const bool face_orientation,
1635  const bool face_flip,
1636  const bool face_rotation,
1637  const unsigned int n_quadrature_points,
1638  const internal::SubfaceCase<1> ref_case)
1639 {
1640  return subface(ReferenceCells::Line,
1641  face_no,
1642  subface_no,
1643  face_orientation,
1644  face_flip,
1645  face_rotation,
1646  n_quadrature_points,
1647  ref_case);
1648 }
1649 
1650 
1651 
1652 template <>
1653 QProjector<2>::DataSetDescriptor
1654 QProjector<2>::DataSetDescriptor::subface(
1655  const ReferenceCell reference_cell,
1656  const unsigned int face_no,
1657  const unsigned int subface_no,
1658  const bool,
1659  const bool,
1660  const bool,
1661  const unsigned int n_quadrature_points,
1662  const internal::SubfaceCase<2>)
1663 {
1664  Assert(reference_cell == ReferenceCells::Quadrilateral, ExcNotImplemented());
1665  (void)reference_cell;
1666 
1667  Assert(face_no < GeometryInfo<2>::faces_per_cell, ExcInternalError());
1668  Assert(subface_no < GeometryInfo<2>::max_children_per_face,
1669  ExcInternalError());
1670 
1671  return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1672  n_quadrature_points);
1673 }
1674 
1675 
1676 
1677 template <>
1678 QProjector<2>::DataSetDescriptor
1679 QProjector<2>::DataSetDescriptor::subface(
1680  const unsigned int face_no,
1681  const unsigned int subface_no,
1682  const bool face_orientation,
1683  const bool face_flip,
1684  const bool face_rotation,
1685  const unsigned int n_quadrature_points,
1686  const internal::SubfaceCase<2> ref_case)
1687 {
1688  return subface(ReferenceCells::Quadrilateral,
1689  face_no,
1690  subface_no,
1691  face_orientation,
1692  face_flip,
1693  face_rotation,
1694  n_quadrature_points,
1695  ref_case);
1696 }
1697 
1698 
1699 template <>
1700 QProjector<3>::DataSetDescriptor
1701 QProjector<3>::DataSetDescriptor::subface(
1702  const ReferenceCell reference_cell,
1703  const unsigned int face_no,
1704  const unsigned int subface_no,
1705  const bool face_orientation,
1706  const bool face_flip,
1707  const bool face_rotation,
1708  const unsigned int n_quadrature_points,
1709  const internal::SubfaceCase<3> ref_case)
1710 {
1711  const unsigned int dim = 3;
1712 
1713  Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
1714  (void)reference_cell;
1715 
1716  Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
1717  Assert(subface_no < GeometryInfo<dim>::max_children_per_face,
1718  ExcInternalError());
1719 
1720  // As the quadrature points created by
1721  // QProjector are on subfaces in their
1722  // "standard location" we have to use a
1723  // permutation of the equivalent subface
1724  // number in order to respect face
1725  // orientation, flip and rotation. The
1726  // information we need here is exactly the
1727  // same as the
1728  // GeometryInfo<3>::child_cell_on_face info
1729  // for the bottom face (face 4) of a hex, as
1730  // on this the RefineCase of the cell matches
1731  // that of the face and the subfaces are
1732  // numbered in the same way as the child
1733  // cells.
1734 
1735  // in 3d, we have to account for faces that
1736  // have non-standard face orientation, flip
1737  // and rotation. thus, we have to store
1738  // _eight_ data sets per face or subface
1739  // already for the isotropic
1740  // case. Additionally, we have three
1741  // different refinement cases, resulting in
1742  // <tt>4 + 2 + 2 = 8</tt> different subfaces
1743  // for each face.
1744  const unsigned int total_subfaces_per_face = 8;
1745 
1746  // set up a table with the according offsets
1747  // for non-standard orientation, first index:
1748  // face_orientation (standard true=1), second
1749  // index: face_flip (standard false=0), third
1750  // index: face_rotation (standard false=0)
1751  //
1752  // note, that normally we should use the
1753  // obvious offsets 0,1,2,3,4,5,6,7. However,
1754  // prior to the changes enabling flipped and
1755  // rotated faces, in many places of the
1756  // library the convention was used, that the
1757  // first dataset with offset 0 corresponds to
1758  // a face in standard orientation. therefore
1759  // we use the offsets 4,5,6,7,0,1,2,3 here to
1760  // stick to that (implicit) convention
1761  static const unsigned int orientation_offset[2][2][2] = {
1762  {// face_orientation=false; face_flip=false; face_rotation=false and true
1763  {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1764  5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1765  // face_orientation=false; face_flip=true; face_rotation=false and true
1766  {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1767  7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
1768  {// face_orientation=true; face_flip=false; face_rotation=false and true
1769  {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1770  1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1771  // face_orientation=true; face_flip=true; face_rotation=false and true
1772  {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1773  3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};
1774 
1775  // set up a table with the offsets for a
1776  // given refinement case respecting the
1777  // corresponding number of subfaces. the
1778  // index corresponds to (RefineCase::Type - 1)
1779 
1780  // note, that normally we should use the
1781  // obvious offsets 0,2,6. However, prior to
1782  // the implementation of anisotropic
1783  // refinement, in many places of the library
1784  // the convention was used, that the first
1785  // dataset with offset 0 corresponds to a
1786  // standard (isotropic) face
1787  // refinement. therefore we use the offsets
1788  // 6,4,0 here to stick to that (implicit)
1789  // convention
1790  static const unsigned int ref_case_offset[3] = {
1791  6, // cut_x
1792  4, // cut_y
1793  0 // cut_xy
1794  };
1795 
1796 
1797  // for each subface of a given FaceRefineCase
1798  // there is a corresponding equivalent
1799  // subface number of one of the "standard"
1800  // RefineCases (cut_x, cut_y, cut_xy). Map
1801  // the given values to those equivalent
1802  // ones.
1803 
1804  // first, define an invalid number
1805  static const unsigned int e = numbers::invalid_unsigned_int;
1806 
1807  static const RefinementCase<dim - 1>
1808  equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
1809  [GeometryInfo<3>::max_children_per_face] = {
1810  // case_none. there should be only
1811  // invalid values here. However, as
1812  // this function is also called (in
1813  // tests) for cells which have no
1814  // refined faces, use isotropic
1815  // refinement instead
1816  {RefinementCase<dim - 1>::cut_xy,
1817  RefinementCase<dim - 1>::cut_xy,
1818  RefinementCase<dim - 1>::cut_xy,
1819  RefinementCase<dim - 1>::cut_xy},
1820  // case_x
1821  {RefinementCase<dim - 1>::cut_x,
1822  RefinementCase<dim - 1>::cut_x,
1823  RefinementCase<dim - 1>::no_refinement,
1824  RefinementCase<dim - 1>::no_refinement},
1825  // case_x1y
1826  {RefinementCase<dim - 1>::cut_xy,
1827  RefinementCase<dim - 1>::cut_xy,
1828  RefinementCase<dim - 1>::cut_x,
1829  RefinementCase<dim - 1>::no_refinement},
1830  // case_x2y
1831  {RefinementCase<dim - 1>::cut_x,
1832  RefinementCase<dim - 1>::cut_xy,
1833  RefinementCase<dim - 1>::cut_xy,
1834  RefinementCase<dim - 1>::no_refinement},
1835  // case_x1y2y
1836  {RefinementCase<dim - 1>::cut_xy,
1837  RefinementCase<dim - 1>::cut_xy,
1838  RefinementCase<dim - 1>::cut_xy,
1839  RefinementCase<dim - 1>::cut_xy},
1840  // case_y
1841  {RefinementCase<dim - 1>::cut_y,
1842  RefinementCase<dim - 1>::cut_y,
1843  RefinementCase<dim - 1>::no_refinement,
1844  RefinementCase<dim - 1>::no_refinement},
1845  // case_y1x
1846  {RefinementCase<dim - 1>::cut_xy,
1847  RefinementCase<dim - 1>::cut_xy,
1848  RefinementCase<dim - 1>::cut_y,
1849  RefinementCase<dim - 1>::no_refinement},
1850  // case_y2x
1851  {RefinementCase<dim - 1>::cut_y,
1852  RefinementCase<dim - 1>::cut_xy,
1853  RefinementCase<dim - 1>::cut_xy,
1854  RefinementCase<dim - 1>::no_refinement},
1855  // case_y1x2x
1856  {RefinementCase<dim - 1>::cut_xy,
1857  RefinementCase<dim - 1>::cut_xy,
1858  RefinementCase<dim - 1>::cut_xy,
1859  RefinementCase<dim - 1>::cut_xy},
1860  // case_xy (case_isotropic)
1861  {RefinementCase<dim - 1>::cut_xy,
1862  RefinementCase<dim - 1>::cut_xy,
1863  RefinementCase<dim - 1>::cut_xy,
1864  RefinementCase<dim - 1>::cut_xy}};
1865 
1866  static const unsigned int
1867  equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
1868  [GeometryInfo<3>::max_children_per_face] = {
1869  // case_none, see above
1870  {0, 1, 2, 3},
1871  // case_x
1872  {0, 1, e, e},
1873  // case_x1y
1874  {0, 2, 1, e},
1875  // case_x2y
1876  {0, 1, 3, e},
1877  // case_x1y2y
1878  {0, 2, 1, 3},
1879  // case_y
1880  {0, 1, e, e},
1881  // case_y1x
1882  {0, 1, 1, e},
1883  // case_y2x
1884  {0, 2, 3, e},
1885  // case_y1x2x
1886  {0, 1, 2, 3},
1887  // case_xy (case_isotropic)
1888  {0, 1, 2, 3}};
1889 
1890  // If face-orientation or face_rotation are
1891  // non-standard, cut_x and cut_y have to be
1892  // exchanged.
1893  static const RefinementCase<dim - 1> ref_case_permutation[4] = {
1894  RefinementCase<dim - 1>::no_refinement,
1895  RefinementCase<dim - 1>::cut_y,
1896  RefinementCase<dim - 1>::cut_x,
1897  RefinementCase<dim - 1>::cut_xy};
1898 
1899  // set a corresponding (equivalent)
1900  // RefineCase and subface number
1901  const RefinementCase<dim - 1> equ_ref_case =
1902  equivalent_refine_case[ref_case][subface_no];
1903  const unsigned int equ_subface_no =
1904  equivalent_subface_number[ref_case][subface_no];
1905  // make sure, that we got a valid subface and RefineCase
1906  Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
1907  ExcInternalError());
1908  Assert(equ_subface_no != e, ExcInternalError());
1909  // now, finally respect non-standard faces
1910  const RefinementCase<dim - 1> final_ref_case =
1911  (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1912  equ_ref_case);
1913 
1914  // what we have now is the number of
1915  // the subface in the natural
1916  // orientation of the *face*. what we
1917  // need to know is the number of the
1918  // subface concerning the standard face
1919  // orientation as seen from the *cell*.
1920 
1921  // this mapping is not trivial, but we
1922  // have done exactly this stuff in the
1923  // child_cell_on_face function. in
1924  // order to reduce the amount of code
1925  // as well as to make maintaining the
1926  // functionality easier we want to
1927  // reuse that information. So we note
1928  // that on the bottom face (face 4) of
1929  // a hex cell the local x and y
1930  // coordinates of the face and the cell
1931  // coincide, thus also the refinement
1932  // case of the face corresponds to the
1933  // refinement case of the cell
1934  // (ignoring cell refinement along the
1935  // z direction). Using this knowledge
1936  // we can (ab)use the
1937  // child_cell_on_face function to do
1938  // exactly the transformation we are in
1939  // need of now
1940  const unsigned int final_subface_no =
1941  GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(final_ref_case),
1942  4,
1943  equ_subface_no,
1944  face_orientation,
1945  face_flip,
1946  face_rotation,
1947  equ_ref_case);
1948 
1949  return (((face_no * total_subfaces_per_face +
1950  ref_case_offset[final_ref_case - 1] + final_subface_no) +
1951  orientation_offset[face_orientation][face_flip][face_rotation]) *
1952  n_quadrature_points);
1953 }
1954 
1955 
1956 template <>
1957 QProjector<3>::DataSetDescriptor
1958 QProjector<3>::DataSetDescriptor::subface(
1959  const unsigned int face_no,
1960  const unsigned int subface_no,
1961  const bool face_orientation,
1962  const bool face_flip,
1963  const bool face_rotation,
1964  const unsigned int n_quadrature_points,
1965  const internal::SubfaceCase<3> ref_case)
1966 {
1967  return subface(ReferenceCells::Hexahedron,
1968  face_no,
1969  subface_no,
1970  face_orientation,
1971  face_flip,
1972  face_rotation,
1973  n_quadrature_points,
1974  ref_case);
1975 }
1976 
1977 
1978 
1979 template <int dim>
1980 Quadrature<dim>
1981 QProjector<dim>::project_to_face(const SubQuadrature &quadrature,
1982  const unsigned int face_no)
1983 {
1984  return project_to_face(ReferenceCells::get_hypercube<dim>(),
1985  quadrature,
1986  face_no);
1987 }
1988 
1989 
1990 
1991 template <int dim>
1992 Quadrature<dim>
1993 QProjector<dim>::project_to_face(const ReferenceCell reference_cell,
1994  const SubQuadrature &quadrature,
1995  const unsigned int face_no)
1996 {
1997  Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
1998  ExcNotImplemented());
1999  (void)reference_cell;
2000 
2001  std::vector<Point<dim>> points(quadrature.size());
2002  project_to_face(quadrature, face_no, points);
2003  return Quadrature<dim>(points, quadrature.get_weights());
2004 }
2005 
2006 
2007 
2008 template <int dim>
2009 Quadrature<dim>
2010 QProjector<dim>::project_to_subface(const SubQuadrature & quadrature,
2011  const unsigned int face_no,
2012  const unsigned int subface_no,
2013  const RefinementCase<dim - 1> &ref_case)
2014 {
2015  return project_to_subface(ReferenceCells::get_hypercube<dim>(),
2016  quadrature,
2017  face_no,
2018  subface_no,
2019  ref_case);
2020 }
2021 
2022 
2023 
2024 template <int dim>
2025 Quadrature<dim>
2026 QProjector<dim>::project_to_subface(const ReferenceCell reference_cell,
2027  const SubQuadrature &quadrature,
2028  const unsigned int face_no,
2029  const unsigned int subface_no,
2030  const RefinementCase<dim - 1> &ref_case)
2031 {
2032  Assert(reference_cell == ReferenceCells::get_hypercube<dim>(),
2033  ExcNotImplemented());
2034  (void)reference_cell;
2035 
2036  std::vector<Point<dim>> points(quadrature.size());
2037  project_to_subface(quadrature, face_no, subface_no, points, ref_case);
2038  return Quadrature<dim>(points, quadrature.get_weights());
2039 }
2040 
2041 
2042 // explicit instantiations; note: we need them all for all dimensions
2043 template class QProjector<1>;
2044 template class QProjector<2>;
2045 template class QProjector<3>;
2046 
2047 DEAL_II_NAMESPACE_CLOSE
BarycentricPolynomials::get_fe_p_basis
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition: polynomials_barycentric.cc:58
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:2001
SymmetricTensor::determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2622
VectorizedArray::sqrt
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5406
hp::Collection::size
unsigned int size() const
Definition: collection.h:109
hp::QCollection::max_n_quadrature_points
unsigned int max_n_quadrature_points() const
Definition: q_collection.h:181
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:396
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:397
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:175
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:1465
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:2042
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:701
ReferenceCells::Tetrahedron
constexpr const ReferenceCell Tetrahedron
Definition: reference_cell.h:644
ReferenceCells::Quadrilateral
constexpr const ReferenceCell Quadrilateral
Definition: reference_cell.h:642
ReferenceCells::Wedge
constexpr const ReferenceCell Wedge
Definition: reference_cell.h:648
ReferenceCells::Pyramid
constexpr const ReferenceCell Pyramid
Definition: reference_cell.h:646
ReferenceCells::Triangle
constexpr const ReferenceCell Triangle
Definition: reference_cell.h:640
ReferenceCells::Hexahedron
constexpr const ReferenceCell Hexahedron
Definition: reference_cell.h:650
ReferenceCells::Line
constexpr const ReferenceCell Line
Definition: reference_cell.h:638
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:196
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