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