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