Reference documentation for deal.II version Git c336d204ee 2020-06-05 23:44:40 -0400
\(\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\}}\)
quadrature.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
20 #include <deal.II/base/utilities.h>
21 
22 #include <cmath>
23 #include <cstdlib>
24 #include <iterator>
25 #include <memory>
26 
28 
29 
30 template <>
31 Quadrature<0>::Quadrature(const unsigned int n_q)
32  : quadrature_points(n_q)
33  , weights(n_q, 0)
34  , is_tensor_product_flag(false)
35 {}
36 
37 
38 
39 template <int dim>
40 Quadrature<dim>::Quadrature(const unsigned int n_q)
41  : quadrature_points(n_q, Point<dim>())
42  , weights(n_q, 0)
43  , is_tensor_product_flag(dim == 1)
44 {}
45 
46 
47 
48 template <int dim>
49 void
51  const std::vector<double> & w)
52 {
53  AssertDimension(w.size(), p.size());
55  weights = w;
56 }
57 
58 
59 
60 template <int dim>
61 Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points,
62  const std::vector<double> & weights)
63  : quadrature_points(points)
64  , weights(weights)
65  , is_tensor_product_flag(dim == 1)
66 {
67  Assert(weights.size() == points.size(),
68  ExcDimensionMismatch(weights.size(), points.size()));
69 }
70 
71 
72 
73 template <int dim>
74 Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points)
75  : quadrature_points(points)
76  , weights(points.size(), std::numeric_limits<double>::infinity())
77  , is_tensor_product_flag(dim == 1)
78 {
79  Assert(weights.size() == points.size(),
80  ExcDimensionMismatch(weights.size(), points.size()));
81 }
82 
83 
84 
85 template <int dim>
87  : quadrature_points(std::vector<Point<dim>>(1, point))
88  , weights(std::vector<double>(1, 1.))
90  , tensor_basis(new std::array<Quadrature<1>, dim>())
91 {
92  for (unsigned int i = 0; i < dim; ++i)
93  {
94  const std::vector<Point<1>> quad_vec_1d(1, Point<1>(point[i]));
95  (*tensor_basis)[i] = Quadrature<1>(quad_vec_1d, weights);
96  }
97 }
98 
99 
100 #ifndef DOXYGEN
101 template <>
103  : quadrature_points(std::vector<Point<1>>(1, point))
104  , weights(std::vector<double>(1, 1.))
105  , is_tensor_product_flag(true)
106 {}
107 
108 
109 
110 template <>
112 {
113  Assert(false, ExcImpossibleInDim(0));
114 }
115 #endif // DOXYGEN
116 
117 
118 
119 template <int dim>
121  : quadrature_points(q1.size() * q2.size())
122  , weights(q1.size() * q2.size())
124 {
125  unsigned int present_index = 0;
126  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
127  for (unsigned int i1 = 0; i1 < q1.size(); ++i1)
128  {
129  // compose coordinates of
130  // new quadrature point by tensor
131  // product in the last component
132  for (unsigned int d = 0; d < dim - 1; ++d)
133  quadrature_points[present_index](d) = q1.point(i1)(d);
134  quadrature_points[present_index](dim - 1) = q2.point(i2)(0);
135 
136  weights[present_index] = q1.weight(i1) * q2.weight(i2);
137 
138  ++present_index;
139  }
140 
141 #ifdef DEBUG
142  if (size() > 0)
143  {
144  double sum = 0;
145  for (unsigned int i = 0; i < size(); ++i)
146  sum += weights[i];
147  // we cannot guarantee the sum of weights
148  // to be exactly one, but it should be
149  // near that.
150  Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
151  }
152 #endif
153 
155  {
156  tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
157  for (unsigned int i = 0; i < dim - 1; ++i)
158  (*tensor_basis)[i] = q1.get_tensor_basis()[i];
159  (*tensor_basis)[dim - 1] = q2;
160  }
161 }
162 
163 
164 #ifndef DOXYGEN
165 template <>
167  : quadrature_points(q2.size())
168  , weights(q2.size())
169  , is_tensor_product_flag(true)
170 {
171  unsigned int present_index = 0;
172  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
173  {
174  // compose coordinates of
175  // new quadrature point by tensor
176  // product in the last component
177  quadrature_points[present_index](0) = q2.point(i2)(0);
178 
179  weights[present_index] = q2.weight(i2);
180 
181  ++present_index;
182  }
183 
184 # ifdef DEBUG
185  if (size() > 0)
186  {
187  double sum = 0;
188  for (unsigned int i = 0; i < size(); ++i)
189  sum += weights[i];
190  // we cannot guarantee the sum of weights
191  // to be exactly one, but it should be
192  // near that.
193  Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
194  }
195 # endif
196 }
197 
198 
199 
200 template <>
202  : Subscriptor()
203  ,
204  // quadrature_points(1),
205  weights(1, 1.)
206  , is_tensor_product_flag(false)
207 {}
208 
209 
210 template <>
212  : Subscriptor()
213 {
214  // this function should never be
215  // called -- this should be the
216  // copy constructor in 1d...
217  Assert(false, ExcImpossibleInDim(1));
218 }
219 #endif // DOXYGEN
220 
221 
222 
223 template <int dim>
225  : Subscriptor()
227  , weights(Utilities::fixed_power<dim>(q.size()))
228  , is_tensor_product_flag(true)
229 {
230  Assert(dim <= 3, ExcNotImplemented());
231 
232  const unsigned int n0 = q.size();
233  const unsigned int n1 = (dim > 1) ? n0 : 1;
234  const unsigned int n2 = (dim > 2) ? n0 : 1;
235 
236  unsigned int k = 0;
237  for (unsigned int i2 = 0; i2 < n2; ++i2)
238  for (unsigned int i1 = 0; i1 < n1; ++i1)
239  for (unsigned int i0 = 0; i0 < n0; ++i0)
240  {
241  quadrature_points[k](0) = q.point(i0)(0);
242  if (dim > 1)
243  quadrature_points[k](1) = q.point(i1)(0);
244  if (dim > 2)
245  quadrature_points[k](2) = q.point(i2)(0);
246  weights[k] = q.weight(i0);
247  if (dim > 1)
248  weights[k] *= q.weight(i1);
249  if (dim > 2)
250  weights[k] *= q.weight(i2);
251  ++k;
252  }
253 
254  tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
255  for (unsigned int i = 0; i < dim; ++i)
256  (*tensor_basis)[i] = q;
257 }
258 
259 
260 
261 template <int dim>
263  : Subscriptor()
265  , weights(q.weights)
267 {
268  if (dim > 1 && is_tensor_product_flag)
269  tensor_basis =
270  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
271 }
272 
273 
274 
275 template <int dim>
278 {
279  weights = q.weights;
282  if (dim > 1 && is_tensor_product_flag)
283  {
284  if (tensor_basis == nullptr)
285  tensor_basis =
286  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
287  else
289  }
290  return *this;
291 }
292 
293 
294 
295 template <int dim>
296 bool
298 {
299  return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
300 }
301 
302 
303 
304 template <int dim>
305 std::size_t
307 {
310 }
311 
312 
313 
314 template <int dim>
315 typename std::conditional<dim == 1,
316  std::array<Quadrature<1>, dim>,
317  const std::array<Quadrature<1>, dim> &>::type
319 {
320  Assert(this->is_tensor_product_flag == true,
321  ExcMessage("This function only makes sense if "
322  "this object represents a tensor product!"));
323  Assert(tensor_basis != nullptr, ExcInternalError());
324 
325  return *tensor_basis;
326 }
327 
328 
329 #ifndef DOXYGEN
330 template <>
331 std::array<Quadrature<1>, 1>
333 {
334  Assert(this->is_tensor_product_flag == true,
335  ExcMessage("This function only makes sense if "
336  "this object represents a tensor product!"));
337 
338  return std::array<Quadrature<1>, 1>{{*this}};
339 }
340 #endif
341 
342 
343 
344 //---------------------------------------------------------------------------
345 template <int dim>
347  : Quadrature<dim>(qx.size())
348 {
349  Assert(dim == 1, ExcImpossibleInDim(dim));
350  unsigned int k = 0;
351  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
352  {
353  this->quadrature_points[k](0) = qx.point(k1)(0);
354  this->weights[k++] = qx.weight(k1);
355  }
356  Assert(k == this->size(), ExcInternalError());
357  this->is_tensor_product_flag = true;
358 }
359 
360 
361 
362 template <int dim>
364  const Quadrature<1> &qy)
365  : Quadrature<dim>(qx.size() * qy.size())
366 {
367  Assert(dim == 2, ExcImpossibleInDim(dim));
368 }
369 
370 
371 
372 template <>
374  : Quadrature<2>(qx.size() * qy.size())
375 {
376  unsigned int k = 0;
377  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
378  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
379  {
380  this->quadrature_points[k](0) = qx.point(k1)(0);
381  this->quadrature_points[k](1) = qy.point(k2)(0);
382  this->weights[k++] = qx.weight(k1) * qy.weight(k2);
383  }
384  Assert(k == this->size(), ExcInternalError());
385  this->is_tensor_product_flag = true;
386  const std::array<Quadrature<1>, 2> q_array{{qx, qy}};
387  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 2>>(q_array);
388 }
389 
390 
391 
392 template <int dim>
394  const Quadrature<1> &qy,
395  const Quadrature<1> &qz)
396  : Quadrature<dim>(qx.size() * qy.size() * qz.size())
397 {
398  Assert(dim == 3, ExcImpossibleInDim(dim));
399 }
400 
401 
402 
403 template <>
405  const Quadrature<1> &qy,
406  const Quadrature<1> &qz)
407  : Quadrature<3>(qx.size() * qy.size() * qz.size())
408 {
409  unsigned int k = 0;
410  for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
411  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
412  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
413  {
414  this->quadrature_points[k](0) = qx.point(k1)(0);
415  this->quadrature_points[k](1) = qy.point(k2)(0);
416  this->quadrature_points[k](2) = qz.point(k3)(0);
417  this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
418  }
419  Assert(k == this->size(), ExcInternalError());
420  this->is_tensor_product_flag = true;
421  const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
422  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 3>>(q_array);
423 }
424 
425 
426 
427 //---------------------------------------------------------------------------
428 
429 
430 
431 template <int dim>
434 {
435  std::vector<Point<2>> q_points(q.size());
436  std::vector<double> weights(q.size());
437  for (unsigned int i = 0; i < q.size(); ++i)
438  {
439  q_points[i][0] = q.point(i)[1];
440  q_points[i][1] = q.point(i)[0];
441 
442  weights[i] = q.weight(i);
443  }
444 
445  return Quadrature<2>(q_points, weights);
446 }
447 
448 
449 template <int dim>
451 QProjector<dim>::rotate(const Quadrature<2> &q, const unsigned int n_times)
452 {
453  std::vector<Point<2>> q_points(q.size());
454  std::vector<double> weights(q.size());
455  for (unsigned int i = 0; i < q.size(); ++i)
456  {
457  switch (n_times % 4)
458  {
459  case 0:
460  // 0 degree
461  q_points[i][0] = q.point(i)[0];
462  q_points[i][1] = q.point(i)[1];
463  break;
464  case 1:
465  // 90 degree counterclockwise
466  q_points[i][0] = 1.0 - q.point(i)[1];
467  q_points[i][1] = q.point(i)[0];
468  break;
469  case 2:
470  // 180 degree counterclockwise
471  q_points[i][0] = 1.0 - q.point(i)[0];
472  q_points[i][1] = 1.0 - q.point(i)[1];
473  break;
474  case 3:
475  // 270 degree counterclockwise
476  q_points[i][0] = q.point(i)[1];
477  q_points[i][1] = 1.0 - q.point(i)[0];
478  break;
479  }
480 
481  weights[i] = q.weight(i);
482  }
483 
484  return Quadrature<2>(q_points, weights);
485 }
486 
487 
488 template <>
489 void
491  const unsigned int face_no,
492  std::vector<Point<1>> &q_points)
493 {
494  const unsigned int dim = 1;
496  AssertDimension(q_points.size(), 1);
497 
498  q_points[0] = Point<dim>(static_cast<double>(face_no));
499 }
500 
501 
502 
503 template <>
504 void
506  const unsigned int face_no,
507  std::vector<Point<2>> &q_points)
508 {
509  const unsigned int dim = 2;
511  Assert(q_points.size() == quadrature.size(),
512  ExcDimensionMismatch(q_points.size(), quadrature.size()));
513 
514  for (unsigned int p = 0; p < quadrature.size(); ++p)
515  switch (face_no)
516  {
517  case 0:
518  q_points[p] = Point<dim>(0, quadrature.point(p)(0));
519  break;
520  case 1:
521  q_points[p] = Point<dim>(1, quadrature.point(p)(0));
522  break;
523  case 2:
524  q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
525  break;
526  case 3:
527  q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
528  break;
529  default:
530  Assert(false, ExcInternalError());
531  }
532 }
533 
534 
535 
536 template <>
537 void
539  const unsigned int face_no,
540  std::vector<Point<3>> &q_points)
541 {
542  const unsigned int dim = 3;
544  Assert(q_points.size() == quadrature.size(),
545  ExcDimensionMismatch(q_points.size(), quadrature.size()));
546 
547  for (unsigned int p = 0; p < quadrature.size(); ++p)
548  switch (face_no)
549  {
550  case 0:
551  q_points[p] =
552  Point<dim>(0, quadrature.point(p)(0), quadrature.point(p)(1));
553  break;
554  case 1:
555  q_points[p] =
556  Point<dim>(1, quadrature.point(p)(0), quadrature.point(p)(1));
557  break;
558  case 2:
559  q_points[p] =
560  Point<dim>(quadrature.point(p)(1), 0, quadrature.point(p)(0));
561  break;
562  case 3:
563  q_points[p] =
564  Point<dim>(quadrature.point(p)(1), 1, quadrature.point(p)(0));
565  break;
566  case 4:
567  q_points[p] =
568  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 0);
569  break;
570  case 5:
571  q_points[p] =
572  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 1);
573  break;
574 
575  default:
576  Assert(false, ExcInternalError());
577  }
578 }
579 
580 
581 
582 template <>
583 void
585  const unsigned int face_no,
586  const unsigned int,
587  std::vector<Point<1>> &q_points,
588  const RefinementCase<0> &)
589 {
590  const unsigned int dim = 1;
592  AssertDimension(q_points.size(), 1);
593 
594  q_points[0] = Point<dim>(static_cast<double>(face_no));
595 }
596 
597 
598 
599 template <>
600 void
602  const unsigned int face_no,
603  const unsigned int subface_no,
604  std::vector<Point<2>> &q_points,
605  const RefinementCase<1> &)
606 {
607  const unsigned int dim = 2;
610 
611  Assert(q_points.size() == quadrature.size(),
612  ExcDimensionMismatch(q_points.size(), quadrature.size()));
613 
614  for (unsigned int p = 0; p < quadrature.size(); ++p)
615  switch (face_no)
616  {
617  case 0:
618  switch (subface_no)
619  {
620  case 0:
621  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
622  break;
623  case 1:
624  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
625  break;
626  default:
627  Assert(false, ExcInternalError());
628  }
629  break;
630  case 1:
631  switch (subface_no)
632  {
633  case 0:
634  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
635  break;
636  case 1:
637  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
638  break;
639  default:
640  Assert(false, ExcInternalError());
641  }
642  break;
643  case 2:
644  switch (subface_no)
645  {
646  case 0:
647  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
648  break;
649  case 1:
650  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
651  break;
652  default:
653  Assert(false, ExcInternalError());
654  }
655  break;
656  case 3:
657  switch (subface_no)
658  {
659  case 0:
660  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
661  break;
662  case 1:
663  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
664  break;
665  default:
666  Assert(false, ExcInternalError());
667  }
668  break;
669 
670  default:
671  Assert(false, ExcInternalError());
672  }
673 }
674 
675 
676 
677 template <>
678 void
680  const unsigned int face_no,
681  const unsigned int subface_no,
682  std::vector<Point<3>> & q_points,
683  const RefinementCase<2> &ref_case)
684 {
685  const unsigned int dim = 3;
688  Assert(q_points.size() == quadrature.size(),
689  ExcDimensionMismatch(q_points.size(), quadrature.size()));
690 
691  // one coordinate is at a const value. for
692  // faces 0, 2 and 4 this value is 0.0, for
693  // faces 1, 3 and 5 it is 1.0
694  double const_value = face_no % 2;
695  // local 2d coordinates are xi and eta,
696  // global 3d coordinates are x, y and
697  // z. those have to be mapped. the following
698  // indices tell, which global coordinate
699  // (0->x, 1->y, 2->z) corresponds to which
700  // local one
701  unsigned int xi_index = numbers::invalid_unsigned_int,
702  eta_index = numbers::invalid_unsigned_int,
703  const_index = face_no / 2;
704  // the xi and eta values have to be scaled
705  // (by factor 0.5 or factor 1.0) depending on
706  // the refinement case and translated (by 0.0
707  // or 0.5) depending on the refinement case
708  // and subface_no.
709  double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
710  eta_translation = 0.0;
711  // set the index mapping between local and
712  // global coordinates
713  switch (face_no / 2)
714  {
715  case 0:
716  xi_index = 1;
717  eta_index = 2;
718  break;
719  case 1:
720  xi_index = 2;
721  eta_index = 0;
722  break;
723  case 2:
724  xi_index = 0;
725  eta_index = 1;
726  break;
727  }
728  // set the scale and translation parameter
729  // for individual subfaces
730  switch (ref_case)
731  {
732  case RefinementCase<dim - 1>::cut_x:
733  xi_scale = 0.5;
734  xi_translation = subface_no % 2 * 0.5;
735  break;
736  case RefinementCase<dim - 1>::cut_y:
737  eta_scale = 0.5;
738  eta_translation = subface_no % 2 * 0.5;
739  break;
740  case RefinementCase<dim - 1>::cut_xy:
741  xi_scale = 0.5;
742  eta_scale = 0.5;
743  xi_translation = int(subface_no % 2) * 0.5;
744  eta_translation = int(subface_no / 2) * 0.5;
745  break;
746  default:
747  Assert(false, ExcInternalError());
748  break;
749  }
750  // finally, compute the scaled, translated,
751  // projected quadrature points
752  for (unsigned int p = 0; p < quadrature.size(); ++p)
753  {
754  q_points[p][xi_index] =
755  xi_scale * quadrature.point(p)(0) + xi_translation;
756  q_points[p][eta_index] =
757  eta_scale * quadrature.point(p)(1) + eta_translation;
758  q_points[p][const_index] = const_value;
759  }
760 }
761 
762 
763 template <>
766 {
767  const unsigned int dim = 1;
768 
769  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;
770 
771  // first fix quadrature points
772  std::vector<Point<dim>> q_points;
773  q_points.reserve(n_points * n_faces);
774  std::vector<Point<dim>> help(n_points);
775 
776 
777  // project to each face and append
778  // results
779  for (unsigned int face = 0; face < n_faces; ++face)
780  {
781  project_to_face(quadrature, face, help);
782  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
783  }
784 
785  // next copy over weights
786  std::vector<double> weights;
787  weights.reserve(n_points * n_faces);
788  for (unsigned int face = 0; face < n_faces; ++face)
789  std::copy(quadrature.get_weights().begin(),
790  quadrature.get_weights().end(),
791  std::back_inserter(weights));
792 
793  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
794  Assert(weights.size() == n_points * n_faces, ExcInternalError());
795 
796  return Quadrature<dim>(q_points, weights);
797 }
798 
799 
800 
801 template <>
804 {
805  const unsigned int dim = 2;
806 
807  const unsigned int n_points = quadrature.size(),
809 
810  // first fix quadrature points
811  std::vector<Point<dim>> q_points;
812  q_points.reserve(n_points * n_faces);
813  std::vector<Point<dim>> help(n_points);
814 
815  // project to each face and append
816  // results
817  for (unsigned int face = 0; face < n_faces; ++face)
818  {
819  project_to_face(quadrature, face, help);
820  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
821  }
822 
823  // next copy over weights
824  std::vector<double> weights;
825  weights.reserve(n_points * n_faces);
826  for (unsigned int face = 0; face < n_faces; ++face)
827  std::copy(quadrature.get_weights().begin(),
828  quadrature.get_weights().end(),
829  std::back_inserter(weights));
830 
831  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
832  Assert(weights.size() == n_points * n_faces, ExcInternalError());
833 
834  return Quadrature<dim>(q_points, weights);
835 }
836 
837 
838 
839 template <>
842 {
843  const unsigned int dim = 3;
844 
845  SubQuadrature q_reflected = reflect(quadrature);
846  SubQuadrature q[8] = {quadrature,
847  rotate(quadrature, 1),
848  rotate(quadrature, 2),
849  rotate(quadrature, 3),
850  q_reflected,
851  rotate(q_reflected, 3),
852  rotate(q_reflected, 2),
853  rotate(q_reflected, 1)};
854 
855 
856 
857  const unsigned int n_points = quadrature.size(),
859 
860  // first fix quadrature points
861  std::vector<Point<dim>> q_points;
862  q_points.reserve(n_points * n_faces * 8);
863  std::vector<Point<dim>> help(n_points);
864 
865  std::vector<double> weights;
866  weights.reserve(n_points * n_faces * 8);
867 
868  // do the following for all possible
869  // mutations of a face (mutation==0
870  // corresponds to a face with standard
871  // orientation, no flip and no rotation)
872  for (const auto &mutation : q)
873  {
874  // project to each face and append
875  // results
876  for (unsigned int face = 0; face < n_faces; ++face)
877  {
878  project_to_face(mutation, face, help);
879  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
880  }
881 
882  // next copy over weights
883  for (unsigned int face = 0; face < n_faces; ++face)
884  std::copy(mutation.get_weights().begin(),
885  mutation.get_weights().end(),
886  std::back_inserter(weights));
887  }
888 
889 
890  Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError());
891  Assert(weights.size() == n_points * n_faces * 8, ExcInternalError());
892 
893  return Quadrature<dim>(q_points, weights);
894 }
895 
896 
897 
898 template <>
901 {
902  const unsigned int dim = 1;
903 
904  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
905  subfaces_per_face =
907 
908  // first fix quadrature points
909  std::vector<Point<dim>> q_points;
910  q_points.reserve(n_points * n_faces * subfaces_per_face);
911  std::vector<Point<dim>> help(n_points);
912 
913  // project to each face and copy
914  // results
915  for (unsigned int face = 0; face < n_faces; ++face)
916  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
917  {
918  project_to_subface(quadrature, face, subface, help);
919  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
920  }
921 
922  // next copy over weights
923  std::vector<double> weights;
924  weights.reserve(n_points * n_faces * subfaces_per_face);
925  for (unsigned int face = 0; face < n_faces; ++face)
926  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
927  std::copy(quadrature.get_weights().begin(),
928  quadrature.get_weights().end(),
929  std::back_inserter(weights));
930 
931  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
932  ExcInternalError());
933  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
934  ExcInternalError());
935 
936  return Quadrature<dim>(q_points, weights);
937 }
938 
939 
940 
941 template <>
944 {
945  const unsigned int dim = 2;
946 
947  const unsigned int n_points = quadrature.size(),
949  subfaces_per_face =
951 
952  // first fix quadrature points
953  std::vector<Point<dim>> q_points;
954  q_points.reserve(n_points * n_faces * subfaces_per_face);
955  std::vector<Point<dim>> help(n_points);
956 
957  // project to each face and copy
958  // results
959  for (unsigned int face = 0; face < n_faces; ++face)
960  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
961  {
962  project_to_subface(quadrature, face, subface, help);
963  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
964  }
965 
966  // next copy over weights
967  std::vector<double> weights;
968  weights.reserve(n_points * n_faces * subfaces_per_face);
969  for (unsigned int face = 0; face < n_faces; ++face)
970  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
971  std::copy(quadrature.get_weights().begin(),
972  quadrature.get_weights().end(),
973  std::back_inserter(weights));
974 
975  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
976  ExcInternalError());
977  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
978  ExcInternalError());
979 
980  return Quadrature<dim>(q_points, weights);
981 }
982 
983 
984 
985 template <>
988 {
989  const unsigned int dim = 3;
990  SubQuadrature q_reflected = reflect(quadrature);
991  SubQuadrature q[8] = {quadrature,
992  rotate(quadrature, 1),
993  rotate(quadrature, 2),
994  rotate(quadrature, 3),
995  q_reflected,
996  rotate(q_reflected, 3),
997  rotate(q_reflected, 2),
998  rotate(q_reflected, 1)};
999 
1000  const unsigned int n_points = quadrature.size(),
1002  total_subfaces_per_face = 2 + 2 + 4;
1003 
1004  // first fix quadrature points
1005  std::vector<Point<dim>> q_points;
1006  q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1007  std::vector<Point<dim>> help(n_points);
1008 
1009  std::vector<double> weights;
1010  weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1011 
1012  // do the following for all possible
1013  // mutations of a face (mutation==0
1014  // corresponds to a face with standard
1015  // orientation, no flip and no rotation)
1016  for (const auto &mutation : q)
1017  {
1018  // project to each face and copy
1019  // results
1020  for (unsigned int face = 0; face < n_faces; ++face)
1021  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1022  ref_case >= RefinementCase<dim - 1>::cut_x;
1023  --ref_case)
1024  for (unsigned int subface = 0;
1026  RefinementCase<dim - 1>(ref_case));
1027  ++subface)
1028  {
1029  project_to_subface(mutation,
1030  face,
1031  subface,
1032  help,
1033  RefinementCase<dim - 1>(ref_case));
1034  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1035  }
1036 
1037  // next copy over weights
1038  for (unsigned int face = 0; face < n_faces; ++face)
1039  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1040  ref_case >= RefinementCase<dim - 1>::cut_x;
1041  --ref_case)
1042  for (unsigned int subface = 0;
1044  RefinementCase<dim - 1>(ref_case));
1045  ++subface)
1046  std::copy(mutation.get_weights().begin(),
1047  mutation.get_weights().end(),
1048  std::back_inserter(weights));
1049  }
1050 
1051  Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1052  ExcInternalError());
1053  Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1054  ExcInternalError());
1055 
1056  return Quadrature<dim>(q_points, weights);
1057 }
1058 
1059 
1060 
1061 // This function is not used in the library
1062 template <int dim>
1065  const unsigned int child_no)
1066 {
1068 
1069  const unsigned int n_q_points = quadrature.size();
1070 
1071  std::vector<Point<dim>> q_points(n_q_points);
1072  for (unsigned int i = 0; i < n_q_points; ++i)
1073  q_points[i] =
1075  child_no);
1076 
1077  // for the weights, things are
1078  // equally simple: copy them and
1079  // scale them
1080  std::vector<double> weights = quadrature.get_weights();
1081  for (unsigned int i = 0; i < n_q_points; ++i)
1082  weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1083 
1084  return Quadrature<dim>(q_points, weights);
1085 }
1086 
1087 
1088 template <int dim>
1091 {
1092  const unsigned int n_points = quadrature.size(),
1094 
1095  std::vector<Point<dim>> q_points(n_points * n_children);
1096  std::vector<double> weights(n_points * n_children);
1097 
1098  // project to each child and copy
1099  // results
1100  for (unsigned int child = 0; child < n_children; ++child)
1101  {
1102  Quadrature<dim> help = project_to_child(quadrature, child);
1103  for (unsigned int i = 0; i < n_points; ++i)
1104  {
1105  q_points[child * n_points + i] = help.point(i);
1106  weights[child * n_points + i] = help.weight(i);
1107  }
1108  }
1109  return Quadrature<dim>(q_points, weights);
1110 }
1111 
1112 
1113 
1114 template <int dim>
1117  const Point<dim> & p1,
1118  const Point<dim> & p2)
1119 {
1120  const unsigned int n = quadrature.size();
1121  std::vector<Point<dim>> points(n);
1122  std::vector<double> weights(n);
1123  const double length = p1.distance(p2);
1124 
1125  for (unsigned int k = 0; k < n; ++k)
1126  {
1127  const double alpha = quadrature.point(k)(0);
1128  points[k] = alpha * p2;
1129  points[k] += (1. - alpha) * p1;
1130  weights[k] = length * quadrature.weight(k);
1131  }
1132  return Quadrature<dim>(points, weights);
1133 }
1134 
1135 
1136 
1137 template <int dim>
1139 QProjector<dim>::DataSetDescriptor::face(const unsigned int face_no,
1140  const bool face_orientation,
1141  const bool face_flip,
1142  const bool face_rotation,
1143  const unsigned int n_quadrature_points)
1144 {
1146 
1147  switch (dim)
1148  {
1149  case 1:
1150  case 2:
1151  return face_no * n_quadrature_points;
1152 
1153 
1154  case 3:
1155  {
1156  // in 3d, we have to account for faces that
1157  // have non-standard face orientation, flip
1158  // and rotation. thus, we have to store
1159  // _eight_ data sets per face or subface
1160 
1161  // set up a table with the according offsets
1162  // for non-standard orientation, first index:
1163  // face_orientation (standard true=1), second
1164  // index: face_flip (standard false=0), third
1165  // index: face_rotation (standard false=0)
1166  //
1167  // note, that normally we should use the
1168  // obvious offsets 0,1,2,3,4,5,6,7. However,
1169  // prior to the changes enabling flipped and
1170  // rotated faces, in many places of the
1171  // library the convention was used, that the
1172  // first dataset with offset 0 corresponds to
1173  // a face in standard orientation. therefore
1174  // we use the offsets 4,5,6,7,0,1,2,3 here to
1175  // stick to that (implicit) convention
1176  static const unsigned int offset[2][2][2] = {
1178  5 * GeometryInfo<dim>::
1179  faces_per_cell}, // face_orientation=false; face_flip=false;
1180  // face_rotation=false and true
1182  7 * GeometryInfo<dim>::
1183  faces_per_cell}}, // face_orientation=false; face_flip=true;
1184  // face_rotation=false and true
1186  1 * GeometryInfo<dim>::
1187  faces_per_cell}, // face_orientation=true; face_flip=false;
1188  // face_rotation=false and true
1190  3 * GeometryInfo<dim>::
1191  faces_per_cell}}}; // face_orientation=true; face_flip=true;
1192  // face_rotation=false and true
1193 
1194  return (
1195  (face_no + offset[face_orientation][face_flip][face_rotation]) *
1196  n_quadrature_points);
1197  }
1198 
1199  default:
1200  Assert(false, ExcInternalError());
1201  }
1203 }
1204 
1205 
1206 
1207 template <>
1210  const unsigned int face_no,
1211  const unsigned int subface_no,
1212  const bool,
1213  const bool,
1214  const bool,
1215  const unsigned int n_quadrature_points,
1217 {
1220  ExcInternalError());
1221 
1222  return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1223  n_quadrature_points);
1224 }
1225 
1226 
1227 
1228 template <>
1231  const unsigned int face_no,
1232  const unsigned int subface_no,
1233  const bool,
1234  const bool,
1235  const bool,
1236  const unsigned int n_quadrature_points,
1238 {
1241  ExcInternalError());
1242 
1243  return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1244  n_quadrature_points);
1245 }
1246 
1247 
1248 template <>
1251  const unsigned int face_no,
1252  const unsigned int subface_no,
1253  const bool face_orientation,
1254  const bool face_flip,
1255  const bool face_rotation,
1256  const unsigned int n_quadrature_points,
1257  const internal::SubfaceCase<3> ref_case)
1258 {
1259  const unsigned int dim = 3;
1260 
1263  ExcInternalError());
1264 
1265  // As the quadrature points created by
1266  // QProjector are on subfaces in their
1267  // "standard location" we have to use a
1268  // permutation of the equivalent subface
1269  // number in order to respect face
1270  // orientation, flip and rotation. The
1271  // information we need here is exactly the
1272  // same as the
1273  // GeometryInfo<3>::child_cell_on_face info
1274  // for the bottom face (face 4) of a hex, as
1275  // on this the RefineCase of the cell matches
1276  // that of the face and the subfaces are
1277  // numbered in the same way as the child
1278  // cells.
1279 
1280  // in 3d, we have to account for faces that
1281  // have non-standard face orientation, flip
1282  // and rotation. thus, we have to store
1283  // _eight_ data sets per face or subface
1284  // already for the isotropic
1285  // case. Additionally, we have three
1286  // different refinement cases, resulting in
1287  // <tt>4 + 2 + 2 = 8</tt> different subfaces
1288  // for each face.
1289  const unsigned int total_subfaces_per_face = 8;
1290 
1291  // set up a table with the according offsets
1292  // for non-standard orientation, first index:
1293  // face_orientation (standard true=1), second
1294  // index: face_flip (standard false=0), third
1295  // index: face_rotation (standard false=0)
1296  //
1297  // note, that normally we should use the
1298  // obvious offsets 0,1,2,3,4,5,6,7. However,
1299  // prior to the changes enabling flipped and
1300  // rotated faces, in many places of the
1301  // library the convention was used, that the
1302  // first dataset with offset 0 corresponds to
1303  // a face in standard orientation. therefore
1304  // we use the offsets 4,5,6,7,0,1,2,3 here to
1305  // stick to that (implicit) convention
1306  static const unsigned int orientation_offset[2][2][2] = {
1307  {// face_orientation=false; face_flip=false; face_rotation=false and true
1308  {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1309  5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1310  // face_orientation=false; face_flip=true; face_rotation=false and true
1311  {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1312  7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
1313  {// face_orientation=true; face_flip=false; face_rotation=false and true
1314  {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1315  1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1316  // face_orientation=true; face_flip=true; face_rotation=false and true
1317  {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1318  3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};
1319 
1320  // set up a table with the offsets for a
1321  // given refinement case respecting the
1322  // corresponding number of subfaces. the
1323  // index corresponds to (RefineCase::Type - 1)
1324 
1325  // note, that normally we should use the
1326  // obvious offsets 0,2,6. However, prior to
1327  // the implementation of anisotropic
1328  // refinement, in many places of the library
1329  // the convention was used, that the first
1330  // dataset with offset 0 corresponds to a
1331  // standard (isotropic) face
1332  // refinement. therefore we use the offsets
1333  // 6,4,0 here to stick to that (implicit)
1334  // convention
1335  static const unsigned int ref_case_offset[3] = {
1336  6, // cut_x
1337  4, // cut_y
1338  0 // cut_xy
1339  };
1340 
1341 
1342  // for each subface of a given FaceRefineCase
1343  // there is a corresponding equivalent
1344  // subface number of one of the "standard"
1345  // RefineCases (cut_x, cut_y, cut_xy). Map
1346  // the given values to those equivalent
1347  // ones.
1348 
1349  // first, define an invalid number
1350  static const unsigned int e = numbers::invalid_unsigned_int;
1351 
1352  static const RefinementCase<dim - 1>
1353  equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
1355  // case_none. there should be only
1356  // invalid values here. However, as
1357  // this function is also called (in
1358  // tests) for cells which have no
1359  // refined faces, use isotropic
1360  // refinement instead
1361  {RefinementCase<dim - 1>::cut_xy,
1362  RefinementCase<dim - 1>::cut_xy,
1363  RefinementCase<dim - 1>::cut_xy,
1364  RefinementCase<dim - 1>::cut_xy},
1365  // case_x
1366  {RefinementCase<dim - 1>::cut_x,
1367  RefinementCase<dim - 1>::cut_x,
1368  RefinementCase<dim - 1>::no_refinement,
1369  RefinementCase<dim - 1>::no_refinement},
1370  // case_x1y
1371  {RefinementCase<dim - 1>::cut_xy,
1372  RefinementCase<dim - 1>::cut_xy,
1373  RefinementCase<dim - 1>::cut_x,
1374  RefinementCase<dim - 1>::no_refinement},
1375  // case_x2y
1376  {RefinementCase<dim - 1>::cut_x,
1377  RefinementCase<dim - 1>::cut_xy,
1378  RefinementCase<dim - 1>::cut_xy,
1379  RefinementCase<dim - 1>::no_refinement},
1380  // case_x1y2y
1381  {RefinementCase<dim - 1>::cut_xy,
1382  RefinementCase<dim - 1>::cut_xy,
1383  RefinementCase<dim - 1>::cut_xy,
1384  RefinementCase<dim - 1>::cut_xy},
1385  // case_y
1386  {RefinementCase<dim - 1>::cut_y,
1387  RefinementCase<dim - 1>::cut_y,
1388  RefinementCase<dim - 1>::no_refinement,
1389  RefinementCase<dim - 1>::no_refinement},
1390  // case_y1x
1391  {RefinementCase<dim - 1>::cut_xy,
1392  RefinementCase<dim - 1>::cut_xy,
1393  RefinementCase<dim - 1>::cut_y,
1394  RefinementCase<dim - 1>::no_refinement},
1395  // case_y2x
1396  {RefinementCase<dim - 1>::cut_y,
1397  RefinementCase<dim - 1>::cut_xy,
1398  RefinementCase<dim - 1>::cut_xy,
1399  RefinementCase<dim - 1>::no_refinement},
1400  // case_y1x2x
1401  {RefinementCase<dim - 1>::cut_xy,
1402  RefinementCase<dim - 1>::cut_xy,
1403  RefinementCase<dim - 1>::cut_xy,
1404  RefinementCase<dim - 1>::cut_xy},
1405  // case_xy (case_isotropic)
1406  {RefinementCase<dim - 1>::cut_xy,
1407  RefinementCase<dim - 1>::cut_xy,
1408  RefinementCase<dim - 1>::cut_xy,
1409  RefinementCase<dim - 1>::cut_xy}};
1410 
1411  static const unsigned int
1412  equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
1414  // case_none, see above
1415  {0, 1, 2, 3},
1416  // case_x
1417  {0, 1, e, e},
1418  // case_x1y
1419  {0, 2, 1, e},
1420  // case_x2y
1421  {0, 1, 3, e},
1422  // case_x1y2y
1423  {0, 2, 1, 3},
1424  // case_y
1425  {0, 1, e, e},
1426  // case_y1x
1427  {0, 1, 1, e},
1428  // case_y2x
1429  {0, 2, 3, e},
1430  // case_y1x2x
1431  {0, 1, 2, 3},
1432  // case_xy (case_isotropic)
1433  {0, 1, 2, 3}};
1434 
1435  // If face-orientation or face_rotation are
1436  // non-standard, cut_x and cut_y have to be
1437  // exchanged.
1438  static const RefinementCase<dim - 1> ref_case_permutation[4] = {
1439  RefinementCase<dim - 1>::no_refinement,
1440  RefinementCase<dim - 1>::cut_y,
1441  RefinementCase<dim - 1>::cut_x,
1442  RefinementCase<dim - 1>::cut_xy};
1443 
1444  // set a corresponding (equivalent)
1445  // RefineCase and subface number
1446  const RefinementCase<dim - 1> equ_ref_case =
1447  equivalent_refine_case[ref_case][subface_no];
1448  const unsigned int equ_subface_no =
1449  equivalent_subface_number[ref_case][subface_no];
1450  // make sure, that we got a valid subface and RefineCase
1452  ExcInternalError());
1453  Assert(equ_subface_no != e, ExcInternalError());
1454  // now, finally respect non-standard faces
1455  const RefinementCase<dim - 1> final_ref_case =
1456  (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1457  equ_ref_case);
1458 
1459  // what we have now is the number of
1460  // the subface in the natural
1461  // orientation of the *face*. what we
1462  // need to know is the number of the
1463  // subface concerning the standard face
1464  // orientation as seen from the *cell*.
1465 
1466  // this mapping is not trivial, but we
1467  // have done exactly this stuff in the
1468  // child_cell_on_face function. in
1469  // order to reduce the amount of code
1470  // as well as to make maintaining the
1471  // functionality easier we want to
1472  // reuse that information. So we note
1473  // that on the bottom face (face 4) of
1474  // a hex cell the local x and y
1475  // coordinates of the face and the cell
1476  // coincide, thus also the refinement
1477  // case of the face corresponds to the
1478  // refinement case of the cell
1479  // (ignoring cell refinement along the
1480  // z direction). Using this knowledge
1481  // we can (ab)use the
1482  // child_cell_on_face function to do
1483  // exactly the transformation we are in
1484  // need of now
1485  const unsigned int final_subface_no =
1487  4,
1488  equ_subface_no,
1489  face_orientation,
1490  face_flip,
1491  face_rotation,
1492  equ_ref_case);
1493 
1494  return (((face_no * total_subfaces_per_face +
1495  ref_case_offset[final_ref_case - 1] + final_subface_no) +
1496  orientation_offset[face_orientation][face_flip][face_rotation]) *
1497  n_quadrature_points);
1498 }
1499 
1500 
1501 
1502 template <int dim>
1505  const unsigned int face_no)
1506 {
1507  std::vector<Point<dim>> points(quadrature.size());
1508  project_to_face(quadrature, face_no, points);
1509  return Quadrature<dim>(points, quadrature.get_weights());
1510 }
1511 
1512 
1513 template <int dim>
1516  const unsigned int face_no,
1517  const unsigned int subface_no,
1518  const RefinementCase<dim - 1> &ref_case)
1519 {
1520  std::vector<Point<dim>> points(quadrature.size());
1521  project_to_subface(quadrature, face_no, subface_no, points, ref_case);
1522  return Quadrature<dim>(points, quadrature.get_weights());
1523 }
1524 
1525 
1526 // ------------------------------------------------------------ //
1527 
1528 namespace internal
1529 {
1530  namespace QIteratedImplementation
1531  {
1532  namespace
1533  {
1534  bool
1535  uses_both_endpoints(const Quadrature<1> &base_quadrature)
1536  {
1537  const bool at_left =
1538  std::any_of(base_quadrature.get_points().cbegin(),
1539  base_quadrature.get_points().cend(),
1540  [](const Point<1> &p) { return p == Point<1>{0.}; });
1541  const bool at_right =
1542  std::any_of(base_quadrature.get_points().cbegin(),
1543  base_quadrature.get_points().cend(),
1544  [](const Point<1> &p) { return p == Point<1>{1.}; });
1545  return (at_left && at_right);
1546  }
1547  } // namespace
1548  } // namespace QIteratedImplementation
1549 } // namespace internal
1550 
1551 // template <>
1552 // void
1553 // QIterated<1>::fill(Quadrature<1>& dst,
1554 // const Quadrature<1> &base_quadrature,
1555 // const unsigned int n_copies)
1556 // {
1557 // Assert (n_copies > 0, ExcZero());
1558 // Assert (base_quadrature.size() > 0, ExcZero());
1559 
1560 // const unsigned int np =
1561 // uses_both_endpoints(base_quadrature)
1562 // ? (base_quadrature.size()-1) * n_copies + 1
1563 // : base_quadrature.size() * n_copies;
1564 
1565 // dst.quadrature_points.resize(np);
1566 // dst.weights.resize(np);
1567 
1568 // if (!uses_both_endpoints(base_quadrature))
1569 // // we don't have to skip some
1570 // // points in order to get a
1571 // // reasonable quadrature formula
1572 // {
1573 // unsigned int next_point = 0;
1574 // for (unsigned int copy=0; copy<n_copies; ++copy)
1575 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
1576 // {
1577 // dst.quadrature_points[next_point](0)
1578 // = (copy + base_quadrature.point(q_point)(0)) / n_copies;
1579 // dst.weights[next_point]
1580 // = base_quadrature.weight(q_point) / n_copies;
1581 // ++next_point;
1582 // }
1583 // }
1584 // else
1585 // // skip doubly available points
1586 // {
1587 // unsigned int next_point = 0;
1588 
1589 // // first find out the weights of
1590 // // the left and the right boundary
1591 // // points. note that these usually
1592 // // are but need not necessarily be
1593 // // the same
1594 // double double_point_weight = 0;
1595 // unsigned int n_end_points = 0;
1596 // for (unsigned int i=0; i<base_quadrature.size(); ++i)
1597 // // add up the weight if this
1598 // // is an endpoint
1599 // if ((base_quadrature.point(i)(0) == 0.) ||
1600 // (base_quadrature.point(i)(0) == 1.))
1601 // {
1602 // double_point_weight += base_quadrature.weight(i);
1603 // ++n_end_points;
1604 // }
1605 // // scale the weight correctly
1606 // double_point_weight /= n_copies;
1607 
1608 // // make sure the base quadrature formula
1609 // // has only one quadrature point
1610 // // per end point
1611 // Assert (n_end_points == 2, ExcInvalidQuadratureFormula());
1612 
1613 
1614 // for (unsigned int copy=0; copy<n_copies; ++copy)
1615 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
1616 // {
1617 // // skip the left point of
1618 // // this copy since we
1619 // // have already entered
1620 // // it the last time
1621 // if ((copy > 0) &&
1622 // (base_quadrature.point(q_point)(0) == 0.))
1623 // continue;
1624 
1625 // dst.quadrature_points[next_point](0)
1626 // = (copy+base_quadrature.point(q_point)(0)) / n_copies;
1627 
1628 // // if this is the
1629 // // rightmost point of one
1630 // // of the non-last
1631 // // copies: give it the
1632 // // double weight
1633 // if ((copy != n_copies-1) &&
1634 // (base_quadrature.point(q_point)(0) == 1.))
1635 // dst.weights[next_point] = double_point_weight;
1636 // else
1637 // dst.weights[next_point] = base_quadrature.weight(q_point) /
1638 // n_copies;
1639 
1640 // ++next_point;
1641 // }
1642 // }
1643 
1644 // #if DEBUG
1645 // double sum_of_weights = 0;
1646 // for (unsigned int i=0; i<dst.size(); ++i)
1647 // sum_of_weights += dst.weight(i);
1648 // Assert (std::fabs(sum_of_weights-1) < 1e-15,
1649 // ExcInternalError());
1650 // #endif
1651 
1652 // }
1653 
1654 
1655 template <>
1656 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
1657  : Quadrature<0>()
1658 {
1659  Assert(false, ExcNotImplemented());
1660 }
1661 
1662 
1663 
1664 template <>
1666  const unsigned int n_copies)
1667  : Quadrature<1>(
1668  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
1669  (base_quadrature.size() - 1) * n_copies + 1 :
1670  base_quadrature.size() * n_copies)
1671 {
1672  // fill(*this, base_quadrature, n_copies);
1673  Assert(base_quadrature.size() > 0, ExcNotInitialized());
1674  Assert(n_copies > 0, ExcZero());
1675 
1676  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
1677  // we don't have to skip some
1678  // points in order to get a
1679  // reasonable quadrature formula
1680  {
1681  unsigned int next_point = 0;
1682  for (unsigned int copy = 0; copy < n_copies; ++copy)
1683  for (unsigned int q_point = 0; q_point < base_quadrature.size();
1684  ++q_point)
1685  {
1686  this->quadrature_points[next_point] =
1687  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
1688  (1.0 * copy) / n_copies);
1689  this->weights[next_point] =
1690  base_quadrature.weight(q_point) / n_copies;
1691 
1692  ++next_point;
1693  }
1694  }
1695  else
1696  // skip doubly available points
1697  {
1698  unsigned int next_point = 0;
1699 
1700  // first find out the weights of
1701  // the left and the right boundary
1702  // points. note that these usually
1703  // are but need not necessarily be
1704  // the same
1705  double double_point_weight = 0;
1706  unsigned int n_end_points = 0;
1707  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
1708  // add up the weight if this
1709  // is an endpoint
1710  if ((base_quadrature.point(i) == Point<1>(0.0)) ||
1711  (base_quadrature.point(i) == Point<1>(1.0)))
1712  {
1713  double_point_weight += base_quadrature.weight(i);
1714  ++n_end_points;
1715  }
1716  // scale the weight correctly
1717  double_point_weight /= n_copies;
1718 
1719  // make sure the base quadrature formula
1720  // has only one quadrature point
1721  // per end point
1722  Assert(n_end_points == 2, ExcInvalidQuadratureFormula());
1723 
1724 
1725  for (unsigned int copy = 0; copy < n_copies; ++copy)
1726  for (unsigned int q_point = 0; q_point < base_quadrature.size();
1727  ++q_point)
1728  {
1729  // skip the left point of
1730  // this copy since we
1731  // have already entered
1732  // it the last time
1733  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
1734  continue;
1735 
1736  this->quadrature_points[next_point] =
1737  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
1738  (1.0 * copy) / n_copies);
1739 
1740  // if this is the
1741  // rightmost point of one
1742  // of the non-last
1743  // copies: give it the
1744  // double weight
1745  if ((copy != n_copies - 1) &&
1746  (base_quadrature.point(q_point) == Point<1>(1.0)))
1747  this->weights[next_point] = double_point_weight;
1748  else
1749  this->weights[next_point] =
1750  base_quadrature.weight(q_point) / n_copies;
1751 
1752  ++next_point;
1753  }
1754  }
1755 
1756 #if DEBUG
1757  double sum_of_weights = 0;
1758  for (unsigned int i = 0; i < this->size(); ++i)
1759  sum_of_weights += this->weight(i);
1760  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
1761 #endif
1762 }
1763 
1764 
1765 // template <int dim>
1766 // void
1767 // QIterated<dim>::fill(Quadrature<dim>&, const Quadrature<1>&, unsigned int)
1768 // {
1769 // Assert(false, ExcNotImplemented());
1770 // }
1771 
1772 
1773 // construct higher dimensional quadrature formula by tensor product
1774 // of lower dimensional iterated quadrature formulae
1775 template <int dim>
1777  const unsigned int N)
1778  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
1779  QIterated<1>(base_quadrature, N))
1780 {}
1781 
1782 
1783 
1784 // explicit instantiations; note: we need them all for all dimensions
1785 template class Quadrature<0>;
1786 template class Quadrature<1>;
1787 template class Quadrature<2>;
1788 template class Quadrature<3>;
1789 template class QAnisotropic<1>;
1790 template class QAnisotropic<2>;
1791 template class QAnisotropic<3>;
1792 template class QIterated<1>;
1793 template class QIterated<2>;
1794 template class QIterated<3>;
1795 template class QProjector<1>;
1796 template class QProjector<2>;
1797 template class QProjector<3>;
1798 
static const unsigned int invalid_unsigned_int
Definition: types.h:191
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: quadrature.cc:1064
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
std::vector< double > weights
Definition: quadrature.h:288
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping, const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:444
const std::vector< Point< dim > > & get_points() const
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim >> &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
const std::vector< double > & get_weights() const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Quadrature(const unsigned int n_quadrature_points=0)
Definition: quadrature.cc:40
static Quadrature< 2 > rotate(const Quadrature< 2 > &q, const unsigned int n_times)
Definition: quadrature.cc:451
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
const Point< dim > & point(const unsigned int i) const
STL namespace.
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:318
static ::ExceptionBase & ExcNotInitialized()
Definition: point.h:110
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:277
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
T fixed_power(const T t)
Definition: utilities.h:1045
static Quadrature< dim > project_to_line(const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
Definition: quadrature.cc:1116
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:346
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1403
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
std::size_t memory_consumption() const
Definition: quadrature.cc:306
static ::ExceptionBase & ExcInvalidQuadratureFormula()
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
bool is_tensor_product_flag
Definition: quadrature.h:297
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:282
void rotate(const double angle, Triangulation< dim > &triangulation)
unsigned int size() const
static Quadrature< 2 > reflect(const Quadrature< 2 > &q)
Definition: quadrature.cc:433
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
Definition: cuda.h:31
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
bool operator==(const Quadrature< dim > &p) const
Definition: quadrature.cc:297
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:1776
static const char N
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:303
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static ::ExceptionBase & ExcNotImplemented()
void initialize(const std::vector< Point< dim >> &points, const std::vector< double > &weights)
Definition: quadrature.cc:50
static ::ExceptionBase & ExcZero()
void copy(const T *begin, const T *end, U *dest)
bool is_tensor_product() const
static Quadrature< dim > project_to_all_children(const Quadrature< dim > &quadrature)
Definition: quadrature.cc:1090
double weight(const unsigned int i) const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: quadrature.cc:1139