]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
implement quadrature of dimension zero
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 30 Oct 2009 03:58:31 +0000 (03:58 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 30 Oct 2009 03:58:31 +0000 (03:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@20007 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/quadrature.h
deal.II/base/source/quadrature.cc

index 4db976bc853f35424049fae36481f3a2d542ed34..d9e944abeb8c3b804fdff0775c20e10a63a3f4f2 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -33,10 +33,20 @@ DEAL_II_NAMESPACE_OPEN
  * integration formulæ. Their names names prefixed by
  * <tt>Q</tt>. Refer to the list of derived classes for more details.
  *
- * The schemes for higher dimensions are tensor products of the
- * one-dimansional formulæ. Therefore, a three-dimensional 5-point
- * Gauss formula has 125 quadrature points.
+ * The schemes for higher dimensions are typically tensor products of the
+ * one-dimensional formulæ, but refer to the section on implementation
+ * detail below.
  *
+ * In order to allow for dimension independent programming, a
+ * quadrature formula of dimension zero exists. Since an integral over
+ * zero dimensions is the evaluation at a single point, any
+ * constructor of such a formula initializes to a single quadrature
+ * point with weight one. Access to the weight is possible, while
+ * access to the quadrature point is not permitted, since a Point of
+ * dimension zero contains no information. The main purpose of these
+ * formulæ is their use in QProjector, which will create a useful
+ * formula of dimension one out of them.
+ * 
  * <h3>Mathematical background</h3>
  *
  * For each quadrature formula we denote by <tt>m</tt>, the maximal
@@ -62,22 +72,10 @@ DEAL_II_NAMESPACE_OPEN
  * points in <tt>dim</tt> dimensions, where N is the constructor
  * parameter of QGauss.
  *
- * For some programs it is necessary to have a quadrature object for
- * faces.  These programs fail to link if compiled for only one space
- * dimension, since there quadrature rules for faces just don't make
- * no sense. In order to allow these programs to be linked anyway, for
- * class Quadrature@<0@> all functions are provided in the
- * <tt>quadrature.cc</tt> file, but they will throw exceptions if actually
- * called. The only function which is allowed to be called is the
- * constructor taking one integer, which in this case ignores its
- * parameter, and of course the destructor. Besides this, it is
- * necessary to provide a class Point@<0@> to make the compiler
- * happy. This class also does nothing.
+ * @note Instantiations for this template are provided for dimensions
+ * 0, 1, 2, and 3 (see the section on @ref Instantiations).
  *
- * @note Instantiations for this template are provided for dimensions 1, 2,
- * and 3 (see the section on @ref Instantiations in the manual).
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2005
+ * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2005, 2009
  */
 template <int dim>
 class Quadrature : public Subscriptor
@@ -413,11 +411,11 @@ Quadrature<dim>::weight (const unsigned int i) const
 template <>
 Quadrature<0>::Quadrature (const unsigned int);
 template <>
-Quadrature<0>::Quadrature (const Quadrature<0> &);
-template <>
 Quadrature<0>::Quadrature (const Quadrature<-1> &,
                           const Quadrature<1> &);
 template <>
+Quadrature<0>::Quadrature (const Quadrature<1> &);
+template <>
 Quadrature<0>::~Quadrature ();
 
 template <>
@@ -429,8 +427,6 @@ Quadrature<1>::Quadrature (const Quadrature<0> &);
 
 template <>
 const std::vector<Point<0> > & Quadrature<0>::get_points () const;
-template <>
-const std::vector<double> & Quadrature<0>::get_weights () const;
 
 #endif // DOXYGEN
 DEAL_II_NAMESPACE_CLOSE
index ea0443553894d8c27210bc3e7158809b82756c11..80d46e8449654a61b7798047e55148b43067d6c4 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -43,15 +43,8 @@ namespace
 
 template <>
 Quadrature<0>::Quadrature (const unsigned int)
-              : n_quadrature_points(0)
-{}
-
-
-
-template <>
-Quadrature<0>::Quadrature (const Quadrature<0> &)
-               : Subscriptor(),
-                 n_quadrature_points(0)
+               : n_quadrature_points(1),
+                 weights (1, 1.)
 {}
 
 
@@ -111,21 +104,9 @@ template <>
 Quadrature<0>::Quadrature (const Quadrature<-1> &,
                           const Quadrature<1> &)
                :
-                n_quadrature_points (0)
-{
-  Assert (false, ExcInternalError());
-}
-
-
-
-template <>
-Quadrature<1>::Quadrature (const Quadrature<0> &,
-                          const Quadrature<1> &)
-                :
-                n_quadrature_points (0)
-{
-  Assert (false, ExcInternalError());
-}
+                n_quadrature_points (1),
+               weights(1, 1.)
+{}
 
 
 
@@ -172,10 +153,55 @@ Quadrature<dim>::Quadrature (const SubQuadrature &q1,
 
 
 
+template <>
+Quadrature<1>::Quadrature (const SubQuadrature&,
+                          const Quadrature<1>& q2)
+               :
+               n_quadrature_points (q2.size()),
+               quadrature_points (n_quadrature_points),
+               weights (n_quadrature_points, 0)
+{
+  unsigned int present_index = 0;
+  for (unsigned int i2=0; i2<q2.size(); ++i2)
+    {
+                                        // compose coordinates of
+                                        // new quadrature point by tensor
+                                        // product in the last component
+      quadrature_points[present_index](0)
+       = q2.point(i2)(0);
+      
+      weights[present_index] = q2.weight(i2);
+      
+      ++present_index;
+    }
+
+#ifdef DEBUG
+  if (size() > 0)
+    {
+      double sum = 0;
+      for (unsigned int i=0; i<size(); ++i)
+       sum += weights[i];
+                                      // we cant guarantee the sum of weights
+                                      // to be exactly one, but it should be
+                                      // near that. 
+      Assert ((sum>0.999999) && (sum<1.000001), ExcInternalError());
+    }
+#endif
+}
+
+
+
+template <>
+Quadrature<0>::Quadrature (const Quadrature<1> &)
+               :
+               n_quadrature_points (1),
+               weights (1, 1.)
+{}
+
+
 template <>
 Quadrature<1>::Quadrature (const Quadrature<0> &)
                :
-               Subscriptor(),
                n_quadrature_points (numbers::invalid_unsigned_int),
                quadrature_points (),
                weights ()
@@ -279,16 +305,6 @@ Quadrature<dim>::get_weights () const
 
 
 
-template <>
-const std::vector<double> &
-Quadrature<0>::get_weights () const
-{
-  Assert (false, ExcInternalError());
-  return weights;
-}
-
-
-
 template <int dim>
 unsigned int
 Quadrature<dim>::memory_consumption () const
@@ -424,10 +440,14 @@ QProjector<dim>::rotate (const Quadrature<2> &q,
 template <>
 void
 QProjector<1>::project_to_face (const Quadrature<0> &,
-                                const unsigned int,
-                                std::vector<Point<1> > &)
+                                const unsigned int face_no,
+                                std::vector<Point<1> > &q_points)
 {
-  Assert (false, ExcNotImplemented());
+  const unsigned int dim=1;
+  Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
+  AssertDimension (q_points.size(), 1);
+  
+  q_points[0] = Point<dim>((double) face_no);
 }
 
 
@@ -520,12 +540,16 @@ QProjector<3>::project_to_face (const Quadrature<2>    &quadrature,
 template <>
 void
 QProjector<1>::project_to_subface (const Quadrature<0> &,
+                                   const unsigned int face_no,
                                    const unsigned int,
-                                   const unsigned int,
-                                   std::vector<Point<1> > &,
+                                   std::vector<Point<1> >& q_points,
                                   const RefinementCase<0> &)
 {
-  Assert(false, ExcNotImplemented());
+  const unsigned int dim=1;
+  Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
+  AssertDimension (q_points.size(), 1);
+  
+  q_points[0] = Point<dim>((double) face_no);
 }
 
 
@@ -697,10 +721,42 @@ QProjector<3>::project_to_subface (const Quadrature<2>    &quadrature,
 
 template <>
 Quadrature<1>
-QProjector<1>::project_to_all_faces (const Quadrature<0> &)
+QProjector<1>::project_to_all_faces (const Quadrature<0>& quadrature)
 {
-  Assert (false, ExcImpossibleInDim(1));
-  return Quadrature<1>(0);
+  const unsigned int dim = 1;
+  
+  const unsigned int n_points = 1,
+                    n_faces  = GeometryInfo<dim>::faces_per_cell;
+
+                                  // first fix quadrature points
+  std::vector<Point<dim> > q_points;
+  q_points.reserve(n_points * n_faces);
+  std::vector <Point<dim> > help(n_points);
+  
+  
+                                  // project to each face and append
+                                  // results
+  for (unsigned int face=0; face<n_faces; ++face)
+    {
+      project_to_face(quadrature, face, help);
+      std::copy (help.begin(), help.end(),
+                 std::back_inserter (q_points));
+    }
+
+                                  // next copy over weights
+  std::vector<double> weights;
+  weights.reserve (n_points * n_faces);
+  for (unsigned int face=0; face<n_faces; ++face)
+    std::copy (quadrature.get_weights().begin(),
+               quadrature.get_weights().end(),
+               std::back_inserter (weights));
+
+  Assert (q_points.size() == n_points * n_faces,
+          ExcInternalError());
+  Assert (weights.size() == n_points * n_faces,
+          ExcInternalError());  
+  
+  return Quadrature<dim>(q_points, weights);
 }
 
 
@@ -811,10 +867,44 @@ QProjector<3>::project_to_all_faces (const SubQuadrature &quadrature)
 
 template <>
 Quadrature<1>
-QProjector<1>::project_to_all_subfaces (const Quadrature<0> &)
+QProjector<1>::project_to_all_subfaces (const Quadrature<0>& quadrature)
 {
-  Assert (false, ExcImpossibleInDim(1));
-  return Quadrature<1>(0);
+  const unsigned int dim = 1;
+  
+  const unsigned int n_points          = 1,
+                    n_faces           = GeometryInfo<dim>::faces_per_cell,
+                    subfaces_per_face = GeometryInfo<dim>::max_children_per_face;
+  
+                                  // first fix quadrature points
+  std::vector<Point<dim> > q_points;
+  q_points.reserve (n_points * n_faces * subfaces_per_face);
+  std::vector <Point<dim> > help(n_points);
+  
+                                  // 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, face, subface, help);
+       std::copy (help.begin(), help.end(),
+                   std::back_inserter (q_points));
+      };
+
+                                  // next copy over weights
+  std::vector<double> weights;
+  weights.reserve (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(),
+                 std::back_inserter (weights));
+
+  Assert (q_points.size() == n_points * n_faces * subfaces_per_face,
+          ExcInternalError());
+  Assert (weights.size() == n_points * n_faces * subfaces_per_face,
+          ExcInternalError());
+  
+  return Quadrature<dim>(q_points, weights);
 }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.