QProjector::project_to_child

git-svn-id: https://svn.dealii.org/trunk@4906 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h
index 56bb445fbcf08963277c942a3501b3bcadbf3127..5a9d61f3f5f8241eb4b069882f6c7749daba7eb3 100644 (file)
--- a/deal.II/base/include/base/quadrature.h
+++ b/deal.II/base/include/base/quadrature.h
@@ -279,7 +279,6 @@ template <int dim>
 class QProjector
 {
   public:
-    
                                     /**
                                      * Compute the quadrature points
                                      * on the cell if the given
@@ -347,6 +346,30 @@ class QProjector
                                      */
     static Quadrature<dim>
     project_to_all_subfaces (const Quadrature<dim-1> &quadrature);
+
+                                    /**
+                                     * Project a give quadrature
+                                     * formula to a child of a
+                                     * cell. You may want to use this
+                                     * function in case you want to
+                                     * extend an integral only over
+                                     * the area which a potential
+                                     * child would occupy. The child
+                                     * numbering is the same as the
+                                     * children would be numbered
+                                     * upon refinement of the cell.
+                                     *
+                                     * As integration using this
+                                     * quadrature formula now only
+                                     * extends over a fraction of the
+                                     * cell, the weights of the
+                                     * resulting object are scaled by
+                                     * @p{1./GeometryInfo<dim>::children_per_cell}.
+                                     */
+    static
+    Quadrature<dim>
+    project_to_child (const Quadrature<dim>  &quadrature,
+                     const unsigned int      child_no);
 };
 
 
@@ -369,6 +392,7 @@ template <>
 double Quadrature<0>::weight (const unsigned int) const;
 template <>
 const std::vector<double> & Quadrature<0>::get_weights () const;
+
 template <>
 void QProjector<1>::project_to_face (const Quadrature<0> &,
                                     const unsigned int,

diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc
index 62ab8d2fd9eed242332e7985c3b1e80e5b6bec9d..63e222cbc40588298d71410a42fa81855d640683 100644 (file)
--- a/deal.II/base/source/quadrature.cc
+++ b/deal.II/base/source/quadrature.cc
@@ -240,6 +240,8 @@ void QProjector<2>::project_to_face (const Quadrature<1>      &quadrature,
 {
   const unsigned int dim=2;
   Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
+  Assert (q_points.size() == quadrature.n_quadrature_points,
+         ExcDimensionMismatch (q_points.size(), quadrature.n_quadrature_points));
   
   for (unsigned int p=0; p<quadrature.n_quadrature_points; ++p)
     switch (face_no)
@@ -270,6 +272,8 @@ void QProjector<3>::project_to_face (const Quadrature<2>    &quadrature,
 {
   const unsigned int dim=3;
   Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
+  Assert (q_points.size() == quadrature.n_quadrature_points,
+         ExcDimensionMismatch (q_points.size(), quadrature.n_quadrature_points));
   
   for (unsigned int p=0; p<quadrature.n_quadrature_points; ++p)
     switch (face_no)
@@ -332,6 +336,8 @@ void QProjector<2>::project_to_subface (const Quadrature<1>    &quadrature,
   const unsigned int dim=2;
   Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
   Assert (subface_no<(1<<(dim-1)), ExcIndexRange (face_no, 0, 1<<(dim-1)));
+  Assert (q_points.size() == quadrature.n_quadrature_points,
+         ExcDimensionMismatch (q_points.size(), quadrature.n_quadrature_points));
   
   for (unsigned int p=0; p<quadrature.n_quadrature_points; ++p)
     switch (face_no)
@@ -406,6 +412,8 @@ void QProjector<3>::project_to_subface (const Quadrature<2>    &quadrature,
   const unsigned int dim=3;
   Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
   Assert (subface_no<(1<<(dim-1)), ExcIndexRange (face_no, 0, 1<<(dim-1)));
+  Assert (q_points.size() == quadrature.n_quadrature_points,
+         ExcDimensionMismatch (q_points.size(), quadrature.n_quadrature_points));
 
 
                                   // for all faces and subfaces:
@@ -641,6 +649,60 @@ QProjector<dim>::project_to_all_subfaces (const Quadrature<dim-1> &quadrature)
 
 
 
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_child (const Quadrature<dim>    &quadrature,
+                                  const unsigned int        child_no)
+{
+  Assert (child_no < GeometryInfo<dim>::children_per_cell,
+         ExcIndexRange (child_no, 0, GeometryInfo<dim>::children_per_cell));
+  
+  const unsigned int n_q_points = quadrature.n_quadrature_points;  
+    
+                                  // projection is simple: copy over
+                                  // old points, scale them, and if
+                                  // necessary shift them
+  std::vector<Point<dim> > q_points = quadrature.get_points ();
+  for (unsigned int i=0; i<n_q_points; ++i)
+    q_points[i] /= 2;
+
+  switch (dim)
+    {
+      case 1:
+      {
+       if (child_no == 1)
+         for (unsigned int i=0; i<n_q_points; ++i)
+           q_points[i][0] += 0.5;
+       break;
+      };
+
+      case 2:
+      {
+       if ((child_no == 1) || (child_no == 2))
+         for (unsigned int i=0; i<n_q_points; ++i)
+           q_points[i][0] += 0.5;
+       if ((child_no == 2) || (child_no == 3))
+         for (unsigned int i=0; i<n_q_points; ++i)
+           q_points[i][1] += 0.5;
+       break;
+      };
+
+      default:
+           Assert (false, ExcNotImplemented());
+    };
+
+                                  // for the weights, things are
+                                  // equally simple: copy them and
+                                  // scale them
+  std::vector<double> weights = quadrature.get_weights ();
+  for (unsigned int i=0; i<n_q_points; ++i)
+    weights[i] *= (1./GeometryInfo<dim>::children_per_cell);
+
+  return Quadrature<dim> (q_points, weights);
+};
+
+
+
 // ------------------------------------------------------------ //
 
 

diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index 77bb3a4f4474473fe761dbd48633670ee7a71da2..6da9f8736d440753c07f7fb0d7be6d62e9d0c70e 100644 (file)
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -158,6 +158,13 @@ documentation, etc</a>.
 <h3>base</h3>
 
 <ol>
+  <li> <p> 
+       New: Function <code class="member">QProjector::project_to_child</code>
+       generates quadrature formulae which act on the area which a
+       child would occupy.
+       <br>
+       (WB 2001/08/23)
+
   <li> <p>
        Changed: The examples classes in the base directory are now
        moved into a namespace <code class="class">Functions</code> of