Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
21 #include <deal.II/base/utilities.h>
22 
23 #include <cmath>
24 #include <cstdlib>
25 #include <iterator>
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.))
89  , is_tensor_product_flag(true)
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 <>
111 Quadrature<0>::Quadrature(const SubQuadrature &, const Quadrature<1> &)
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())
123  , is_tensor_product_flag(q1.is_tensor_product())
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_cxx14::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 <>
166 Quadrature<1>::Quadrature(const SubQuadrature &, const Quadrature<1> &q2)
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()
226  , quadrature_points(Utilities::fixed_power<dim>(q.size()))
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_cxx14::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)
266  , is_tensor_product_flag(q.is_tensor_product_flag)
267 {
268  if (dim > 1 && is_tensor_product_flag)
269  tensor_basis =
270  std_cxx14::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
271 }
272 
273 
274 
275 template <int dim>
278 {
279  weights = q.weights;
281  is_tensor_product_flag = q.is_tensor_product_flag;
282  if (dim > 1 && is_tensor_product_flag)
283  {
284  if (tensor_basis == nullptr)
285  tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>(
286  *q.tensor_basis);
287  else
288  *tensor_basis = *q.tensor_basis;
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 =
388  std_cxx14::make_unique<std::array<Quadrature<1>, 2>>(q_array);
389 }
390 
391 
392 
393 template <int dim>
395  const Quadrature<1> &qy,
396  const Quadrature<1> &qz)
397  : Quadrature<dim>(qx.size() * qy.size() * qz.size())
398 {
399  Assert(dim == 3, ExcImpossibleInDim(dim));
400 }
401 
402 
403 
404 template <>
406  const Quadrature<1> &qy,
407  const Quadrature<1> &qz)
408  : Quadrature<3>(qx.size() * qy.size() * qz.size())
409 {
410  unsigned int k = 0;
411  for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
412  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
413  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
414  {
415  this->quadrature_points[k](0) = qx.point(k1)(0);
416  this->quadrature_points[k](1) = qy.point(k2)(0);
417  this->quadrature_points[k](2) = qz.point(k3)(0);
418  this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
419  }
420  Assert(k == this->size(), ExcInternalError());
421  this->is_tensor_product_flag = true;
422  const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
423  this->tensor_basis =
424  std_cxx14::make_unique<std::array<Quadrature<1>, 3>>(q_array);
425 }
426 
427 
428 
429 //---------------------------------------------------------------------------
430 
431 
432 
433 template <int dim>
436 {
437  std::vector<Point<2>> q_points(q.size());
438  std::vector<double> weights(q.size());
439  for (unsigned int i = 0; i < q.size(); ++i)
440  {
441  q_points[i][0] = q.point(i)[1];
442  q_points[i][1] = q.point(i)[0];
443 
444  weights[i] = q.weight(i);
445  }
446 
447  return Quadrature<2>(q_points, weights);
448 }
449 
450 
451 template <int dim>
453 QProjector<dim>::rotate(const Quadrature<2> &q, const unsigned int n_times)
454 {
455  std::vector<Point<2>> q_points(q.size());
456  std::vector<double> weights(q.size());
457  for (unsigned int i = 0; i < q.size(); ++i)
458  {
459  switch (n_times % 4)
460  {
461  case 0:
462  // 0 degree
463  q_points[i][0] = q.point(i)[0];
464  q_points[i][1] = q.point(i)[1];
465  break;
466  case 1:
467  // 90 degree counterclockwise
468  q_points[i][0] = 1.0 - q.point(i)[1];
469  q_points[i][1] = q.point(i)[0];
470  break;
471  case 2:
472  // 180 degree counterclockwise
473  q_points[i][0] = 1.0 - q.point(i)[0];
474  q_points[i][1] = 1.0 - q.point(i)[1];
475  break;
476  case 3:
477  // 270 degree counterclockwise
478  q_points[i][0] = q.point(i)[1];
479  q_points[i][1] = 1.0 - q.point(i)[0];
480  break;
481  }
482 
483  weights[i] = q.weight(i);
484  }
485 
486  return Quadrature<2>(q_points, weights);
487 }
488 
489 
490 template <>
491 void
493  const unsigned int face_no,
494  std::vector<Point<1>> &q_points)
495 {
496  const unsigned int dim = 1;
498  AssertDimension(q_points.size(), 1);
499 
500  q_points[0] = Point<dim>(static_cast<double>(face_no));
501 }
502 
503 
504 
505 template <>
506 void
508  const unsigned int face_no,
509  std::vector<Point<2>> &q_points)
510 {
511  const unsigned int dim = 2;
513  Assert(q_points.size() == quadrature.size(),
514  ExcDimensionMismatch(q_points.size(), quadrature.size()));
515 
516  for (unsigned int p = 0; p < quadrature.size(); ++p)
517  switch (face_no)
518  {
519  case 0:
520  q_points[p] = Point<dim>(0, quadrature.point(p)(0));
521  break;
522  case 1:
523  q_points[p] = Point<dim>(1, quadrature.point(p)(0));
524  break;
525  case 2:
526  q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
527  break;
528  case 3:
529  q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
530  break;
531  default:
532  Assert(false, ExcInternalError());
533  }
534 }
535 
536 
537 
538 template <>
539 void
541  const unsigned int face_no,
542  std::vector<Point<3>> &q_points)
543 {
544  const unsigned int dim = 3;
546  Assert(q_points.size() == quadrature.size(),
547  ExcDimensionMismatch(q_points.size(), quadrature.size()));
548 
549  for (unsigned int p = 0; p < quadrature.size(); ++p)
550  switch (face_no)
551  {
552  case 0:
553  q_points[p] =
554  Point<dim>(0, quadrature.point(p)(0), quadrature.point(p)(1));
555  break;
556  case 1:
557  q_points[p] =
558  Point<dim>(1, quadrature.point(p)(0), quadrature.point(p)(1));
559  break;
560  case 2:
561  q_points[p] =
562  Point<dim>(quadrature.point(p)(1), 0, quadrature.point(p)(0));
563  break;
564  case 3:
565  q_points[p] =
566  Point<dim>(quadrature.point(p)(1), 1, quadrature.point(p)(0));
567  break;
568  case 4:
569  q_points[p] =
570  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 0);
571  break;
572  case 5:
573  q_points[p] =
574  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 1);
575  break;
576 
577  default:
578  Assert(false, ExcInternalError());
579  }
580 }
581 
582 
583 
584 template <>
585 void
587  const unsigned int face_no,
588  const unsigned int,
589  std::vector<Point<1>> &q_points,
590  const RefinementCase<0> &)
591 {
592  const unsigned int dim = 1;
594  AssertDimension(q_points.size(), 1);
595 
596  q_points[0] = Point<dim>(static_cast<double>(face_no));
597 }
598 
599 
600 
601 template <>
602 void
604  const unsigned int face_no,
605  const unsigned int subface_no,
606  std::vector<Point<2>> &q_points,
607  const RefinementCase<1> &)
608 {
609  const unsigned int dim = 2;
612 
613  Assert(q_points.size() == quadrature.size(),
614  ExcDimensionMismatch(q_points.size(), quadrature.size()));
615 
616  for (unsigned int p = 0; p < quadrature.size(); ++p)
617  switch (face_no)
618  {
619  case 0:
620  switch (subface_no)
621  {
622  case 0:
623  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
624  break;
625  case 1:
626  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
627  break;
628  default:
629  Assert(false, ExcInternalError());
630  }
631  break;
632  case 1:
633  switch (subface_no)
634  {
635  case 0:
636  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
637  break;
638  case 1:
639  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
640  break;
641  default:
642  Assert(false, ExcInternalError());
643  }
644  break;
645  case 2:
646  switch (subface_no)
647  {
648  case 0:
649  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
650  break;
651  case 1:
652  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
653  break;
654  default:
655  Assert(false, ExcInternalError());
656  }
657  break;
658  case 3:
659  switch (subface_no)
660  {
661  case 0:
662  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
663  break;
664  case 1:
665  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
666  break;
667  default:
668  Assert(false, ExcInternalError());
669  }
670  break;
671 
672  default:
673  Assert(false, ExcInternalError());
674  }
675 }
676 
677 
678 
679 template <>
680 void
682  const unsigned int face_no,
683  const unsigned int subface_no,
684  std::vector<Point<3>> & q_points,
685  const RefinementCase<2> &ref_case)
686 {
687  const unsigned int dim = 3;
690  Assert(q_points.size() == quadrature.size(),
691  ExcDimensionMismatch(q_points.size(), quadrature.size()));
692 
693  // one coordinate is at a const value. for
694  // faces 0, 2 and 4 this value is 0.0, for
695  // faces 1, 3 and 5 it is 1.0
696  double const_value = face_no % 2;
697  // local 2d coordinates are xi and eta,
698  // global 3d coordinates are x, y and
699  // z. those have to be mapped. the following
700  // indices tell, which global coordinate
701  // (0->x, 1->y, 2->z) corresponds to which
702  // local one
703  unsigned int xi_index = numbers::invalid_unsigned_int,
704  eta_index = numbers::invalid_unsigned_int,
705  const_index = face_no / 2;
706  // the xi and eta values have to be scaled
707  // (by factor 0.5 or factor 1.0) depending on
708  // the refinement case and translated (by 0.0
709  // or 0.5) depending on the refinement case
710  // and subface_no.
711  double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
712  eta_translation = 0.0;
713  // set the index mapping between local and
714  // global coordinates
715  switch (face_no / 2)
716  {
717  case 0:
718  xi_index = 1;
719  eta_index = 2;
720  break;
721  case 1:
722  xi_index = 2;
723  eta_index = 0;
724  break;
725  case 2:
726  xi_index = 0;
727  eta_index = 1;
728  break;
729  }
730  // set the scale and translation parameter
731  // for individual subfaces
732  switch (ref_case)
733  {
734  case RefinementCase<dim - 1>::cut_x:
735  xi_scale = 0.5;
736  xi_translation = subface_no % 2 * 0.5;
737  break;
738  case RefinementCase<dim - 1>::cut_y:
739  eta_scale = 0.5;
740  eta_translation = subface_no % 2 * 0.5;
741  break;
742  case RefinementCase<dim - 1>::cut_xy:
743  xi_scale = 0.5;
744  eta_scale = 0.5;
745  xi_translation = int(subface_no % 2) * 0.5;
746  eta_translation = int(subface_no / 2) * 0.5;
747  break;
748  default:
749  Assert(false, ExcInternalError());
750  break;
751  }
752  // finally, compute the scaled, translated,
753  // projected quadrature points
754  for (unsigned int p = 0; p < quadrature.size(); ++p)
755  {
756  q_points[p][xi_index] =
757  xi_scale * quadrature.point(p)(0) + xi_translation;
758  q_points[p][eta_index] =
759  eta_scale * quadrature.point(p)(1) + eta_translation;
760  q_points[p][const_index] = const_value;
761  }
762 }
763 
764 
765 template <>
768 {
769  const unsigned int dim = 1;
770 
771  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;
772 
773  // first fix quadrature points
774  std::vector<Point<dim>> q_points;
775  q_points.reserve(n_points * n_faces);
776  std::vector<Point<dim>> help(n_points);
777 
778 
779  // project to each face and append
780  // results
781  for (unsigned int face = 0; face < n_faces; ++face)
782  {
783  project_to_face(quadrature, face, help);
784  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
785  }
786 
787  // next copy over weights
788  std::vector<double> weights;
789  weights.reserve(n_points * n_faces);
790  for (unsigned int face = 0; face < n_faces; ++face)
791  std::copy(quadrature.get_weights().begin(),
792  quadrature.get_weights().end(),
793  std::back_inserter(weights));
794 
795  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
796  Assert(weights.size() == n_points * n_faces, ExcInternalError());
797 
798  return Quadrature<dim>(q_points, weights);
799 }
800 
801 
802 
803 template <>
806 {
807  const unsigned int dim = 2;
808 
809  const unsigned int n_points = quadrature.size(),
811 
812  // first fix quadrature points
813  std::vector<Point<dim>> q_points;
814  q_points.reserve(n_points * n_faces);
815  std::vector<Point<dim>> help(n_points);
816 
817  // project to each face and append
818  // results
819  for (unsigned int face = 0; face < n_faces; ++face)
820  {
821  project_to_face(quadrature, face, help);
822  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
823  }
824 
825  // next copy over weights
826  std::vector<double> weights;
827  weights.reserve(n_points * n_faces);
828  for (unsigned int face = 0; face < n_faces; ++face)
829  std::copy(quadrature.get_weights().begin(),
830  quadrature.get_weights().end(),
831  std::back_inserter(weights));
832 
833  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
834  Assert(weights.size() == n_points * n_faces, ExcInternalError());
835 
836  return Quadrature<dim>(q_points, weights);
837 }
838 
839 
840 
841 template <>
844 {
845  const unsigned int dim = 3;
846 
847  SubQuadrature q_reflected = reflect(quadrature);
848  SubQuadrature q[8] = {quadrature,
849  rotate(quadrature, 1),
850  rotate(quadrature, 2),
851  rotate(quadrature, 3),
852  q_reflected,
853  rotate(q_reflected, 3),
854  rotate(q_reflected, 2),
855  rotate(q_reflected, 1)};
856 
857 
858 
859  const unsigned int n_points = quadrature.size(),
861 
862  // first fix quadrature points
863  std::vector<Point<dim>> q_points;
864  q_points.reserve(n_points * n_faces * 8);
865  std::vector<Point<dim>> help(n_points);
866 
867  std::vector<double> weights;
868  weights.reserve(n_points * n_faces * 8);
869 
870  // do the following for all possible
871  // mutations of a face (mutation==0
872  // corresponds to a face with standard
873  // orientation, no flip and no rotation)
874  for (const auto &mutation : q)
875  {
876  // project to each face and append
877  // results
878  for (unsigned int face = 0; face < n_faces; ++face)
879  {
880  project_to_face(mutation, face, help);
881  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
882  }
883 
884  // next copy over weights
885  for (unsigned int face = 0; face < n_faces; ++face)
886  std::copy(mutation.get_weights().begin(),
887  mutation.get_weights().end(),
888  std::back_inserter(weights));
889  }
890 
891 
892  Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError());
893  Assert(weights.size() == n_points * n_faces * 8, ExcInternalError());
894 
895  return Quadrature<dim>(q_points, weights);
896 }
897 
898 
899 
900 template <>
903 {
904  const unsigned int dim = 1;
905 
906  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
907  subfaces_per_face =
909 
910  // first fix quadrature points
911  std::vector<Point<dim>> q_points;
912  q_points.reserve(n_points * n_faces * subfaces_per_face);
913  std::vector<Point<dim>> help(n_points);
914 
915  // project to each face and copy
916  // results
917  for (unsigned int face = 0; face < n_faces; ++face)
918  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
919  {
920  project_to_subface(quadrature, face, subface, help);
921  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
922  }
923 
924  // next copy over weights
925  std::vector<double> weights;
926  weights.reserve(n_points * n_faces * subfaces_per_face);
927  for (unsigned int face = 0; face < n_faces; ++face)
928  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
929  std::copy(quadrature.get_weights().begin(),
930  quadrature.get_weights().end(),
931  std::back_inserter(weights));
932 
933  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
934  ExcInternalError());
935  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
936  ExcInternalError());
937 
938  return Quadrature<dim>(q_points, weights);
939 }
940 
941 
942 
943 template <>
946 {
947  const unsigned int dim = 2;
948 
949  const unsigned int n_points = quadrature.size(),
951  subfaces_per_face =
953 
954  // first fix quadrature points
955  std::vector<Point<dim>> q_points;
956  q_points.reserve(n_points * n_faces * subfaces_per_face);
957  std::vector<Point<dim>> help(n_points);
958 
959  // project to each face and copy
960  // results
961  for (unsigned int face = 0; face < n_faces; ++face)
962  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
963  {
964  project_to_subface(quadrature, face, subface, help);
965  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
966  }
967 
968  // next copy over weights
969  std::vector<double> weights;
970  weights.reserve(n_points * n_faces * subfaces_per_face);
971  for (unsigned int face = 0; face < n_faces; ++face)
972  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
973  std::copy(quadrature.get_weights().begin(),
974  quadrature.get_weights().end(),
975  std::back_inserter(weights));
976 
977  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
978  ExcInternalError());
979  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
980  ExcInternalError());
981 
982  return Quadrature<dim>(q_points, weights);
983 }
984 
985 
986 
987 template <>
990 {
991  const unsigned int dim = 3;
992  SubQuadrature q_reflected = reflect(quadrature);
993  SubQuadrature q[8] = {quadrature,
994  rotate(quadrature, 1),
995  rotate(quadrature, 2),
996  rotate(quadrature, 3),
997  q_reflected,
998  rotate(q_reflected, 3),
999  rotate(q_reflected, 2),
1000  rotate(q_reflected, 1)};
1001 
1002  const unsigned int n_points = quadrature.size(),
1004  total_subfaces_per_face = 2 + 2 + 4;
1005 
1006  // first fix quadrature points
1007  std::vector<Point<dim>> q_points;
1008  q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1009  std::vector<Point<dim>> help(n_points);
1010 
1011  std::vector<double> weights;
1012  weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1013 
1014  // do the following for all possible
1015  // mutations of a face (mutation==0
1016  // corresponds to a face with standard
1017  // orientation, no flip and no rotation)
1018  for (const auto &mutation : q)
1019  {
1020  // project to each face and copy
1021  // results
1022  for (unsigned int face = 0; face < n_faces; ++face)
1023  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1024  ref_case >= RefinementCase<dim - 1>::cut_x;
1025  --ref_case)
1026  for (unsigned int subface = 0;
1028  RefinementCase<dim - 1>(ref_case));
1029  ++subface)
1030  {
1031  project_to_subface(mutation,
1032  face,
1033  subface,
1034  help,
1035  RefinementCase<dim - 1>(ref_case));
1036  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1037  }
1038 
1039  // next copy over weights
1040  for (unsigned int face = 0; face < n_faces; ++face)
1041  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1042  ref_case >= RefinementCase<dim - 1>::cut_x;
1043  --ref_case)
1044  for (unsigned int subface = 0;
1046  RefinementCase<dim - 1>(ref_case));
1047  ++subface)
1048  std::copy(mutation.get_weights().begin(),
1049  mutation.get_weights().end(),
1050  std::back_inserter(weights));
1051  }
1052 
1053  Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1054  ExcInternalError());
1055  Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1056  ExcInternalError());
1057 
1058  return Quadrature<dim>(q_points, weights);
1059 }
1060 
1061 
1062 
1063 // This function is not used in the library
1064 template <int dim>
1067  const unsigned int child_no)
1068 {
1070 
1071  const unsigned int n_q_points = quadrature.size();
1072 
1073  std::vector<Point<dim>> q_points(n_q_points);
1074  for (unsigned int i = 0; i < n_q_points; ++i)
1075  q_points[i] =
1077  child_no);
1078 
1079  // for the weights, things are
1080  // equally simple: copy them and
1081  // scale them
1082  std::vector<double> weights = quadrature.get_weights();
1083  for (unsigned int i = 0; i < n_q_points; ++i)
1084  weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);
1085 
1086  return Quadrature<dim>(q_points, weights);
1087 }
1088 
1089 
1090 template <int dim>
1093 {
1094  const unsigned int n_points = quadrature.size(),
1096 
1097  std::vector<Point<dim>> q_points(n_points * n_children);
1098  std::vector<double> weights(n_points * n_children);
1099 
1100  // project to each child and copy
1101  // results
1102  for (unsigned int child = 0; child < n_children; ++child)
1103  {
1104  Quadrature<dim> help = project_to_child(quadrature, child);
1105  for (unsigned int i = 0; i < n_points; ++i)
1106  {
1107  q_points[child * n_points + i] = help.point(i);
1108  weights[child * n_points + i] = help.weight(i);
1109  }
1110  }
1111  return Quadrature<dim>(q_points, weights);
1112 }
1113 
1114 
1115 
1116 template <int dim>
1119  const Point<dim> & p1,
1120  const Point<dim> & p2)
1121 {
1122  const unsigned int n = quadrature.size();
1123  std::vector<Point<dim>> points(n);
1124  std::vector<double> weights(n);
1125  const double length = p1.distance(p2);
1126 
1127  for (unsigned int k = 0; k < n; ++k)
1128  {
1129  const double alpha = quadrature.point(k)(0);
1130  points[k] = alpha * p2;
1131  points[k] += (1. - alpha) * p1;
1132  weights[k] = length * quadrature.weight(k);
1133  }
1134  return Quadrature<dim>(points, weights);
1135 }
1136 
1137 
1138 
1139 template <int dim>
1141 QProjector<dim>::DataSetDescriptor::face(const unsigned int face_no,
1142  const bool face_orientation,
1143  const bool face_flip,
1144  const bool face_rotation,
1145  const unsigned int n_quadrature_points)
1146 {
1148 
1149  switch (dim)
1150  {
1151  case 1:
1152  case 2:
1153  return face_no * n_quadrature_points;
1154 
1155 
1156  case 3:
1157  {
1158  // in 3d, we have to account for faces that
1159  // have non-standard face orientation, flip
1160  // and rotation. thus, we have to store
1161  // _eight_ data sets per face or subface
1162 
1163  // set up a table with the according offsets
1164  // for non-standard orientation, first index:
1165  // face_orientation (standard true=1), second
1166  // index: face_flip (standard false=0), third
1167  // index: face_rotation (standard false=0)
1168  //
1169  // note, that normally we should use the
1170  // obvious offsets 0,1,2,3,4,5,6,7. However,
1171  // prior to the changes enabling flipped and
1172  // rotated faces, in many places of the
1173  // library the convention was used, that the
1174  // first dataset with offset 0 corresponds to
1175  // a face in standard orientation. therefore
1176  // we use the offsets 4,5,6,7,0,1,2,3 here to
1177  // stick to that (implicit) convention
1178  static const unsigned int offset[2][2][2] = {
1180  5 * GeometryInfo<dim>::
1181  faces_per_cell}, // face_orientation=false; face_flip=false;
1182  // face_rotation=false and true
1184  7 * GeometryInfo<dim>::
1185  faces_per_cell}}, // face_orientation=false; face_flip=true;
1186  // face_rotation=false and true
1188  1 * GeometryInfo<dim>::
1189  faces_per_cell}, // face_orientation=true; face_flip=false;
1190  // face_rotation=false and true
1192  3 * GeometryInfo<dim>::
1193  faces_per_cell}}}; // face_orientation=true; face_flip=true;
1194  // face_rotation=false and true
1195 
1196  return (
1197  (face_no + offset[face_orientation][face_flip][face_rotation]) *
1198  n_quadrature_points);
1199  }
1200 
1201  default:
1202  Assert(false, ExcInternalError());
1203  }
1205 }
1206 
1207 
1208 
1209 template <>
1212  const unsigned int face_no,
1213  const unsigned int subface_no,
1214  const bool,
1215  const bool,
1216  const bool,
1217  const unsigned int n_quadrature_points,
1219 {
1222  ExcInternalError());
1223 
1224  return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
1225  n_quadrature_points);
1226 }
1227 
1228 
1229 
1230 template <>
1233  const unsigned int face_no,
1234  const unsigned int subface_no,
1235  const bool,
1236  const bool,
1237  const bool,
1238  const unsigned int n_quadrature_points,
1240 {
1243  ExcInternalError());
1244 
1245  return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
1246  n_quadrature_points);
1247 }
1248 
1249 
1250 template <>
1253  const unsigned int face_no,
1254  const unsigned int subface_no,
1255  const bool face_orientation,
1256  const bool face_flip,
1257  const bool face_rotation,
1258  const unsigned int n_quadrature_points,
1259  const internal::SubfaceCase<3> ref_case)
1260 {
1261  const unsigned int dim = 3;
1262 
1265  ExcInternalError());
1266 
1267  // As the quadrature points created by
1268  // QProjector are on subfaces in their
1269  // "standard location" we have to use a
1270  // permutation of the equivalent subface
1271  // number in order to respect face
1272  // orientation, flip and rotation. The
1273  // information we need here is exactly the
1274  // same as the
1275  // GeometryInfo<3>::child_cell_on_face info
1276  // for the bottom face (face 4) of a hex, as
1277  // on this the RefineCase of the cell matches
1278  // that of the face and the subfaces are
1279  // numbered in the same way as the child
1280  // cells.
1281 
1282  // in 3d, we have to account for faces that
1283  // have non-standard face orientation, flip
1284  // and rotation. thus, we have to store
1285  // _eight_ data sets per face or subface
1286  // already for the isotropic
1287  // case. Additionally, we have three
1288  // different refinement cases, resulting in
1289  // <tt>4 + 2 + 2 = 8</tt> different subfaces
1290  // for each face.
1291  const unsigned int total_subfaces_per_face = 8;
1292 
1293  // set up a table with the according offsets
1294  // for non-standard orientation, first index:
1295  // face_orientation (standard true=1), second
1296  // index: face_flip (standard false=0), third
1297  // index: face_rotation (standard false=0)
1298  //
1299  // note, that normally we should use the
1300  // obvious offsets 0,1,2,3,4,5,6,7. However,
1301  // prior to the changes enabling flipped and
1302  // rotated faces, in many places of the
1303  // library the convention was used, that the
1304  // first dataset with offset 0 corresponds to
1305  // a face in standard orientation. therefore
1306  // we use the offsets 4,5,6,7,0,1,2,3 here to
1307  // stick to that (implicit) convention
1308  static const unsigned int orientation_offset[2][2][2] = {
1309  {// face_orientation=false; face_flip=false; face_rotation=false and true
1310  {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1311  5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1312  // face_orientation=false; face_flip=true; face_rotation=false and true
1313  {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1314  7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
1315  {// face_orientation=true; face_flip=false; face_rotation=false and true
1316  {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1317  1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
1318  // face_orientation=true; face_flip=true; face_rotation=false and true
1319  {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
1320  3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};
1321 
1322  // set up a table with the offsets for a
1323  // given refinement case respecting the
1324  // corresponding number of subfaces. the
1325  // index corresponds to (RefineCase::Type - 1)
1326 
1327  // note, that normally we should use the
1328  // obvious offsets 0,2,6. However, prior to
1329  // the implementation of anisotropic
1330  // refinement, in many places of the library
1331  // the convention was used, that the first
1332  // dataset with offset 0 corresponds to a
1333  // standard (isotropic) face
1334  // refinement. therefore we use the offsets
1335  // 6,4,0 here to stick to that (implicit)
1336  // convention
1337  static const unsigned int ref_case_offset[3] = {
1338  6, // cut_x
1339  4, // cut_y
1340  0 // cut_xy
1341  };
1342 
1343 
1344  // for each subface of a given FaceRefineCase
1345  // there is a corresponding equivalent
1346  // subface number of one of the "standard"
1347  // RefineCases (cut_x, cut_y, cut_xy). Map
1348  // the given values to those equivalent
1349  // ones.
1350 
1351  // first, define an invalid number
1352  static const unsigned int e = numbers::invalid_unsigned_int;
1353 
1354  static const RefinementCase<dim - 1>
1355  equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
1357  // case_none. there should be only
1358  // invalid values here. However, as
1359  // this function is also called (in
1360  // tests) for cells which have no
1361  // refined faces, use isotropic
1362  // refinement instead
1363  {RefinementCase<dim - 1>::cut_xy,
1364  RefinementCase<dim - 1>::cut_xy,
1365  RefinementCase<dim - 1>::cut_xy,
1366  RefinementCase<dim - 1>::cut_xy},
1367  // case_x
1368  {RefinementCase<dim - 1>::cut_x,
1369  RefinementCase<dim - 1>::cut_x,
1370  RefinementCase<dim - 1>::no_refinement,
1371  RefinementCase<dim - 1>::no_refinement},
1372  // case_x1y
1373  {RefinementCase<dim - 1>::cut_xy,
1374  RefinementCase<dim - 1>::cut_xy,
1375  RefinementCase<dim - 1>::cut_x,
1376  RefinementCase<dim - 1>::no_refinement},
1377  // case_x2y
1378  {RefinementCase<dim - 1>::cut_x,
1379  RefinementCase<dim - 1>::cut_xy,
1380  RefinementCase<dim - 1>::cut_xy,
1381  RefinementCase<dim - 1>::no_refinement},
1382  // case_x1y2y
1383  {RefinementCase<dim - 1>::cut_xy,
1384  RefinementCase<dim - 1>::cut_xy,
1385  RefinementCase<dim - 1>::cut_xy,
1386  RefinementCase<dim - 1>::cut_xy},
1387  // case_y
1388  {RefinementCase<dim - 1>::cut_y,
1389  RefinementCase<dim - 1>::cut_y,
1390  RefinementCase<dim - 1>::no_refinement,
1391  RefinementCase<dim - 1>::no_refinement},
1392  // case_y1x
1393  {RefinementCase<dim - 1>::cut_xy,
1394  RefinementCase<dim - 1>::cut_xy,
1395  RefinementCase<dim - 1>::cut_y,
1396  RefinementCase<dim - 1>::no_refinement},
1397  // case_y2x
1398  {RefinementCase<dim - 1>::cut_y,
1399  RefinementCase<dim - 1>::cut_xy,
1400  RefinementCase<dim - 1>::cut_xy,
1401  RefinementCase<dim - 1>::no_refinement},
1402  // case_y1x2x
1403  {RefinementCase<dim - 1>::cut_xy,
1404  RefinementCase<dim - 1>::cut_xy,
1405  RefinementCase<dim - 1>::cut_xy,
1406  RefinementCase<dim - 1>::cut_xy},
1407  // case_xy (case_isotropic)
1408  {RefinementCase<dim - 1>::cut_xy,
1409  RefinementCase<dim - 1>::cut_xy,
1410  RefinementCase<dim - 1>::cut_xy,
1411  RefinementCase<dim - 1>::cut_xy}};
1412 
1413  static const unsigned int
1414  equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
1416  // case_none, see above
1417  {0, 1, 2, 3},
1418  // case_x
1419  {0, 1, e, e},
1420  // case_x1y
1421  {0, 2, 1, e},
1422  // case_x2y
1423  {0, 1, 3, e},
1424  // case_x1y2y
1425  {0, 2, 1, 3},
1426  // case_y
1427  {0, 1, e, e},
1428  // case_y1x
1429  {0, 1, 1, e},
1430  // case_y2x
1431  {0, 2, 3, e},
1432  // case_y1x2x
1433  {0, 1, 2, 3},
1434  // case_xy (case_isotropic)
1435  {0, 1, 2, 3}};
1436 
1437  // If face-orientation or face_rotation are
1438  // non-standard, cut_x and cut_y have to be
1439  // exchanged.
1440  static const RefinementCase<dim - 1> ref_case_permutation[4] = {
1441  RefinementCase<dim - 1>::no_refinement,
1442  RefinementCase<dim - 1>::cut_y,
1443  RefinementCase<dim - 1>::cut_x,
1444  RefinementCase<dim - 1>::cut_xy};
1445 
1446  // set a corresponding (equivalent)
1447  // RefineCase and subface number
1448  const RefinementCase<dim - 1> equ_ref_case =
1449  equivalent_refine_case[ref_case][subface_no];
1450  const unsigned int equ_subface_no =
1451  equivalent_subface_number[ref_case][subface_no];
1452  // make sure, that we got a valid subface and RefineCase
1454  ExcInternalError());
1455  Assert(equ_subface_no != e, ExcInternalError());
1456  // now, finally respect non-standard faces
1457  const RefinementCase<dim - 1> final_ref_case =
1458  (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
1459  equ_ref_case);
1460 
1461  // what we have now is the number of
1462  // the subface in the natural
1463  // orientation of the *face*. what we
1464  // need to know is the number of the
1465  // subface concerning the standard face
1466  // orientation as seen from the *cell*.
1467 
1468  // this mapping is not trivial, but we
1469  // have done exactly this stuff in the
1470  // child_cell_on_face function. in
1471  // order to reduce the amount of code
1472  // as well as to make maintaining the
1473  // functionality easier we want to
1474  // reuse that information. So we note
1475  // that on the bottom face (face 4) of
1476  // a hex cell the local x and y
1477  // coordinates of the face and the cell
1478  // coincide, thus also the refinement
1479  // case of the face corresponds to the
1480  // refinement case of the cell
1481  // (ignoring cell refinement along the
1482  // z direction). Using this knowledge
1483  // we can (ab)use the
1484  // child_cell_on_face function to do
1485  // exactly the transformation we are in
1486  // need of now
1487  const unsigned int final_subface_no =
1489  4,
1490  equ_subface_no,
1491  face_orientation,
1492  face_flip,
1493  face_rotation,
1494  equ_ref_case);
1495 
1496  return (((face_no * total_subfaces_per_face +
1497  ref_case_offset[final_ref_case - 1] + final_subface_no) +
1498  orientation_offset[face_orientation][face_flip][face_rotation]) *
1499  n_quadrature_points);
1500 }
1501 
1502 
1503 
1504 template <int dim>
1507  const unsigned int face_no)
1508 {
1509  std::vector<Point<dim>> points(quadrature.size());
1510  project_to_face(quadrature, face_no, points);
1511  return Quadrature<dim>(points, quadrature.get_weights());
1512 }
1513 
1514 
1515 template <int dim>
1518  const unsigned int face_no,
1519  const unsigned int subface_no,
1520  const RefinementCase<dim - 1> &ref_case)
1521 {
1522  std::vector<Point<dim>> points(quadrature.size());
1523  project_to_subface(quadrature, face_no, subface_no, points, ref_case);
1524  return Quadrature<dim>(points, quadrature.get_weights());
1525 }
1526 
1527 
1528 // ------------------------------------------------------------ //
1529 
1530 namespace internal
1531 {
1532  namespace QIteratedImplementation
1533  {
1534  namespace
1535  {
1536  bool
1537  uses_both_endpoints(const Quadrature<1> &base_quadrature)
1538  {
1539  const bool at_left =
1540  std::any_of(base_quadrature.get_points().cbegin(),
1541  base_quadrature.get_points().cend(),
1542  [](const Point<1> &p) { return p == Point<1>{0.}; });
1543  const bool at_right =
1544  std::any_of(base_quadrature.get_points().cbegin(),
1545  base_quadrature.get_points().cend(),
1546  [](const Point<1> &p) { return p == Point<1>{1.}; });
1547  return (at_left && at_right);
1548  }
1549  } // namespace
1550  } // namespace QIteratedImplementation
1551 } // namespace internal
1552 
1553 // template <>
1554 // void
1555 // QIterated<1>::fill(Quadrature<1>& dst,
1556 // const Quadrature<1> &base_quadrature,
1557 // const unsigned int n_copies)
1558 // {
1559 // Assert (n_copies > 0, ExcZero());
1560 // Assert (base_quadrature.size() > 0, ExcZero());
1561 
1562 // const unsigned int np =
1563 // uses_both_endpoints(base_quadrature)
1564 // ? (base_quadrature.size()-1) * n_copies + 1
1565 // : base_quadrature.size() * n_copies;
1566 
1567 // dst.quadrature_points.resize(np);
1568 // dst.weights.resize(np);
1569 
1570 // if (!uses_both_endpoints(base_quadrature))
1571 // // we don't have to skip some
1572 // // points in order to get a
1573 // // reasonable quadrature formula
1574 // {
1575 // unsigned int next_point = 0;
1576 // for (unsigned int copy=0; copy<n_copies; ++copy)
1577 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
1578 // {
1579 // dst.quadrature_points[next_point](0)
1580 // = (copy + base_quadrature.point(q_point)(0)) / n_copies;
1581 // dst.weights[next_point]
1582 // = base_quadrature.weight(q_point) / n_copies;
1583 // ++next_point;
1584 // }
1585 // }
1586 // else
1587 // // skip doubly available points
1588 // {
1589 // unsigned int next_point = 0;
1590 
1591 // // first find out the weights of
1592 // // the left and the right boundary
1593 // // points. note that these usually
1594 // // are but need not necessarily be
1595 // // the same
1596 // double double_point_weight = 0;
1597 // unsigned int n_end_points = 0;
1598 // for (unsigned int i=0; i<base_quadrature.size(); ++i)
1599 // // add up the weight if this
1600 // // is an endpoint
1601 // if ((base_quadrature.point(i)(0) == 0.) ||
1602 // (base_quadrature.point(i)(0) == 1.))
1603 // {
1604 // double_point_weight += base_quadrature.weight(i);
1605 // ++n_end_points;
1606 // }
1607 // // scale the weight correctly
1608 // double_point_weight /= n_copies;
1609 
1610 // // make sure the base quadrature formula
1611 // // has only one quadrature point
1612 // // per end point
1613 // Assert (n_end_points == 2, ExcInvalidQuadratureFormula());
1614 
1615 
1616 // for (unsigned int copy=0; copy<n_copies; ++copy)
1617 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
1618 // {
1619 // // skip the left point of
1620 // // this copy since we
1621 // // have already entered
1622 // // it the last time
1623 // if ((copy > 0) &&
1624 // (base_quadrature.point(q_point)(0) == 0.))
1625 // continue;
1626 
1627 // dst.quadrature_points[next_point](0)
1628 // = (copy+base_quadrature.point(q_point)(0)) / n_copies;
1629 
1630 // // if this is the
1631 // // rightmost point of one
1632 // // of the non-last
1633 // // copies: give it the
1634 // // double weight
1635 // if ((copy != n_copies-1) &&
1636 // (base_quadrature.point(q_point)(0) == 1.))
1637 // dst.weights[next_point] = double_point_weight;
1638 // else
1639 // dst.weights[next_point] = base_quadrature.weight(q_point) /
1640 // n_copies;
1641 
1642 // ++next_point;
1643 // }
1644 // }
1645 
1646 // #if DEBUG
1647 // double sum_of_weights = 0;
1648 // for (unsigned int i=0; i<dst.size(); ++i)
1649 // sum_of_weights += dst.weight(i);
1650 // Assert (std::fabs(sum_of_weights-1) < 1e-15,
1651 // ExcInternalError());
1652 // #endif
1653 
1654 // }
1655 
1656 
1657 template <>
1658 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
1659  : Quadrature<0>()
1660 {
1661  Assert(false, ExcNotImplemented());
1662 }
1663 
1664 
1665 
1666 template <>
1668  const unsigned int n_copies)
1669  : Quadrature<1>(
1670  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
1671  (base_quadrature.size() - 1) * n_copies + 1 :
1672  base_quadrature.size() * n_copies)
1673 {
1674  // fill(*this, base_quadrature, n_copies);
1675  Assert(base_quadrature.size() > 0, ExcNotInitialized());
1676  Assert(n_copies > 0, ExcZero());
1677 
1678  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
1679  // we don't have to skip some
1680  // points in order to get a
1681  // reasonable quadrature formula
1682  {
1683  unsigned int next_point = 0;
1684  for (unsigned int copy = 0; copy < n_copies; ++copy)
1685  for (unsigned int q_point = 0; q_point < base_quadrature.size();
1686  ++q_point)
1687  {
1688  this->quadrature_points[next_point] =
1689  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
1690  (1.0 * copy) / n_copies);
1691  this->weights[next_point] =
1692  base_quadrature.weight(q_point) / n_copies;
1693 
1694  ++next_point;
1695  }
1696  }
1697  else
1698  // skip doubly available points
1699  {
1700  unsigned int next_point = 0;
1701 
1702  // first find out the weights of
1703  // the left and the right boundary
1704  // points. note that these usually
1705  // are but need not necessarily be
1706  // the same
1707  double double_point_weight = 0;
1708  unsigned int n_end_points = 0;
1709  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
1710  // add up the weight if this
1711  // is an endpoint
1712  if ((base_quadrature.point(i) == Point<1>(0.0)) ||
1713  (base_quadrature.point(i) == Point<1>(1.0)))
1714  {
1715  double_point_weight += base_quadrature.weight(i);
1716  ++n_end_points;
1717  }
1718  // scale the weight correctly
1719  double_point_weight /= n_copies;
1720 
1721  // make sure the base quadrature formula
1722  // has only one quadrature point
1723  // per end point
1724  Assert(n_end_points == 2, ExcInvalidQuadratureFormula());
1725 
1726 
1727  for (unsigned int copy = 0; copy < n_copies; ++copy)
1728  for (unsigned int q_point = 0; q_point < base_quadrature.size();
1729  ++q_point)
1730  {
1731  // skip the left point of
1732  // this copy since we
1733  // have already entered
1734  // it the last time
1735  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
1736  continue;
1737 
1738  this->quadrature_points[next_point] =
1739  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
1740  (1.0 * copy) / n_copies);
1741 
1742  // if this is the
1743  // rightmost point of one
1744  // of the non-last
1745  // copies: give it the
1746  // double weight
1747  if ((copy != n_copies - 1) &&
1748  (base_quadrature.point(q_point) == Point<1>(1.0)))
1749  this->weights[next_point] = double_point_weight;
1750  else
1751  this->weights[next_point] =
1752  base_quadrature.weight(q_point) / n_copies;
1753 
1754  ++next_point;
1755  }
1756  }
1757 
1758 #if DEBUG
1759  double sum_of_weights = 0;
1760  for (unsigned int i = 0; i < this->size(); ++i)
1761  sum_of_weights += this->weight(i);
1762  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
1763 #endif
1764 }
1765 
1766 
1767 // template <int dim>
1768 // void
1769 // QIterated<dim>::fill(Quadrature<dim>&, const Quadrature<1>&, unsigned int)
1770 // {
1771 // Assert(false, ExcNotImplemented());
1772 // }
1773 
1774 
1775 // construct higher dimensional quadrature formula by tensor product
1776 // of lower dimensional iterated quadrature formulae
1777 template <int dim>
1779  const unsigned int N)
1780  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
1781  QIterated<1>(base_quadrature, N))
1782 {}
1783 
1784 
1785 
1786 // explicit instantiations; note: we need them all for all dimensions
1787 template class Quadrature<0>;
1788 template class Quadrature<1>;
1789 template class Quadrature<2>;
1790 template class Quadrature<3>;
1791 template class QAnisotropic<1>;
1792 template class QAnisotropic<2>;
1793 template class QAnisotropic<3>;
1794 template class QIterated<1>;
1795 template class QIterated<2>;
1796 template class QIterated<3>;
1797 template class QProjector<1>;
1798 template class QProjector<2>;
1799 template class QProjector<3>;
1800 
QAnisotropic
Definition: quadrature.h:320
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
Quadrature::get_tensor_basis
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:318
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Quadrature::memory_consumption
std::size_t memory_consumption() const
Definition: quadrature.cc:306
memory_consumption.h
Quadrature::quadrature_points
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:284
Quadrature::tensor_basis
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:305
utilities.h
QProjector::reflect
static Quadrature< 2 > reflect(const Quadrature< 2 > &q)
Definition: quadrature.cc:435
GeometryInfo
Definition: geometry_info.h:1224
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
Quadrature::weight
double weight(const unsigned int i) const
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Quadrature::operator==
bool operator==(const Quadrature< dim > &p) const
Definition: quadrature.cc:297
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
QProjector::project_to_all_children
static Quadrature< dim > project_to_all_children(const Quadrature< dim > &quadrature)
Definition: quadrature.cc:1092
QProjector::project_to_all_faces
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
QIterated
Definition: quadrature.h:369
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
QProjector::project_to_face
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition: symengine_math.cc:273
Subscriptor
Definition: subscriptor.h:62
QProjector
Definition: qprojector.h:78
Quadrature::point
const Point< dim > & point(const unsigned int i) const
Quadrature::get_weights
const std::vector< double > & get_weights() const
Physics::Elasticity::Kinematics::w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
geometry_info.h
QProjector::project_to_subface
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim >> &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
internal::SubfaceCase
Definition: geometry_info.h:1172
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1072
QProjector::DataSetDescriptor
Definition: qprojector.h:220
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
QProjector::project_to_child
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: quadrature.cc:1066
QProjector::DataSetDescriptor::subface
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
RefinementCase
Definition: geometry_info.h:795
QIterated::QIterated
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:1778
GeometryInfo::child_cell_on_face
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:268
QProjector::rotate
static Quadrature< 2 > rotate(const Quadrature< 2 > &q, const unsigned int n_times)
Definition: quadrature.cc:453
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
Quadrature::weights
std::vector< double > weights
Definition: quadrature.h:290
Quadrature::Quadrature
Quadrature(const unsigned int n_quadrature_points=0)
Definition: quadrature.cc:40
Quadrature::is_tensor_product_flag
bool is_tensor_product_flag
Definition: quadrature.h:299
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
QProjector::DataSetDescriptor::face
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: quadrature.cc:1141
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Particles::Generators::quadrature_points
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)
Definition: generators.cc:442
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
Quadrature::operator=
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:277
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
StandardExceptions::ExcZero
static ::ExceptionBase & ExcZero()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point
Definition: point.h:111
QProjector::project_to_all_subfaces
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
quadrature.h
QAnisotropic::QAnisotropic
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:346
QIterated::ExcInvalidQuadratureFormula
static ::ExceptionBase & ExcInvalidQuadratureFormula()
Quadrature::size
unsigned int size() const
internal
Definition: aligned_vector.h:369
GridTools::rotate
void rotate(const double angle, Triangulation< dim > &triangulation)
LAPACKSupport::N
static const char N
Definition: lapack_support.h:159
memory.h
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Quadrature
Definition: quadrature.h:85
QProjector::project_to_line
static Quadrature< dim > project_to_line(const Quadrature< 1 > &quadrature, const Point< dim > &p1, const Point< dim > &p2)
Definition: quadrature.cc:1118
qprojector.h
Utilities
Definition: cuda.h:31
Quadrature::initialize
void initialize(const std::vector< Point< dim >> &points, const std::vector< double > &weights)
Definition: quadrature.cc:50
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
int