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