Allow to measure vector-valued functions

Unfortunately this involves moving the method bodies inside the class
declaration (I couldn't get the template code to compile otherwise).
I tried to get a patch with as little noise as possible.  Unfortunately
this means that indentation is now off in this file.

dune/fufem/discretizationerror.hh

@@ -19,24 +19,11 @@ public:
 
     /** \brief Compute L2 error between a grid function and an arbitrary function
      */
-    static double computeL2Error(const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >* a,
-               const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,1> >* b,
-               QuadratureRuleKey quadKey);
-
-    static double computeH1HalfNormDifferenceSquared(const GridView& gridView,
-                                                     const VirtualDifferentiableFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,1> >* u,
-                                                     const VirtualDifferentiableFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,1> >* v,
-                                                     QuadratureRuleKey quadKey);
-
-
-};
-
-template <class GridView>
-double DiscretizationError<GridView>::
-computeL2Error(const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >* a,
-               const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,1> >* b,
+    template <int blocksize>
+    static double computeL2Error(const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >* a,
+               const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,blocksize> >* b,
                QuadratureRuleKey quadKey)
-{
+    {
 
     // The error to be computed
     double error = 0;
@@ -55,13 +42,13 @@ computeL2Error(const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1
         for (size_t i=0; i<quad.size(); i++) {
 
             // Evaluate function a
-            Dune::FieldVector<double,1> aValue;
+            Dune::FieldVector<double,blocksize> aValue;
             a->evaluateLocal(*eIt, quad[i].position(),aValue);
 
             // Evaluate function b.  If it is a grid function use that to speed up the evaluation
-            Dune::FieldVector<double,1> bValue;
-            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >*>(b))
-                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >*>(b)->evaluateLocal(*eIt,
+            Dune::FieldVector<double,blocksize> bValue;
+            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >*>(b))
+                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >*>(b)->evaluateLocal(*eIt,
                                                                                                                       quad[i].position(),
                                                                                                                       bValue
                                                                                                                      );
@@ -77,15 +64,15 @@ computeL2Error(const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1
 
     return std::sqrt(error);
 
-}
+    }
 
-template <class GridView>
-double DiscretizationError<GridView>::
-computeH1HalfNormDifferenceSquared(const GridView& gridView,
-                                   const VirtualDifferentiableFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,1> >* u,
-                                   const VirtualDifferentiableFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,1> >* v,
-                                   QuadratureRuleKey quadKey)
-{
+
+    template <int blocksize>
+    static double computeH1HalfNormDifferenceSquared(const GridView& gridView,
+                                                     const VirtualDifferentiableFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,blocksize> >* u,
+                                                     const VirtualDifferentiableFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,blocksize> >* v,
+                                                     QuadratureRuleKey quadKey)
+    {
     double norm = 0.0;
     typedef typename GridView::template Codim<0>::Iterator ElementIterator;
 
@@ -106,16 +93,16 @@ computeH1HalfNormDifferenceSquared(const GridView& gridView,
             const double integrationElement = eIt->geometry().integrationElement(quadPos);
 
             // Evaluate integral
-            Dune::FieldMatrix<double,1,dim> u_di;
-            Dune::FieldMatrix<double,1,dim> v_di;
+            Dune::FieldMatrix<double,blocksize,dim> u_di;
+            Dune::FieldMatrix<double,blocksize,dim> v_di;
 
-            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >*>(u))
-                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >*>(u)->evaluateDerivativeLocal(*eIt, quadPos, u_di);
+            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >*>(u))
+                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >*>(u)->evaluateDerivativeLocal(*eIt, quadPos, u_di);
             else
                 u->evaluateDerivative(eIt->geometry().global(quadPos), u_di);
 
-            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >*>(v))
-                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,1> >*>(v)->evaluateDerivativeLocal(*eIt, quadPos, v_di);
+            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >*>(v))
+                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,blocksize> >*>(v)->evaluateDerivativeLocal(*eIt, quadPos, v_di);
             else
                 v->evaluateDerivative(eIt->geometry().global(quadPos), v_di);
 
@@ -124,8 +111,10 @@ computeH1HalfNormDifferenceSquared(const GridView& gridView,
         }
     }
     return norm;
-}
+    }
 
 
+};
+
 #endif