From dc8a40fc0506f39dce3e307d0cca890c77505248 Mon Sep 17 00:00:00 2001 From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Wed, 8 Aug 2001 08:31:19 +0000 Subject: [PATCH] Remove the kludge that the QProjector happened to become a real class derived from Quadrature. Rather, make the respective constructor call two functions in their own that return a Quadrature object. git-svn-id: https://svn.dealii.org/trunk@4869 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature.h | 111 ++++++++++------- deal.II/base/source/quadrature.cc | 116 ++++++++++-------- deal.II/deal.II/source/fe/fe.cc | 8 +- .../deal.II/source/fe/mapping_cartesian.cc | 8 +- deal.II/deal.II/source/fe/mapping_q.cc | 4 +- deal.II/deal.II/source/fe/mapping_q1.cc | 12 +- deal.II/deal.II/source/grid/grid_out.cc | 6 +- tests/base/quadrature_test.cc | 4 +- 8 files changed, 155 insertions(+), 114 deletions(-) diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index 53802bce32..56bb445fbc 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -276,55 +276,30 @@ class QIterated : public Quadrature<dim> * for a description of the orientation of the different faces. */ template <int dim> -class QProjector : public Quadrature<dim> +class QProjector { public: -//TODO:[GK] remove this constructor again and return to old state with static-only functions. move this function to something in the finite element classes which is the only place where it is used. - /** - * Constructor for a quadrature rule on all (sub)faces. - * The quadrature rule - * @p{quadrature} is applied to - * each face or subface, according - * to the @{sub} flag. - * - * The weights of the new rule are - * replications of the original - * weights. This is not a proper - * handling, but it is - * consistent with later use. If - * the (sub)face rule is applied to - * the unity function, the result - * is the number of (sub)faces. - */ - QProjector (const Quadrature<dim-1>& quadrature, - const bool sub); - + /** - * Compute the quadrature points on the - * cell if the given quadrature formula - * is used on face @p{face_no}. For further - * details, see the general doc for - * this class. + * Compute the quadrature points + * on the cell if the given + * quadrature formula is used on + * face @p{face_no}. For further + * details, see the general doc + * for this class. */ static void project_to_face (const Quadrature<dim-1> &quadrature, const unsigned int face_no, typename std::vector<Point<dim> > &q_points); - /** - * Projection to all faces. - * Generate a formula that integrates - * over all faces at the same time. - */ - static void project_to_faces (const Quadrature<dim-1> &quadrature, - typename std::vector<Point<dim> > &q_points); - /** - * Compute the quadrature points on the - * cell if the given quadrature formula - * is used on face @p{face_no}, subface - * number @p{subface_no}. For further - * details, see the general doc for - * this class. + * Compute the quadrature points + * on the cell if the given + * quadrature formula is used on + * face @p{face_no}, subface + * number @p{subface_no}. For + * further details, see the + * general doc for this class. */ static void project_to_subface (const Quadrature<dim-1> &quadrature, const unsigned int face_no, @@ -332,14 +307,46 @@ class QProjector : public Quadrature<dim> typename std::vector<Point<dim> > &q_points); /** - * Projection to all child faces. - * Project to the children of all - * faces at the same time. The - * ordering is first by face, - * then by subface + * Take a face quadrature formula + * and generate a cell quadrature + * formula from it where the + * quadrature points of the given + * argument are projected on all + * faces. + * + * The weights of the new rule + * are replications of the + * original weights. This is not + * a proper handling, in that the + * sum of weights does not equal + * one, but it is consistent with + * the use of this function, + * namely to generate sets of + * face quadrature points on a + * cell, one set of which will + * then be selected at each + * time. This is used in the + * @ref{FEFaceValues} class, + * where we initialize the + * values, derivatives, etc on + * all faces at once, while + * selecting the data of one + * particular face only happens + * later. + */ + static Quadrature<dim> + project_to_all_faces (const Quadrature<dim-1> &quadrature); + + /** + * This function is alike the + * previous one, but projects the + * given face quadrature formula + * to the subfaces of a cell, + * i.e. to the children of the + * faces of the unit cell. */ - static void project_to_subfaces (const Quadrature<dim-1> &quadrature, - typename std::vector<Point<dim> > &q_points); + static Quadrature<dim> + project_to_all_subfaces (const Quadrature<dim-1> &quadrature); }; @@ -374,6 +381,11 @@ template <> void QProjector<3>::project_to_face (const Quadrature<2> &quadrature, const unsigned int face_no, std::vector<Point<3> > &q_points); + +template <> +Quadrature<1> QProjector<1>::project_to_all_faces (const Quadrature<0> &quadrature); + + template <> void QProjector<1>::project_to_subface (const Quadrature<0> &, const unsigned int, @@ -389,6 +401,11 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector<Point<3> > &q_points); + +template <> +Quadrature<1> QProjector<1>::project_to_all_subfaces (const Quadrature<0> &quadrature); + + template <> bool QIterated<1>::uses_both_endpoints (const Quadrature<1> &base_quadrature); diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index 18ae378a93..62ab8d2fd9 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -222,27 +222,6 @@ Quadrature<dim>::memory_consumption () const //----------------------------------------------------------------------// -template <int dim> -QProjector<dim>::QProjector (const Quadrature<dim-1> &quadrature, - const bool sub) - : - Quadrature<dim> (quadrature.n_quadrature_points - * GeometryInfo<dim>::faces_per_cell - * (sub ? GeometryInfo<dim>::subfaces_per_face : 1)) -{ - if (sub) - project_to_subfaces (quadrature, quadrature_points); - else - project_to_faces (quadrature, quadrature_points); - - const unsigned int n = GeometryInfo<dim>::faces_per_cell - * (sub ? GeometryInfo<dim>::subfaces_per_face : 1); - unsigned int k=0; - for (unsigned int i=0; i<n; ++i) - for (unsigned int j=0; j<quadrature.n_quadrature_points; ++j) - weights[k++] = quadrature.weight(j); -} - template <> void QProjector<1>::project_to_face (const Quadrature<0> &, @@ -576,49 +555,88 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, +template <> +Quadrature<1> +QProjector<1>::project_to_all_faces (const Quadrature<0> &) +{ + Assert (false, ExcImpossibleInDim(1)); + return Quadrature<1>(0); +}; + + + template <int dim> -void -QProjector<dim>::project_to_faces (const Quadrature<dim-1> &quadrature, - typename std::vector<Point<dim> > &q_points) +Quadrature<dim> +QProjector<dim>::project_to_all_faces (const Quadrature<dim-1> &quadrature) { - unsigned int npt = quadrature.n_quadrature_points; - unsigned int nf = GeometryInfo<dim>::faces_per_cell; - - q_points.resize (npt*nf); - std::vector <Point<dim> > help(npt); + const unsigned int n_points = quadrature.n_quadrature_points, + n_faces = GeometryInfo<dim>::faces_per_cell; + + // first fix quadrature points + typename std::vector<Point<dim> > q_points (n_points * n_faces); + std::vector <Point<dim> > help(n_points); - unsigned k=0; - for (unsigned int i=0;i<nf;++i) + // project to each face and copy + // results + for (unsigned int face=0; face<n_faces; ++face) { - project_to_face(quadrature, i, help); - for (unsigned int j=0;j<npt;++j) - q_points[k++] = help[j]; + project_to_face(quadrature, face, help); + std::copy (help.begin(), help.end(), q_points.begin()+n_points*face); } + + // next copy over weights + typename std::vector<double> weights (n_points * n_faces); + for (unsigned int face=0; face<n_faces; ++face) + std::copy (quadrature.get_weights().begin(), + quadrature.get_weights().end(), + weights.begin()+n_points*face); + + return Quadrature<dim>(q_points, weights); } +template <> +Quadrature<1> +QProjector<1>::project_to_all_subfaces (const Quadrature<0> &) +{ + Assert (false, ExcImpossibleInDim(1)); + return Quadrature<1>(0); +}; + + template <int dim> -void -QProjector<dim>::project_to_subfaces (const Quadrature<dim-1> &quadrature, - typename std::vector<Point<dim> > &q_points) +Quadrature<dim> +QProjector<dim>::project_to_all_subfaces (const Quadrature<dim-1> &quadrature) { - unsigned int npt = quadrature.n_quadrature_points; - unsigned int nf = GeometryInfo<dim>::faces_per_cell; - unsigned int nc = GeometryInfo<dim>::subfaces_per_face; + const unsigned int n_points = quadrature.n_quadrature_points, + n_faces = GeometryInfo<dim>::faces_per_cell, + subfaces_per_face = GeometryInfo<dim>::subfaces_per_face; - q_points.resize (npt*nf*nc); - std::vector <Point<dim> > help(npt); + // first fix quadrature points + typename std::vector<Point<dim> > q_points (n_points * n_faces * subfaces_per_face); + std::vector <Point<dim> > help(n_points); - unsigned k=0; - for (unsigned int i=0;i<nf;++i) - for (unsigned int c=0;c<nc;++c) + // project to each face and copy + // results + for (unsigned int face=0; face<n_faces; ++face) + for (unsigned int subface=0; subface<subfaces_per_face; ++subface) { - project_to_subface(quadrature, i, c, help); - for (unsigned int j=0;j<npt;++j) - q_points[k++] = help[j]; - } + project_to_subface(quadrature, face, subface, help); + std::copy (help.begin(), help.end(), + q_points.begin()+n_points*(face*subfaces_per_face+subface)); + }; + + // next copy over weights + typename std::vector<double> weights (n_points * n_faces * subfaces_per_face); + for (unsigned int face=0; face<n_faces; ++face) + for (unsigned int subface=0; subface<subfaces_per_face; ++subface) + std::copy (quadrature.get_weights().begin(), + quadrature.get_weights().end(), + weights.begin()+n_points*(face*subfaces_per_face+subface)); + + return Quadrature<dim>(q_points, weights); } diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index d7b070aa13..8b16e29b7a 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -424,8 +424,8 @@ FiniteElement<dim>::get_face_data (const UpdateFlags flags, const Mapping<dim> &mapping, const Quadrature<dim-1> &quadrature) const { - QProjector<dim> q(quadrature, false); - return get_data (flags, mapping, q); + return get_data (flags, mapping, + QProjector<dim>::project_to_all_faces(quadrature)); } @@ -436,8 +436,8 @@ FiniteElement<dim>::get_subface_data (const UpdateFlags flags, const Mapping<dim> &mapping, const Quadrature<dim-1> &quadrature) const { - QProjector<dim> q(quadrature, true); - return get_data (flags, mapping, q); + return get_data (flags, mapping, + QProjector<dim>::project_to_all_subfaces(quadrature)); } diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 4516ad50e4..45de9489d5 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -99,8 +99,8 @@ typename Mapping<dim>::InternalDataBase * MappingCartesian<dim>::get_face_data (const UpdateFlags update_flags, const Quadrature<dim-1>& quadrature) const { - QProjector<dim> q (quadrature, false); - InternalData* data = new InternalData (q); + InternalData* data + = new InternalData (QProjector<dim>::project_to_all_faces(quadrature)); data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); @@ -116,8 +116,8 @@ typename Mapping<dim>::InternalDataBase * MappingCartesian<dim>::get_subface_data (const UpdateFlags update_flags, const Quadrature<dim-1> &quadrature) const { - QProjector<dim> q (quadrature, true); - InternalData* data = new InternalData (q); + InternalData* data + = new InternalData (QProjector<dim>::project_to_all_subfaces(quadrature)); data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index ae159135f5..926440c822 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -221,7 +221,7 @@ MappingQ<dim>::get_face_data (const UpdateFlags update_flags, const Quadrature<dim-1>& quadrature) const { InternalData *data = new InternalData(n_shape_functions); - QProjector<dim> q (quadrature, false); + Quadrature<dim> q (QProjector<dim>::project_to_all_faces(quadrature)); compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); if (!use_mapping_q_on_all_cells) @@ -239,7 +239,7 @@ MappingQ<dim>::get_subface_data (const UpdateFlags update_flags, const Quadrature<dim-1>& quadrature) const { InternalData *data = new InternalData(n_shape_functions); - QProjector<dim> q (quadrature, true); + Quadrature<dim> q (QProjector<dim>::project_to_all_subfaces(quadrature)); compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); if (!use_mapping_q_on_all_cells) diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 2b5837bea7..1a0e665d11 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -462,8 +462,10 @@ MappingQ1<dim>::get_face_data (const UpdateFlags update_flags, const Quadrature<dim-1> &quadrature) const { InternalData* data = new InternalData(n_shape_functions); - QProjector<dim> q (quadrature, false); - compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); + compute_face_data (update_flags, + QProjector<dim>::project_to_all_faces(quadrature), + quadrature.n_quadrature_points, + *data); return data; } @@ -476,8 +478,10 @@ MappingQ1<dim>::get_subface_data (const UpdateFlags update_flags, const Quadrature<dim-1>& quadrature) const { InternalData* data = new InternalData(n_shape_functions); - QProjector<dim> q (quadrature, true); - compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); + compute_face_data (update_flags, + QProjector<dim>::project_to_all_subfaces(quadrature), + quadrature.n_quadrature_points, + *data); return data; } diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index 9a02433b87..58424c7641 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -279,7 +279,7 @@ void GridOut::write_gnuplot (const Triangulation<dim> &tria, // quadrature formula which will be // used to probe boundary points at // curved faces - QProjector<dim> *q_projector=0; + Quadrature<dim> *q_projector=0; if (mapping!=0) { typename std::vector<Point<dim-1> > boundary_points(n_points); @@ -287,7 +287,7 @@ void GridOut::write_gnuplot (const Triangulation<dim> &tria, boundary_points[i](0)= 1.*(i+1)/(n_points+1); Quadrature<dim-1> quadrature(boundary_points); - q_projector = new QProjector<dim> (quadrature, false); + q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature)); } for (; cell!=endc; ++cell) @@ -602,7 +602,7 @@ void GridOut::write_eps (const Triangulation<dim> &tria, boundary_points[i](0) = 1.*(i+1)/(n_points+1); Quadrature<dim-1> quadrature (boundary_points); - QProjector<dim> q_projector (quadrature, false); + Quadrature<dim> q_projector (QProjector<dim>::project_to_all_faces(quadrature)); // next loop over all // boundary faces and diff --git a/tests/base/quadrature_test.cc b/tests/base/quadrature_test.cc index c7d3d887ad..2ab4caa994 100644 --- a/tests/base/quadrature_test.cc +++ b/tests/base/quadrature_test.cc @@ -118,7 +118,9 @@ check_faces (const std::vector<Quadrature<dim-1>*>& quadratures, const bool sub) for (unsigned int n=0; n<quadratures.size(); ++n) { - QProjector<dim> quadrature(*quadratures[n], sub); + Quadrature<dim> quadrature (sub == false? + QProjector<dim>::project_to_all_faces(*quadratures[n]) : + QProjector<dim>::project_to_all_subfaces(*quadratures[n])); const std::vector<Point<dim> > &points=quadrature.get_points(); const std::vector<double> &weights=quadrature.get_weights(); -- 2.39.5