]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Check in a class which is able to iterate a given quadrature formula n times in each...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 5 Jun 1999 14:05:57 +0000 (14:05 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 5 Jun 1999 14:05:57 +0000 (14:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@1369 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 5c78399d41a6e48482c86bd93d43a17d2999f9e5..6a17905eff6a8f9bc6f13d398086a0418f0637ba 100644 (file)
@@ -126,6 +126,59 @@ class Quadrature {
 
 
 
+/**
+ * Quadrature formula constructed by iteration of another quadrature formula in
+ * each direction. In more than one space dimension, the resulting quadrature
+ * formula is constructed in the usual way by building the tensor product of
+ * the respective iterated quadrature formula in one space dimension.
+ *
+ * In one space dimension, the given base formula is copied and scaled onto
+ * a given number of subintervals of length #1/n_copies#. If the quadrature
+ * formula uses both end points of theunit interval, then in the interior
+ * of the iterated quadrature formula there would be quadrature points which
+ * are used twice; we merge them into one with a weight which is the sum
+ * of the weights of the left- and the rightmost quadrature point.
+ *
+ * Since all dimensions higher than one are built up by tensor products of
+ * one dimensional and #dim-1# dimensional quadrature formulae, the
+ * argument given to the constructor needs to be a quadrature formula in
+ * one space dimension, rather than in #dim# dimensions.
+ *
+ * @author Wolfgang Bangerth 1999
+ */
+template <int dim>
+class QIterated : public Quadrature<dim>
+{
+  public:
+                                    /**
+                                     * Constructor. Iterate the given
+                                     * quadrature formula $n_copies$ times in
+                                     * each direction.
+                                     */
+    QIterated (const Quadrature<1> &base_quadrature,
+              const unsigned int   n_copies);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcSumOfWeightsNotOne);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidQuadratureFormula);
+    
+  private:
+                                    /**
+                                     * check whether the given
+                                     * quadrature formula has quadrature
+                                     * points at the left and right end points
+                                     * of the interval
+                                     */
+    static bool uses_both_endpoints (const Quadrature<1> &base_quadrature);
+};
+
+
+
 
 /**
  *  This class is a helper class to facilitate the usage of quadrature formulae
index fe33ddfb4ba0eca463097188bcf48dc5fff281e5..6c70bc76f9f2026326e602ee994418d6744765c5 100644 (file)
@@ -3,6 +3,8 @@
 
 
 #include <base/quadrature.h>
+#include <cmath>
+
 
 
 template <>
@@ -448,7 +450,134 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
 
 
 
+template <>
+bool
+QIterated<1>::uses_both_endpoints (const Quadrature<1> &base_quadrature)
+{
+  bool at_left = false,
+      at_right = false;
+  for (unsigned int i=0; i<base_quadrature.n_quadrature_points; ++i)
+    {
+      if (base_quadrature.quad_point(i) == Point<1>(0.0))
+       at_left = true;
+      if (base_quadrature.quad_point(i) == Point<1>(1.0))
+       at_right = true;
+    };
+
+  return (at_left && at_right);
+};
+
+
+
+template <>
+QIterated<1>::QIterated (const Quadrature<1> &base_quadrature,
+                        const unsigned int   n_copies) :
+               Quadrature<1> (uses_both_endpoints(base_quadrature) ?
+                              (base_quadrature.n_quadrature_points-1) * n_copies + 1 :
+                              base_quadrature.n_quadrature_points * n_copies) 
+{
+                                  // we don't have to skip some points
+                                  // in order to get a reasonable quadrature
+                                  // formula
+  if (!uses_both_endpoints(base_quadrature))
+    {
+      unsigned int next_point = 0;
+      for (unsigned int copy=0; copy<n_copies; ++copy)
+       for (unsigned int q_point=0; q_point<base_quadrature.n_quadrature_points; ++q_point)
+         {
+           quadrature_points[next_point] = Point<1>(base_quadrature.quad_point(q_point)(0) / n_copies
+                                                    +
+                                                    (1.0*copy)/n_copies);
+           weights[next_point]           = base_quadrature.weight(q_point) / n_copies;
+
+           ++next_point;
+         };
+    }
+  else
+                                    // skip doubly available points
+    {
+      unsigned int next_point = 0;
+
+                                      // first find out the weights of
+                                      // the left and the right boundary
+                                      // points. note that these usually
+                                      // are but need not necessarily be
+                                      // the same
+      double double_point_weight = 0;
+      unsigned int n_end_points = 0;
+      for (unsigned int i=0; i<base_quadrature.n_quadrature_points; ++i)
+                                        // add up the weight if this
+                                        // is an endpoint
+       if ((base_quadrature.quad_point(i) == Point<1>(0.0)) ||
+           (base_quadrature.quad_point(i) == Point<1>(1.0)))
+         {
+           double_point_weight += base_quadrature.weight(i);
+           ++n_end_points;
+         };
+                                      // scale the weight correctly
+      double_point_weight /= n_copies;
+
+                                      // make sure the base quadrature formula
+                                      // has only one quadrature point
+                                      // per end point
+      Assert (n_end_points == 2, ExcInvalidQuadratureFormula());
+      
+      
+      for (unsigned int copy=0; copy<n_copies; ++copy)
+       for (unsigned int q_point=0; q_point<base_quadrature.n_quadrature_points; ++q_point)
+         {
+                                            // skip the left point of this copy
+                                            // since we have already entered it the
+                                            // last time
+           if ((copy > 0) &&
+               (base_quadrature.quad_point(q_point) == Point<1>(0.0)))
+             continue;
+           
+           quadrature_points[next_point] = Point<1>(base_quadrature.quad_point(q_point)(0) / n_copies
+                                                    +
+                                                    (1.0*copy)/n_copies);
+
+                                            // if this is the rightmost point of
+                                            // one of the non-last copies: give
+                                            // it the double weight
+           if ((copy != n_copies-1) &&
+               (base_quadrature.quad_point(q_point) == Point<1>(1.0)))
+             weights[next_point] = double_point_weight;
+           else
+             weights[next_point] = base_quadrature.weight(q_point) / n_copies;
+           
+           ++next_point;
+         };
+    };
+
+#if DEBUG
+  double sum_of_weights = 0;
+  for (unsigned int i=0; i<n_quadrature_points; ++i)
+    sum_of_weights += weight(i);
+  Assert (fabs(sum_of_weights-1) < 1e-15,
+         ExcSumOfWeightsNotOne());
+#endif
+};
+
+
+
+
+// construct higher dimensional quadrature formula by tensor product of lower
+// dimensional iterated quadrature formulae
+template <int dim>
+QIterated<dim>::QIterated (const Quadrature<1> &base_quadrature,
+                          const unsigned int   N):
+               Quadrature<dim> (QIterated<dim-1>(base_quadrature, N),
+                                QIterated<1>(base_quadrature, N))
+{};
+
+
+
+
 // explicit instantiations; note: we need them all for all dimensions
 template class Quadrature<1>;
 template class Quadrature<2>;
 template class Quadrature<3>;
+template class QIterated<1>;
+template class QIterated<2>;
+template class QIterated<3>;

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.