[quadraturerules]

add GaussLobatto rules and update writequad.mac to match the new indentation

README

@@ -70,4 +70,3 @@ for further options.
 The full build-system is described in the dune-common/doc/buildsystem (SVN version) or under share/doc/dune-common/buildsystem if you installed DUNE!
 
 $Id: duneproject 5232 2008-07-03 09:34:11Z christi $
-

dune/geometry/quadraturerules.hh

@@ -75,11 +75,7 @@ namespace Dune {
       Jacobian_1_0 = 1,
       Jacobian_2_0 = 2,
 
-      Simpson    = 3,
-      Trap       = 4,
-      Grid       = 5,
-
-      Clough     = 21,
+      GaussLobatto = 4,
 
       Invalid_Rule = 127
     };
@@ -410,6 +406,73 @@ namespace Dune {
 #define DUNE_INCLUDING_IMPLEMENTATION
 #include "quadraturerules/jacobi_2_0_imp.hh"
 
+namespace Dune {
+
+  //! \internal Helper template for the initialization of the quadrature rules
+  template<typename ct,
+      bool fundamental = std::numeric_limits<ct>::is_specialized>
+  struct GaussLobattoQuadratureInitHelper;
+  template<typename ct>
+  struct GaussLobattoQuadratureInitHelper<ct, true> {
+    static void init(int p,
+                     std::vector< FieldVector<ct, 1> > & _points,
+                     std::vector< ct > & _weight,
+                     int & delivered_order);
+  };
+  template<typename ct>
+  struct GaussLobattoQuadratureInitHelper<ct, false> {
+    static void init(int p,
+                     std::vector< FieldVector<ct, 1> > & _points,
+                     std::vector< ct > & _weight,
+                     int & delivered_order);
+  };
+
+  /** \brief Jacobi-Gauss quadrature for alpha=2, beta=0
+      \ingroup Quadrature
+   */
+  template<typename ct>
+  class GaussLobattoQuadratureRule1D :
+    public QuadratureRule<ct,1>
+  {
+  public:
+    /** \brief The space dimension */
+    enum { dim=1 };
+
+    /** \brief The highest quadrature order available */
+    enum { highest_order=61 };
+
+    ~GaussLobattoQuadratureRule1D(){}
+  private:
+    friend class QuadratureRuleFactory<ct,dim>;
+    GaussLobattoQuadratureRule1D (int p)
+      : QuadratureRule<ct,1>(GeometryType(GeometryType::cube, 1))
+    {
+      //! set up quadrature of given order in d dimensions
+      std::vector< FieldVector<ct, dim> > _points;
+      std::vector< ct > _weight;
+
+      int delivered_order;
+
+      GaussLobattoQuadratureInitHelper<ct>::init
+        (p, _points, _weight, delivered_order);
+
+      this->delivered_order = delivered_order;
+      assert(_points.size() == _weight.size());
+      for (size_t i = 0; i < _points.size(); i++)
+        this->push_back(QuadraturePoint<ct,dim>(_points[i], _weight[i]));
+    }
+  };
+
+#ifndef DOXYGEN
+  extern template GaussLobattoQuadratureRule1D<float>::GaussLobattoQuadratureRule1D(int);
+  extern template GaussLobattoQuadratureRule1D<double>::GaussLobattoQuadratureRule1D(int);
+#endif // !DOXYGEN
+
+} // namespace Dune
+
+#define DUNE_INCLUDING_IMPLEMENTATION
+#include "quadraturerules/gausslobatto_imp.hh"
+
 #include "quadraturerules/genericquadrature.hh"
 
 #include "quadraturerules/simplexquadrature.hh"
@@ -640,6 +703,8 @@ namespace Dune {
           return Jacobi1QuadratureRule1D<ctype>(p);
         case QuadratureType::Jacobian_2_0 :
           return Jacobi2QuadratureRule1D<ctype>(p);
+        case QuadratureType::GaussLobatto :
+          return GaussLobattoQuadratureRule1D<ctype>(p);
         default :
           DUNE_THROW(Exception, "Unknown QuadratureType");
         }

dune/geometry/quadraturerules/gausslobatto.mac

+/*
+ * Run `maxima --batch jacobian.mac'.
+ */
+
+/* adjust search path */
+file_search_maxima: append(["$$$.{mac,max}"],file_search_maxima);
+
+load(orthopoly)$
+load(stringproc)$
+load(writequad)$
+
+orthopoly_returns_intervals : false;
+
+fpprec: 120$
+fpprintprec: 100$
+maxp: 30$
+/* float2bf: false$ */
+
+/* points are given as the roots of (x^2-1)*L_p(x) */
+point(n) := (x^2-1)*diff(legendre_p(n,x),x)$
+weight(n,xx) := ratsimp(2/((n*(n+1))*legendre_p(n,xx)^2))$
+
+/*
+  write files
+*/
+
+/* num-points 4 */
+write_quad("gausslobatto_imp.hh", "GaussLobatto", "gausslobatto.mac",
+		   lambda([i],point(i)),
+           lambda([i,p],weight(i,p)),
+		   maxp)$
+
+/*
+ * Local variables:
+ * mode: maxima
+ * compile-command: "maxima --batch jacobian.mac"
+ * End:
+ */

dune/geometry/quadraturerules/jacobian.mac

@@ -34,15 +34,15 @@ weight(n,a,b,xx) := A(n,a,b)/A(n-1,a,b)*gamma_n(n-1,a,b)/(at(ratexpand(jacobi_p(
 write_quad("jacobi_1_0_imp.hh", "Jacobi1", "jacobian.mac",
 		   lambda([i],jacobi(i,1,0)),
            lambda([i,p],weight(i,1,0,p)),
-		   maxp);
+		   maxp)$
 write_quad("jacobi_2_0_imp.hh", "Jacobi2", "jacobian.mac",
 		   lambda([i],jacobi(i,2,0)),
            lambda([i,p],weight(i,2,0,p)),
-           maxp);
+           maxp)$
 write_quad("gauss_imp.hh", "Gauss", "jacobian.mac",
 		   lambda([i],jacobi(i,0,0)),
            lambda([i,p],weight(i,0,0,p)),
-           maxp);
+           maxp)$
 
 /*
  * Local variables:

dune/geometry/quadraturerules/writequad.mac

+load(to_poly_solve)$
+
+epsilon: 1e-110$
+
 write_quad(file, name, scriptname, polynome_fnkt, weight_fnkt, maxp) :=
 block([i,N,n_,points,p,w,fd],
 fd: openw(file),
-printf(fd, "// -*-c++-*-
+printf(fd, "// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+//
 // WARNING
 // This file is automatically generated by ~a! Don't edit by hand!
 
@@ -25,8 +31,8 @@ namespace Dune {
 ", scriptname, name),
 for i: 1 thru maxp/2+1 do block([],
   printf(fd,"    // order ~d,~d
-    case ~d:
-    case ~d:
+    case ~d :
+    case ~d :
 ", 2*i-2,2*i-1, 2*i-2,2*i-1),
   points: makelist(rhs(p), p, bfallroots(polynome_fnkt(i))),
   weights: makelist(bfloat(weight_fnkt(i, p)), p, points),
@@ -40,19 +46,21 @@ for i: 1 thru maxp/2+1 do block([],
 ",N),
   for n_: 1 thru N do block([n],
     n: S[n_],
-    p: points[n],
+    p: (points[n]+1)/2,
+    p: max(p-epsilon,-epsilon)+epsilon,
+    w: weights[n]/2,
     printf(fd,"      _points[~d] = ~h;
-",n_-1,(p+1)/2),
+",n_-1, p),
     printf(fd,"      _weight[~d] = ~h;
-",n_-1, weights[n]/2)
+",n_-1, w)
   ),
   printf(fd,"      break;
 
 ")
 ),
 
-printf(fd,"    default:
-	      DUNE_THROW(QuadratureOrderOutOfRange, \"Quadrature rule \" << p << \" not supported!\");
+printf(fd,"    default :
+      DUNE_THROW(QuadratureOrderOutOfRange, \"Quadrature rule \" << p << \" not supported!\");
     }
   }
 
@@ -69,8 +77,8 @@ printf(fd,"    default:
 
 for i: 1 thru maxp/2+1 do block([],
   printf(fd,"    // order ~d,~d
-    case ~d:
-    case ~d:
+    case ~d :
+    case ~d :
 ", 2*i-2,2*i-1, 2*i-2,2*i-1),
   points: makelist(rhs(p), p, bfallroots(polynome_fnkt(i))),
   weights: makelist(bfloat(weight_fnkt(i,p)), p, points),
@@ -84,18 +92,20 @@ for i: 1 thru maxp/2+1 do block([],
 ",N),
   for n_: 1 thru N do block([n],
     n: S[n_],
-    p: points[n],
+    p: (points[n]+1)/2,
+    p: max(p-epsilon,-epsilon)+epsilon,
+    w: weights[n]/2,
     printf(fd,"      _points[~d] = \"~h\";
-",n_-1,(p+1)/2),
+",n_-1,p),
     printf(fd,"      _weight[~d] = \"~h\";
-",n_-1,weights[n]/2)
+",n_-1,w)
   ),
   printf(fd,"      break;
 
 ")
 ),
 
-printf(fd,"    default:
+printf(fd,"    default :
       DUNE_THROW(QuadratureOrderOutOfRange, \"Quadrature rule \" << p << \" not supported!\");
     }
   }