Reference documentation for deal.II version Git 0491297983 2019-09-23 09:31:59 +0200
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
quadrature.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/geometry_info.h>
17 #include <deal.II/base/memory_consumption.h>
18 #include <deal.II/base/qprojector.h>
19 #include <deal.II/base/quadrature.h>
20 #include <deal.II/base/std_cxx14/memory.h>
21 #include <deal.II/base/utilities.h>
22 
23 #include <cmath>
24 #include <cstdlib>
25 #include <iterator>
26 
27 DEAL_II_NAMESPACE_OPEN
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 
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 
116 
117 
118 template <int dim>
120  : quadrature_points(q1.size() * q2.size())
121  , weights(q1.size() * q2.size())
123 {
124  unsigned int present_index = 0;
125  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
126  for (unsigned int i1 = 0; i1 < q1.size(); ++i1)
127  {
128  // compose coordinates of
129  // new quadrature point by tensor
130  // product in the last component
131  for (unsigned int d = 0; d < dim - 1; ++d)
132  quadrature_points[present_index](d) = q1.point(i1)(d);
133  quadrature_points[present_index](dim - 1) = q2.point(i2)(0);
134 
135  weights[present_index] = q1.weight(i1) * q2.weight(i2);
136 
137  ++present_index;
138  }
139 
140 #ifdef DEBUG
141  if (size() > 0)
142  {
143  double sum = 0;
144  for (unsigned int i = 0; i < size(); ++i)
145  sum += weights[i];
146  // we cannot guarantee the sum of weights
147  // to be exactly one, but it should be
148  // near that.
149  Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
150  }
151 #endif
152 
154  {
155  tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>();
156  for (unsigned int i = 0; i < dim - 1; ++i)
157  (*tensor_basis)[i] = q1.get_tensor_basis()[i];
158  (*tensor_basis)[dim - 1] = q2;
159  }
160 }
161 
162 
163 
164 template <>
166  : quadrature_points(q2.size())
167  , weights(q2.size())
168  , is_tensor_product_flag(true)
169 {
170  unsigned int present_index = 0;
171  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
172  {
173  // compose coordinates of
174  // new quadrature point by tensor
175  // product in the last component
176  quadrature_points[present_index](0) = q2.point(i2)(0);
177 
178  weights[present_index] = q2.weight(i2);
179 
180  ++present_index;
181  }
182 
183 #ifdef DEBUG
184  if (size() > 0)
185  {
186  double sum = 0;
187  for (unsigned int i = 0; i < size(); ++i)
188  sum += weights[i];
189  // we cannot guarantee the sum of weights
190  // to be exactly one, but it should be
191  // near that.
192  Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
193  }
194 #endif
195 }
196 
197 
198 
199 template <>
201  : Subscriptor()
202  ,
203  // quadrature_points(1),
204  weights(1, 1.)
205  , is_tensor_product_flag(false)
206 {}
207 
208 
209 template <>
211  : Subscriptor()
212 {
213  // this function should never be
214  // called -- this should be the
215  // copy constructor in 1d...
216  Assert(false, ExcImpossibleInDim(1));
217 }
218 
219 
220 
221 template <int dim>
223  : Subscriptor()
224  , quadrature_points(Utilities::fixed_power<dim>(q.size()))
225  , weights(Utilities::fixed_power<dim>(q.size()))
226  , is_tensor_product_flag(true)
227 {
228  Assert(dim <= 3, ExcNotImplemented());
229 
230  const unsigned int n0 = q.size();
231  const unsigned int n1 = (dim > 1) ? n0 : 1;
232  const unsigned int n2 = (dim > 2) ? n0 : 1;
233 
234  unsigned int k = 0;
235  for (unsigned int i2 = 0; i2 < n2; ++i2)
236  for (unsigned int i1 = 0; i1 < n1; ++i1)
237  for (unsigned int i0 = 0; i0 < n0; ++i0)
238  {
239  quadrature_points[k](0) = q.point(i0)(0);
240  if (dim > 1)
241  quadrature_points[k](1) = q.point(i1)(0);
242  if (dim > 2)
243  quadrature_points[k](2) = q.point(i2)(0);
244  weights[k] = q.weight(i0);
245  if (dim > 1)
246  weights[k] *= q.weight(i1);
247  if (dim > 2)
248  weights[k] *= q.weight(i2);
249  ++k;
250  }
251 
252  tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>();
253  for (unsigned int i = 0; i < dim; ++i)
254  (*tensor_basis)[i] = q;
255 }
256 
257 
258 
259 template <int dim>
261  : Subscriptor()
263  , weights(q.weights)
265 {
266  if (dim > 1 && is_tensor_product_flag)
267  tensor_basis =
268  std_cxx14::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
269 }
270 
271 
272 
273 template <int dim>
276 {
277  weights = q.weights;
280  if (dim > 1 && is_tensor_product_flag)
281  {
282  if (tensor_basis == nullptr)
283  tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>(
284  *q.tensor_basis);
285  else
287  }
288  return *this;
289 }
290 
291 
292 
293 template <int dim>
294 bool
296 {
297  return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
298 }
299 
300 
301 
302 template <int dim>
303 std::size_t
305 {
308 }
309 
310 
311 
312 template <int dim>
313 typename std::conditional<dim == 1,
314  std::array<Quadrature<1>, dim>,
315  const std::array<Quadrature<1>, dim> &>::type
317 {
318  Assert(this->is_tensor_product_flag == true,
319  ExcMessage("This function only makes sense if "
320  "this object represents a tensor product!"));
321  Assert(tensor_basis != nullptr, ExcInternalError());
322 
323  return *tensor_basis;
324 }
325 
326 
327 
328 template <>
329 std::array<Quadrature<1>, 1>
331 {
332  Assert(this->is_tensor_product_flag == true,
333  ExcMessage("This function only makes sense if "
334  "this object represents a tensor product!"));
335 
336  return std::array<Quadrature<1>, 1>{{*this}};
337 }
338 
339 
340 
341 //---------------------------------------------------------------------------
342 template <int dim>
344  : Quadrature<dim>(qx.size())
345 {
346  Assert(dim == 1, ExcImpossibleInDim(dim));
347  unsigned int k = 0;
348  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
349  {
350  this->quadrature_points[k](0) = qx.point(k1)(0);
351  this->weights[k++] = qx.weight(k1);
352  }
353  Assert(k == this->size(), ExcInternalError());
354  this->is_tensor_product_flag = true;
355 }
356 
357 
358 
359 template <int dim>
361  const Quadrature<1> &qy)
362  : Quadrature<dim>(qx.size() * qy.size())
363 {
364  Assert(dim == 2, ExcImpossibleInDim(dim));
365 }
366 
367 
368 
369 template <>
371  : Quadrature<2>(qx.size() * qy.size())
372 {
373  unsigned int k = 0;
374  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
375  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
376  {
377  this->quadrature_points[k](0) = qx.point(k1)(0);
378  this->quadrature_points[k](1) = qy.point(k2)(0);
379  this->weights[k++] = qx.weight(k1) * qy.weight(k2);
380  }
381  Assert(k == this->size(), ExcInternalError());
382  this->is_tensor_product_flag = true;
383  const std::array<Quadrature<1>, 2> q_array{{qx, qy}};
384  this->tensor_basis =
385  std_cxx14::make_unique<std::array<Quadrature<1>, 2>>(q_array);
386 }
387 
388 
389 
390 template <int dim>
392  const Quadrature<1> &qy,
393  const Quadrature<1> &qz)
394  : Quadrature<dim>(qx.size() * qy.size() * qz.size())
395 {
396  Assert(dim == 3, ExcImpossibleInDim(dim));
397 }
398 
399 
400 
401 template <>
403  const Quadrature<1> &qy,
404  const Quadrature<1> &qz)
405  : Quadrature<3>(qx.size() * qy.size() * qz.size())
406 {
407  unsigned int k = 0;
408  for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
409  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
410  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
411  {
412  this->quadrature_points[k](0) = qx.point(k1)(0);
413  this->quadrature_points[k](1) = qy.point(k2)(0);
414  this->quadrature_points[k](2) = qz.point(k3)(0);
415  this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
416  }
417  Assert(k == this->size(), ExcInternalError());
418  this->is_tensor_product_flag = true;
419  const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
420  this->tensor_basis =
421  std_cxx14::make_unique<std::array<Quadrature<1>, 3>>(q_array);
422 }
423 
424 
425 
426 //---------------------------------------------------------------------------
427 
428 
429 
430 template <int dim>
433 {
434  std::vector<Point<2>> q_points(q.size());
435  std::vector<double> weights(q.size());
436  for (unsigned int i = 0; i < q.size(); ++i)
437  {
438  q_points[i][0] = q.point(i)[1];
439  q_points[i][1] = q.point(i)[0];
440 
441  weights[i] = q.weight(i);
442  }
443 
444  return Quadrature<2>(q_points, weights);
445 }
446 
447 
448 template <int dim>
450 QProjector<dim>::rotate(const Quadrature<2> &q, const unsigned int n_times)
451 {
452  std::vector<Point<2>> q_points(q.size());
453  std::vector<double> weights(q.size());
454  for (unsigned int i = 0; i < q.size(); ++i)
455  {
456  switch (n_times % 4)
457  {
458  case 0:
459  // 0 degree
460  q_points[i][0] = q.point(i)[0];
461  q_points[i][1] = q.point(i)[1];
462  break;
463  case 1:
464  // 90 degree counterclockwise
465  q_points[i][0] = 1.0 - q.point(i)[1];
466  q_points[i][1] = q.point(i)[0];
467  break;
468  case 2:
469  // 180 degree counterclockwise
470  q_points[i][0] = 1.0 - q.point(i)[0];
471  q_points[i][1] = 1.0 - q.point(i)[1];
472  break;
473  case 3:
474  // 270 degree counterclockwise
475  q_points[i][0] = q.point(i)[1];
476  q_points[i][1] = 1.0 - q.point(i)[0];
477  break;
478  }
479 
480  weights[i] = q.weight(i);
481  }
482 
483  return Quadrature<2>(q_points, weights);
484 }
485 
486 
487 template <>
488 void
490  const unsigned int face_no,
491  std::vector<Point<1>> &q_points)
492 {
493  const unsigned int dim = 1;
495  AssertDimension(q_points.size(), 1);
496 
497  q_points[0] = Point<dim>(static_cast<double>(face_no));
498 }
499 
500 
501 
502 template <>
503 void
505  const unsigned int face_no,
506  std::vector<Point<2>> &q_points)
507 {
508  const unsigned int dim = 2;
510  Assert(q_points.size() == quadrature.size(),
511  ExcDimensionMismatch(q_points.size(), quadrature.size()));
512 
513  for (unsigned int p = 0; p < quadrature.size(); ++p)
514  switch (face_no)
515  {
516  case 0:
517  q_points[p] = Point<dim>(0, quadrature.point(p)(0));
518  break;
519  case 1:
520  q_points[p] = Point<dim>(1, quadrature.point(p)(0));
521  break;
522  case 2:
523  q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
524  break;
525  case 3:
526  q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
527  break;
528  default:
529  Assert(false, ExcInternalError());
530  }
531 }
532 
533 
534 
535 template <>
536 void
538  const unsigned int face_no,
539  std::vector<Point<3>> &q_points)
540 {
541  const unsigned int dim = 3;
543  Assert(q_points.size() == quadrature.size(),
544  ExcDimensionMismatch(q_points.size(), quadrature.size()));
545 
546  for (unsigned int p = 0; p < quadrature.size(); ++p)
547  switch (face_no)
548  {
549  case 0:
550  q_points[p] =
551  Point<dim>(0, quadrature.point(p)(0), quadrature.point(p)(1));
552  break;
553  case 1:
554  q_points[p] =
555  Point<dim>(1, quadrature.point(p)(0), quadrature.point(p)(1));
556  break;
557  case 2:
558  q_points[p] =
559  Point<dim>(quadrature.point(p)(1), 0, quadrature.point(p)(0));
560  break;
561  case 3:
562  q_points[p] =
563  Point<dim>(quadrature.point(p)(1), 1, quadrature.point(p)(0));
564  break;
565  case 4:
566  q_points[p] =
567  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 0);
568  break;
569  case 5:
570  q_points[p] =
571  Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 1);
572  break;
573 
574  default:
575  Assert(false, ExcInternalError());
576  }
577 }
578 
579 
580 
581 template <>
582 void
584  const unsigned int face_no,
585  const unsigned int,
586  std::vector<Point<1>> &q_points,
587  const RefinementCase<0> &)
588 {
589  const unsigned int dim = 1;
591  AssertDimension(q_points.size(), 1);
592 
593  q_points[0] = Point<dim>(static_cast<double>(face_no));
594 }
595 
596 
597 
598 template <>
599 void
601  const unsigned int face_no,
602  const unsigned int subface_no,
603  std::vector<Point<2>> &q_points,
604  const RefinementCase<1> &)
605 {
606  const unsigned int dim = 2;
609 
610  Assert(q_points.size() == quadrature.size(),
611  ExcDimensionMismatch(q_points.size(), quadrature.size()));
612 
613  for (unsigned int p = 0; p < quadrature.size(); ++p)
614  switch (face_no)
615  {
616  case 0:
617  switch (subface_no)
618  {
619  case 0:
620  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
621  break;
622  case 1:
623  q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
624  break;
625  default:
626  Assert(false, ExcInternalError());
627  }
628  break;
629  case 1:
630  switch (subface_no)
631  {
632  case 0:
633  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
634  break;
635  case 1:
636  q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
637  break;
638  default:
639  Assert(false, ExcInternalError());
640  }
641  break;
642  case 2:
643  switch (subface_no)
644  {
645  case 0:
646  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
647  break;
648  case 1:
649  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
650  break;
651  default:
652  Assert(false, ExcInternalError());
653  }
654  break;
655  case 3:
656  switch (subface_no)
657  {
658  case 0:
659  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
660  break;
661  case 1:
662  q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
663  break;
664  default:
665  Assert(false, ExcInternalError());
666  }
667  break;
668 
669  default:
670  Assert(false, ExcInternalError());
671  }
672 }
673 
674 
675 
676 template <>
677 void
679  const unsigned int face_no,
680  const unsigned int subface_no,
681  std::vector<Point<3>> & q_points,
682  const RefinementCase<2> &ref_case)
683 {
684  const unsigned int dim = 3;
687  Assert(q_points.size() == quadrature.size(),
688  ExcDimensionMismatch(q_points.size(), quadrature.size()));
689 
690  // one coordinate is at a const value. for
691  // faces 0, 2 and 4 this value is 0.0, for
692  // faces 1, 3 and 5 it is 1.0
693  double const_value = face_no % 2;
694  // local 2d coordinates are xi and eta,
695  // global 3d coordinates are x, y and
696  // z. those have to be mapped. the following
697  // indices tell, which global coordinate
698  // (0->x, 1->y, 2->z) corresponds to which
699  // local one
700  unsigned int xi_index = numbers::invalid_unsigned_int,
701  eta_index = numbers::invalid_unsigned_int,
702  const_index = face_no / 2;
703  // the xi and eta values have to be scaled
704  // (by factor 0.5 or factor 1.0) depending on
705  // the refinement case and translated (by 0.0
706  // or 0.5) depending on the refinement case
707  // and subface_no.
708  double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
709  eta_translation = 0.0;
710  // set the index mapping between local and
711  // global coordinates
712  switch (face_no / 2)
713  {
714  case 0:
715  xi_index = 1;
716  eta_index = 2;
717  break;
718  case 1:
719  xi_index = 2;
720  eta_index = 0;
721  break;
722  case 2:
723  xi_index = 0;
724  eta_index = 1;
725  break;
726  }
727  // set the scale and translation parameter
728  // for individual subfaces
729  switch (ref_case)
730  {
731  case RefinementCase<dim - 1>::cut_x:
732  xi_scale = 0.5;
733  xi_translation = subface_no % 2 * 0.5;
734  break;
735  case RefinementCase<dim - 1>::cut_y:
736  eta_scale = 0.5;
737  eta_translation = subface_no % 2 * 0.5;
738  break;
739  case RefinementCase<dim - 1>::cut_xy:
740  xi_scale = 0.5;
741  eta_scale = 0.5;
742  xi_translation = int(subface_no % 2) * 0.5;
743  eta_translation = int(subface_no / 2) * 0.5;
744  break;
745  default:
746  Assert(false, ExcInternalError());
747  break;
748  }
749  // finally, compute the scaled, translated,
750  // projected quadrature points
751  for (unsigned int p = 0; p < quadrature.size(); ++p)
752  {
753  q_points[p][xi_index] =
754  xi_scale * quadrature.point(p)(0) + xi_translation;
755  q_points[p][eta_index] =
756  eta_scale * quadrature.point(p)(1) + eta_translation;
757  q_points[p][const_index] = const_value;
758  }
759 }
760 
761 
762 template <>
765 {
766  const unsigned int dim = 1;
767 
768  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;
769 
770  // first fix quadrature points
771  std::vector<Point<dim>> q_points;
772  q_points.reserve(n_points * n_faces);
773  std::vector<Point<dim>> help(n_points);
774 
775 
776  // project to each face and append
777  // results
778  for (unsigned int face = 0; face < n_faces; ++face)
779  {
780  project_to_face(quadrature, face, help);
781  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
782  }
783 
784  // next copy over weights
785  std::vector<double> weights;
786  weights.reserve(n_points * n_faces);
787  for (unsigned int face = 0; face < n_faces; ++face)
788  std::copy(quadrature.get_weights().begin(),
789  quadrature.get_weights().end(),
790  std::back_inserter(weights));
791 
792  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
793  Assert(weights.size() == n_points * n_faces, ExcInternalError());
794 
795  return Quadrature<dim>(q_points, weights);
796 }
797 
798 
799 
800 template <>
803 {
804  const unsigned int dim = 2;
805 
806  const unsigned int n_points = quadrature.size(),
808 
809  // first fix quadrature points
810  std::vector<Point<dim>> q_points;
811  q_points.reserve(n_points * n_faces);
812  std::vector<Point<dim>> help(n_points);
813 
814  // project to each face and append
815  // results
816  for (unsigned int face = 0; face < n_faces; ++face)
817  {
818  project_to_face(quadrature, face, help);
819  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
820  }
821 
822  // next copy over weights
823  std::vector<double> weights;
824  weights.reserve(n_points * n_faces);
825  for (unsigned int face = 0; face < n_faces; ++face)
826  std::copy(quadrature.get_weights().begin(),
827  quadrature.get_weights().end(),
828  std::back_inserter(weights));
829 
830  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
831  Assert(weights.size() == n_points * n_faces, ExcInternalError());
832 
833  return Quadrature<dim>(q_points, weights);
834 }
835 
836 
837 
838 template <>
841 {
842  const unsigned int dim = 3;
843 
844  SubQuadrature q_reflected = reflect(quadrature);
845  SubQuadrature q[8] = {quadrature,
846  rotate(quadrature, 1),
847  rotate(quadrature, 2),
848  rotate(quadrature, 3),
849  q_reflected,
850  rotate(q_reflected, 3),
851  rotate(q_reflected, 2),
852  rotate(q_reflected, 1)};
853 
854 
855 
856  const unsigned int n_points = quadrature.size(),
858 
859  // first fix quadrature points
860  std::vector<Point<dim>> q_points;
861  q_points.reserve(n_points * n_faces * 8);
862  std::vector<Point<dim>> help(n_points);
863 
864  std::vector<double> weights;
865  weights.reserve(n_points * n_faces * 8);
866 
867  // do the following for all possible
868  // mutations of a face (mutation==0
869  // corresponds to a face with standard
870  // orientation, no flip and no rotation)
871  for (const auto &mutation : q)
872  {
873  // project to each face and append
874  // results
875  for (unsigned int face = 0; face < n_faces; ++face)
876  {
877  project_to_face(mutation, face, help);
878  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
879  }
880 
881  // next copy over weights
882  for (unsigned int face = 0; face < n_faces; ++face)
883  std::copy(mutation.get_weights().begin(),
884  mutation.get_weights().end(),
885  std::back_inserter(weights));
886  }
887 
888 
889  Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError());
890  Assert(weights.size() == n_points * n_faces * 8, ExcInternalError());
891 
892  return Quadrature<dim>(q_points, weights);
893 }
894 
895 
896 
897 template <>
900 {
901  const unsigned int dim = 1;
902 
903  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
904  subfaces_per_face =
906 
907  // first fix quadrature points
908  std::vector<Point<dim>> q_points;
909  q_points.reserve(n_points * n_faces * subfaces_per_face);
910  std::vector<Point<dim>> help(n_points);
911 
912  // project to each face and copy
913  // results
914  for (unsigned int face = 0; face < n_faces; ++face)
915  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
916  {
917  project_to_subface(quadrature, face, subface, help);
918  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
919  }
920 
921  // next copy over weights
922  std::vector<double> weights;
923  weights.reserve(n_points * n_faces * subfaces_per_face);
924  for (unsigned int face = 0; face < n_faces; ++face)
925  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
926  std::copy(quadrature.get_weights().begin(),
927  quadrature.get_weights().end(),
928  std::back_inserter(weights));
929 
930  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
931  ExcInternalError());
932  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
933  ExcInternalError());
934 
935  return Quadrature<dim>(q_points, weights);
936 }
937 
938 
939 
940 template <>
943 {
944  const unsigned int dim = 2;
945 
946  const unsigned int n_points = quadrature.size(),
948  subfaces_per_face =
950 
951  // first fix quadrature points
952  std::vector<Point<dim>> q_points;
953  q_points.reserve(n_points * n_faces * subfaces_per_face);
954  std::vector<Point<dim>> help(n_points);
955 
956  // project to each face and copy
957  // results
958  for (unsigned int face = 0; face < n_faces; ++face)
959  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
960  {
961  project_to_subface(quadrature, face, subface, help);
962  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
963  }
964 
965  // next copy over weights
966  std::vector<double> weights;
967  weights.reserve(n_points * n_faces * subfaces_per_face);
968  for (unsigned int face = 0; face < n_faces; ++face)
969  for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
970  std::copy(quadrature.get_weights().begin(),
971  quadrature.get_weights().end(),
972  std::back_inserter(weights));
973 
974  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
975  ExcInternalError());
976  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
977  ExcInternalError());
978 
979  return Quadrature<dim>(q_points, weights);
980 }
981 
982 
983 
984 template <>
987 {
988  const unsigned int dim = 3;
989  SubQuadrature q_reflected = reflect(quadrature);
990  SubQuadrature q[8] = {quadrature,
991  rotate(quadrature, 1),
992  rotate(quadrature, 2),
993  rotate(quadrature, 3),
994  q_reflected,
995  rotate(q_reflected, 3),
996  rotate(q_reflected, 2),
997  rotate(q_reflected, 1)};
998 
999  const unsigned int n_points = quadrature.size(),
1001  total_subfaces_per_face = 2 + 2 + 4;
1002 
1003  // first fix quadrature points
1004  std::vector<Point<dim>> q_points;
1005  q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1006  std::vector<Point<dim>> help(n_points);
1007 
1008  std::vector<double> weights;
1009  weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);
1010 
1011  // do the following for all possible
1012  // mutations of a face (mutation==0
1013  // corresponds to a face with standard
1014  // orientation, no flip and no rotation)
1015  for (const auto &mutation : q)
1016  {
1017  // project to each face and copy
1018  // results
1019  for (unsigned int face = 0; face < n_faces; ++face)
1020  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1021  ref_case >= RefinementCase<dim - 1>::cut_x;
1022  --ref_case)
1023  for (unsigned int subface = 0;
1025  RefinementCase<dim - 1>(ref_case));
1026  ++subface)
1027  {
1028  project_to_subface(mutation,
1029  face,
1030  subface,
1031  help,
1032  RefinementCase<dim - 1>(ref_case));
1033  std::copy(help.begin(), help.end(), std::back_inserter(q_points));
1034  }
1035 
1036  // next copy over weights
1037  for (unsigned int face = 0; face < n_faces; ++face)
1038  for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
1039  ref_case >= RefinementCase<dim - 1>::cut_x;
1040  --ref_case)
1041  for (unsigned int subface = 0;
1043  RefinementCase<dim - 1>(ref_case));
1044  ++subface)
1045  std::copy(mutation.get_weights().begin(),
1046  mutation.get_weights().end(),
1047  std::back_inserter(weights));
1048  }
1049 
1050  Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
1051  ExcInternalError());
1052  Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1053  ExcInternalError());
1054 
1055  return Quadrature<dim>(q_points, weights);
1056 }
1057 
1058 
1059 
1060 // This function is not used in the library
1061 template <int dim>
1064  const unsigned int child_no)
1065 {
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  bool at_left = false, at_right = false;
1538  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
1539  {
1540  if (base_quadrature.point(i) == Point<1>(0.0))
1541  at_left = true;
1542  if (base_quadrature.point(i) == Point<1>(1.0))
1543  at_right = true;
1544  }
1545 
1546  return (at_left && at_right);
1547  }
1548  } // namespace
1549  } // namespace QIteratedImplementation
1550 } // namespace internal
1551 
1552 // template <>
1553 // void
1554 // QIterated<1>::fill(Quadrature<1>& dst,
1555 // const Quadrature<1> &base_quadrature,
1556 // const unsigned int n_copies)
1557 // {
1558 // Assert (n_copies > 0, ExcZero());
1559 // Assert (base_quadrature.size() > 0, ExcZero());
1560 
1561 // const unsigned int np =
1562 // uses_both_endpoints(base_quadrature)
1563 // ? (base_quadrature.size()-1) * n_copies + 1
1564 // : base_quadrature.size() * n_copies;
1565 
1566 // dst.quadrature_points.resize(np);
1567 // dst.weights.resize(np);
1568 
1569 // if (!uses_both_endpoints(base_quadrature))
1570 // // we don't have to skip some
1571 // // points in order to get a
1572 // // reasonable quadrature formula
1573 // {
1574 // unsigned int next_point = 0;
1575 // for (unsigned int copy=0; copy<n_copies; ++copy)
1576 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
1577 // {
1578 // dst.quadrature_points[next_point](0)
1579 // = (copy + base_quadrature.point(q_point)(0)) / n_copies;
1580 // dst.weights[next_point]
1581 // = base_quadrature.weight(q_point) / n_copies;
1582 // ++next_point;
1583 // }
1584 // }
1585 // else
1586 // // skip doubly available points
1587 // {
1588 // unsigned int next_point = 0;
1589 
1590 // // first find out the weights of
1591 // // the left and the right boundary
1592 // // points. note that these usually
1593 // // are but need not necessarily be
1594 // // the same
1595 // double double_point_weight = 0;
1596 // unsigned int n_end_points = 0;
1597 // for (unsigned int i=0; i<base_quadrature.size(); ++i)
1598 // // add up the weight if this
1599 // // is an endpoint
1600 // if ((base_quadrature.point(i)(0) == 0.) ||
1601 // (base_quadrature.point(i)(0) == 1.))
1602 // {
1603 // double_point_weight += base_quadrature.weight(i);
1604 // ++n_end_points;
1605 // }
1606 // // scale the weight correctly
1607 // double_point_weight /= n_copies;
1608 
1609 // // make sure the base quadrature formula
1610 // // has only one quadrature point
1611 // // per end point
1612 // Assert (n_end_points == 2, ExcInvalidQuadratureFormula());
1613 
1614 
1615 // for (unsigned int copy=0; copy<n_copies; ++copy)
1616 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
1617 // {
1618 // // skip the left point of
1619 // // this copy since we
1620 // // have already entered
1621 // // it the last time
1622 // if ((copy > 0) &&
1623 // (base_quadrature.point(q_point)(0) == 0.))
1624 // continue;
1625 
1626 // dst.quadrature_points[next_point](0)
1627 // = (copy+base_quadrature.point(q_point)(0)) / n_copies;
1628 
1629 // // if this is the
1630 // // rightmost point of one
1631 // // of the non-last
1632 // // copies: give it the
1633 // // double weight
1634 // if ((copy != n_copies-1) &&
1635 // (base_quadrature.point(q_point)(0) == 1.))
1636 // dst.weights[next_point] = double_point_weight;
1637 // else
1638 // dst.weights[next_point] = base_quadrature.weight(q_point) /
1639 // n_copies;
1640 
1641 // ++next_point;
1642 // }
1643 // }
1644 
1645 // #if DEBUG
1646 // double sum_of_weights = 0;
1647 // for (unsigned int i=0; i<dst.size(); ++i)
1648 // sum_of_weights += dst.weight(i);
1649 // Assert (std::fabs(sum_of_weights-1) < 1e-15,
1650 // ExcInternalError());
1651 // #endif
1652 
1653 // }
1654 
1655 
1656 template <>
1657 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
1658  : Quadrature<0>()
1659 {
1660  Assert(false, ExcNotImplemented());
1661 }
1662 
1663 
1664 
1665 template <>
1666 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
1667  const unsigned int n_copies)
1668  : Quadrature<1>(
1669  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
1670  (base_quadrature.size() - 1) * n_copies + 1 :
1671  base_quadrature.size() * n_copies)
1672 {
1673  // fill(*this, base_quadrature, n_copies);
1674  Assert(base_quadrature.size() > 0, ExcNotInitialized());
1675  Assert(n_copies > 0, ExcZero());
1676 
1677  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
1678  // we don't have to skip some
1679  // points in order to get a
1680  // reasonable quadrature formula
1681  {
1682  unsigned int next_point = 0;
1683  for (unsigned int copy = 0; copy < n_copies; ++copy)
1684  for (unsigned int q_point = 0; q_point < base_quadrature.size();
1685  ++q_point)
1686  {
1687  this->quadrature_points[next_point] =
1688  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
1689  (1.0 * copy) / n_copies);
1690  this->weights[next_point] =
1691  base_quadrature.weight(q_point) / n_copies;
1692 
1693  ++next_point;
1694  }
1695  }
1696  else
1697  // skip doubly available points
1698  {
1699  unsigned int next_point = 0;
1700 
1701  // first find out the weights of
1702  // the left and the right boundary
1703  // points. note that these usually
1704  // are but need not necessarily be
1705  // the same
1706  double double_point_weight = 0;
1707  unsigned int n_end_points = 0;
1708  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
1709  // add up the weight if this
1710  // is an endpoint
1711  if ((base_quadrature.point(i) == Point<1>(0.0)) ||
1712  (base_quadrature.point(i) == Point<1>(1.0)))
1713  {
1714  double_point_weight += base_quadrature.weight(i);
1715  ++n_end_points;
1716  }
1717  // scale the weight correctly
1718  double_point_weight /= n_copies;
1719 
1720  // make sure the base quadrature formula
1721  // has only one quadrature point
1722  // per end point
1723  Assert(n_end_points == 2, ExcInvalidQuadratureFormula());
1724 
1725 
1726  for (unsigned int copy = 0; copy < n_copies; ++copy)
1727  for (unsigned int q_point = 0; q_point < base_quadrature.size();
1728  ++q_point)
1729  {
1730  // skip the left point of
1731  // this copy since we
1732  // have already entered
1733  // it the last time
1734  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
1735  continue;
1736 
1737  this->quadrature_points[next_point] =
1738  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
1739  (1.0 * copy) / n_copies);
1740 
1741  // if this is the
1742  // rightmost point of one
1743  // of the non-last
1744  // copies: give it the
1745  // double weight
1746  if ((copy != n_copies - 1) &&
1747  (base_quadrature.point(q_point) == Point<1>(1.0)))
1748  this->weights[next_point] = double_point_weight;
1749  else
1750  this->weights[next_point] =
1751  base_quadrature.weight(q_point) / n_copies;
1752 
1753  ++next_point;
1754  }
1755  }
1756 
1757 #if DEBUG
1758  double sum_of_weights = 0;
1759  for (unsigned int i = 0; i < this->size(); ++i)
1760  sum_of_weights += this->weight(i);
1761  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
1762 #endif
1763 }
1764 
1765 
1766 // template <int dim>
1767 // void
1768 // QIterated<dim>::fill(Quadrature<dim>&, const Quadrature<1>&, unsigned int)
1769 // {
1770 // Assert(false, ExcNotImplemented());
1771 // }
1772 
1773 
1774 // construct higher dimensional quadrature formula by tensor product
1775 // of lower dimensional iterated quadrature formulae
1776 template <int dim>
1778  const unsigned int N)
1779  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
1780  QIterated<1>(base_quadrature, N))
1781 {}
1782 
1783 
1784 
1785 // explicit instantiations; note: we need them all for all dimensions
1786 template class Quadrature<0>;
1787 template class Quadrature<1>;
1788 template class Quadrature<2>;
1789 template class Quadrature<3>;
1790 template class QAnisotropic<1>;
1791 template class QAnisotropic<2>;
1792 template class QAnisotropic<3>;
1793 template class QIterated<1>;
1794 template class QIterated<2>;
1795 template class QIterated<3>;
1796 template class QProjector<1>;
1797 template class QProjector<2>;
1798 template class QProjector<3>;
1799 
1800 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
Definition: types.h:187
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: quadrature.cc:1063
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::vector< double > weights
Definition: quadrature.h:290
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
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:450
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
const Point< dim > & point(const unsigned int i) const
STL namespace.
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:316
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Definition: point.h:110
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:275
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
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:343
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::size_t memory_consumption() const
Definition: quadrature.cc:304
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
bool is_tensor_product_flag
Definition: quadrature.h:299
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:284
unsigned int size() const
static Quadrature< 2 > reflect(const Quadrature< 2 > &q)
Definition: quadrature.cc:432
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
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:295
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)
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:1777
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:305
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()
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