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> 26 DEAL_II_NAMESPACE_OPEN
32 quadrature_points (n_q),
47 quadrature_points (n_q,
Point<dim>()),
56 const std::vector<double> &w)
67 const std::vector<double> &
weights)
72 Assert (weights.size() == points.size(),
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)
121 for (
unsigned int d=0; d<dim-1; ++d)
136 for (
unsigned int i=0; i<
size(); ++i)
155 unsigned int present_index = 0;
156 for (
unsigned int i2=0; i2<q2.
size(); ++i2)
173 for (
unsigned int i=0; i<
size(); ++i)
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;
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)
294 for (
unsigned int k1=0; k1<qx.
size(); ++k1)
312 for (
unsigned int k2=0; k2<qy.
size(); ++k2)
313 for (
unsigned int k1=0; k1<qx.
size(); ++k1)
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)
356 std::vector<Point<2> > q_points (q.
size());
358 for (
unsigned int i=0; i<q.
size(); ++i)
360 q_points[i][0] = q.
point(i)[1];
361 q_points[i][1] = q.
point(i)[0];
373 const unsigned int n_times)
375 std::vector<Point<2> > q_points (q.
size());
377 for (
unsigned int i=0; i<q.
size(); ++i)
383 q_points[i][0] = q.
point(i)[0];
384 q_points[i][1] = q.
point(i)[1];
388 q_points[i][0] = 1.0 - q.
point(i)[1];
389 q_points[i][1] = q.
point(i)[0];
393 q_points[i][0] = 1.0 - q.
point(i)[0];
394 q_points[i][1] = 1.0 - q.
point(i)[1];
398 q_points[i][0] = q.
point(i)[1];
399 q_points[i][1] = 1.0 - q.
point(i)[0];
413 const unsigned int face_no,
416 const unsigned int dim=1;
428 const unsigned int face_no,
431 const unsigned int dim=2;
436 for (
unsigned int p=0; p<quadrature.
size(); ++p)
461 const unsigned int face_no,
464 const unsigned int dim=3;
469 for (
unsigned int p=0; p<quadrature.
size(); ++p)
474 quadrature.
point(p)(0),
475 quadrature.
point(p)(1));
479 quadrature.
point(p)(0),
480 quadrature.
point(p)(1));
485 quadrature.
point(p)(0));
490 quadrature.
point(p)(0));
494 quadrature.
point(p)(1),
499 quadrature.
point(p)(1),
513 const unsigned int face_no,
518 const unsigned int dim=1;
530 const unsigned int face_no,
531 const unsigned int subface_no,
535 const unsigned int dim=2;
542 for (
unsigned int p=0; p<quadrature.
size(); ++p)
610 const unsigned int face_no,
611 const unsigned int subface_no,
615 const unsigned int dim=3;
624 double const_value=face_no%2;
633 const_index = face_no/2;
662 switch ((
unsigned char)ref_case)
666 xi_translation = subface_no%2 * 0.5;
670 eta_translation = subface_no%2 * 0.5;
675 xi_translation = int(subface_no%2) * 0.5;
676 eta_translation = int(subface_no/2) * 0.5;
684 for (
unsigned int p=0; p<quadrature.
size(); ++p)
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;
697 const unsigned int dim = 1;
699 const unsigned int n_points = 1,
703 std::vector<Point<dim> > q_points;
704 q_points.reserve(n_points * n_faces);
705 std::vector <Point<dim> > help(n_points);
710 for (
unsigned int face=0; face<n_faces; ++face)
712 project_to_face(quadrature, face, help);
713 std::copy (help.begin(), help.end(),
714 std::back_inserter (q_points));
719 weights.reserve (n_points * n_faces);
720 for (
unsigned int face=0; face<n_faces; ++face)
723 std::back_inserter (weights));
725 Assert (q_points.size() == n_points * n_faces,
727 Assert (weights.size() == n_points * n_faces,
739 const unsigned int dim = 2;
741 const unsigned int n_points = quadrature.
size(),
745 std::vector<Point<dim> > q_points;
746 q_points.reserve(n_points * n_faces);
747 std::vector <Point<dim> > help(n_points);
751 for (
unsigned int face=0; face<n_faces; ++face)
753 project_to_face(quadrature, face, help);
754 std::copy (help.begin(), help.end(),
755 std::back_inserter (q_points));
760 weights.reserve (n_points * n_faces);
761 for (
unsigned int face=0; face<n_faces; ++face)
764 std::back_inserter (weights));
766 Assert (q_points.size() == n_points * n_faces,
768 Assert (weights.size() == n_points * n_faces,
780 const unsigned int dim = 3;
786 rotate (quadrature,1),
787 rotate (quadrature,2),
788 rotate (quadrature,3),
790 rotate (q_reflected,3),
791 rotate (q_reflected,2),
792 rotate (q_reflected,1)
797 const unsigned int n_points = quadrature.
size(),
801 std::vector<Point<dim> > q_points;
802 q_points.reserve(n_points * n_faces * 8);
803 std::vector <Point<dim> > help(n_points);
806 weights.reserve (n_points * n_faces * 8);
812 for (
unsigned int mutation=0; mutation<8; ++mutation)
816 for (
unsigned int face=0; face<n_faces; ++face)
818 project_to_face(q[mutation], face, help);
819 std::copy (help.begin(), help.end(),
820 std::back_inserter (q_points));
824 for (
unsigned int face=0; face<n_faces; ++face)
827 std::back_inserter (weights));
831 Assert (q_points.size() == n_points * n_faces * 8,
833 Assert (weights.size() == n_points * n_faces * 8,
845 const unsigned int dim = 1;
847 const unsigned int n_points = 1,
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);
858 for (
unsigned int face=0; face<n_faces; ++face)
859 for (
unsigned int subface=0; subface<subfaces_per_face; ++subface)
861 project_to_subface(quadrature, face, subface, help);
862 std::copy (help.begin(), help.end(),
863 std::back_inserter (q_points));
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)
873 std::back_inserter (weights));
875 Assert (q_points.size() == n_points * n_faces * subfaces_per_face,
877 Assert (weights.size() == n_points * n_faces * subfaces_per_face,
889 const unsigned int dim = 2;
891 const unsigned int n_points = quadrature.
size(),
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);
902 for (
unsigned int face=0; face<n_faces; ++face)
903 for (
unsigned int subface=0; subface<subfaces_per_face; ++subface)
905 project_to_subface(quadrature, face, subface, help);
906 std::copy (help.begin(), help.end(),
907 std::back_inserter (q_points));
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)
917 std::back_inserter (weights));
919 Assert (q_points.size() == n_points * n_faces * subfaces_per_face,
921 Assert (weights.size() == n_points * n_faces * subfaces_per_face,
933 const unsigned int dim = 3;
938 rotate (quadrature,1),
939 rotate (quadrature,2),
940 rotate (quadrature,3),
942 rotate (q_reflected,3),
943 rotate (q_reflected,2),
944 rotate (q_reflected,1)
947 const unsigned int n_points = quadrature.
size(),
949 total_subfaces_per_face = 2 + 2 + 4;
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);
957 weights.reserve (n_points * n_faces * total_subfaces_per_face * 8);
963 for (
unsigned int mutation=0; mutation<8; ++mutation)
967 for (
unsigned int face=0; face<n_faces; ++face)
971 for (
unsigned int subface=0; subface<GeometryInfo<dim-1>::n_children(
RefinementCase<dim-1>(ref_case)); ++subface)
973 project_to_subface(q[mutation], face, subface, help,
975 std::copy (help.begin(), help.end(),
976 std::back_inserter (q_points));
980 for (
unsigned int face=0; face<n_faces; ++face)
984 for (
unsigned int subface=0; subface<GeometryInfo<dim-1>::n_children(
RefinementCase<dim-1>(ref_case)); ++subface)
987 std::back_inserter (weights));
990 Assert (q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
992 Assert (weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
1004 const unsigned int child_no)
1009 const unsigned int n_q_points = quadrature.
size();
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);
1020 for (
unsigned int i=0; i<n_q_points; ++i)
1031 const unsigned int n_points = quadrature.
size(),
1034 std::vector<Point<dim> > q_points(n_points * n_children);
1035 std::vector<double>
weights(n_points * n_children);
1039 for (
unsigned int child=0; child<n_children; ++child)
1042 for (
unsigned int i=0; i<n_points; ++i)
1044 q_points[child*n_points+i] = help.
point(i);
1045 weights[child*n_points+i] = help.
weight(i);
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);
1066 for (
unsigned int k=0; k<n; ++k)
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);
1082 const bool face_orientation,
1083 const bool face_flip,
1084 const bool face_rotation,
1085 const unsigned int n_quadrature_points)
1094 return face_no * n_quadrature_points;
1119 static const unsigned int offset[2][2][2]=
1122 {6*GeometryInfo<dim>::faces_per_cell, 7*GeometryInfo<dim>::faces_per_cell}
1124 { {0*GeometryInfo<dim>::faces_per_cell, 1*GeometryInfo<dim>::faces_per_cell},
1125 {2*GeometryInfo<dim>::faces_per_cell, 3*GeometryInfo<dim>::faces_per_cell}
1130 + offset[face_orientation][face_flip][face_rotation])
1131 * n_quadrature_points);
1145 subface (
const unsigned int face_no,
1146 const unsigned int subface_no,
1150 const unsigned int n_quadrature_points,
1160 * n_quadrature_points);
1168 subface (
const unsigned int face_no,
1169 const unsigned int subface_no,
1173 const unsigned int n_quadrature_points,
1183 * n_quadrature_points);
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,
1198 const unsigned int dim = 3;
1229 const unsigned int total_subfaces_per_face=8;
1246 static const unsigned int orientation_offset[2][2][2]=
1289 static const unsigned int ref_case_offset[3]=
1388 static const unsigned int 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];
1433 final_ref_case = (face_orientation==face_rotation
1435 ref_case_permutation[equ_ref_case]
1474 return (((face_no * total_subfaces_per_face
1475 + ref_case_offset[final_ref_case-1]
1477 + orientation_offset[face_orientation][face_flip][face_rotation]
1479 * n_quadrature_points);
1487 const unsigned int face_no)
1489 std::vector<Point<dim> > points(quadrature.
size());
1490 project_to_face(quadrature, face_no, points);
1498 const unsigned int face_no,
1499 const unsigned int subface_no,
1502 std::vector<Point<dim> > points(quadrature.
size());
1503 project_to_subface(quadrature, face_no, subface_no, points, ref_case);
1515 bool at_left =
false,
1517 for (
unsigned int i=0; i<base_quadrature.
size(); ++i)
1525 return (at_left && at_right);
1635 const unsigned int )
1646 const unsigned int n_copies)
1649 (base_quadrature.
size()-1) * n_copies + 1 :
1650 base_quadrature.
size() * n_copies)
1656 if (!uses_both_endpoints(base_quadrature))
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)
1668 (1.0*copy)/n_copies);
1670 = base_quadrature.
weight(q_point) / n_copies;
1678 unsigned int next_point = 0;
1685 double double_point_weight = 0;
1686 unsigned int n_end_points = 0;
1687 for (
unsigned int i=0; i<base_quadrature.
size(); ++i)
1693 double_point_weight += base_quadrature.
weight(i);
1697 double_point_weight /= n_copies;
1702 Assert (n_end_points == 2, ExcInvalidQuadratureFormula());
1705 for (
unsigned int copy=0; copy<n_copies; ++copy)
1706 for (
unsigned int q_point=0; q_point<base_quadrature.
size(); ++q_point)
1719 (1.0*copy)/n_copies);
1726 if ((copy != n_copies-1) &&
1728 this->
weights[next_point] = double_point_weight;
1730 this->
weights[next_point] = base_quadrature.
weight(q_point) /
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,
1759 const unsigned int N)
1782 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
#define AssertDimension(dim1, dim2)
std::vector< double > weights
Quadrature(const unsigned int n_quadrature_points=0)
static Quadrature< 2 > rotate(const Quadrature< 2 > &q, const unsigned int n_times)
#define AssertIndexRange(index, range)
static bool uses_both_endpoints(const Quadrature< 1 > &base_quadrature)
static::ExceptionBase & ExcNotInitialized()
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
Quadrature & operator=(const Quadrature< dim > &)
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)
QAnisotropic(const Quadrature< 1 > &qx)
unsigned int size() const
static::ExceptionBase & ExcImpossibleInDim(int arg1)
void initialize(const std::vector< Point< dim > > &points, const std::vector< double > &weights)
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)
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)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double weight(const unsigned int i) const
std::vector< Point< dim > > quadrature_points
static Quadrature< 2 > reflect(const Quadrature< 2 > &q)
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)
std::size_t memory_consumption() const
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)
bool operator==(const Quadrature< dim > &p) const
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()
static::ExceptionBase & ExcZero()
static Quadrature< dim > project_to_all_children(const Quadrature< dim > &quadrature)
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)