Now also face quadratures can be cached

[[Imported from SVN: r1655]]

fem/common/cachingbase.cc

@@ -7,7 +7,7 @@
 template <class FunctionSpaceType>
 CachingBaseFunctionSet<FunctionSpaceType>::
 CachingBaseFunctionSet(FunctionSpaceType& fuspace,
-                       ElementType elType,
+                       GeometryType elType,
                        int nBaseFnc) :
   BaseFunctionSetDefault<FunctionSpaceType,
       CachingBaseFunctionSet<FunctionSpaceType> > (fuspace),
@@ -88,11 +88,23 @@ gradients(int baseFunct, const QuadratureType& quad) const {
   return it->second[baseFunct];
 }
 
+template <class FunctionSpaceType>
+template <class QuadratureType>
+const std::vector<FunctionSpaceType::Range>&
+CachingBaseFunctionSet<FunctionSpaceType>::
+faces(int faceIdx, int baseFct, const QuadratureType& quad) const {
+  ConstFaceMapIterator it = faces_.find(quad.getIdentifier());
+
+  assert(it != faces_.end());
+  assert(baseFct < numOfBaseFct_);
+
+  return it->second[faceIdx][baseFct];
+}
 
 //! if the quadrature is new, then once cache the values
 template <class FunctionSpaceType>
 template <class QuadratureType>
-void CachingBaseFunctionSet<FunctionSpaceType >::
+void CachingBaseFunctionSet<FunctionSpaceType>::
 registerQuadrature( const QuadratureType &quad )
 {
   // check if quadrature is already registered
@@ -107,8 +119,9 @@ registerQuadrature( const QuadratureType &quad )
     // Get iterators to the right map entries
     std::pair<RangeMapIterator, bool> rp =
       vals_.insert(std::make_pair(identifier,
-                                  std::vector<std::vector<Range> >(nBaseFct,
-                                                                   std::vector<Range>(nQuadPts, tmp))));
+                                  std::vector<std::vector<Range> >
+                                    (nBaseFct,
+                                    std::vector<Range>(nQuadPts, tmp))));
 
     std::pair<JacobianMapIterator, bool> jp =
       grads_.insert(std::make_pair(identifier,
@@ -127,13 +140,61 @@ registerQuadrature( const QuadratureType &quad )
   } // end if
 }
 
+template <class FunctionSpaceType>
+template <class QuadratureType, class EntityType>
+void CachingBaseFunctionSet<FunctionSpaceType>::
+registerQuadrature(const QuadratureType &quad, const EntityType& en) {
+  //- Local typedefs
+  typedef typename EntityType::IntersectionIterator IntersectionIterator;
+
+  //- Actual code
+  // Check if quadrature is already present
+  int identifier = quad.getIdentifier();
+  if (faces_.find(identifier) == faces_.end()) {
+
+    // Initialise face values
+    int nBaseFct = getNumberOfBaseFunctions();
+    int nQuadPts = quad.nop();
+    Range tmp(0.0);
+
+    // Count faces
+    int count = 0;
+    IntersectionIterator endit = en.iend();
+    for (IntersectionIterator it = en.ibegin(); it != endit; ++it, ++count) ;
+
+    // Get iterator to the right map entry
+    std::pair<FaceMapIterator, bool> rp =
+      faces_.
+      insert(std::make_pair(identifier,
+                            std::vector<std::vector<std::vector<Range> > >
+                              (count,
+                              std::vector<std::vector<Range> >
+                                (nBaseFct,
+                                std::vector<Range>(nQuadPts, tmp)))));
+
+    for (IntersectionIterator it = en.ibegin(); it != endit; ++it) {
+      for (int i = 0; i < nBaseFct; ++i) {
+        for (int j = 0; j < nQuadPts; ++j) {
+          // * Is this correct? (intersectionSelfLocal().global(quad.point(j)))
+          this->eval(i,
+                     it.intersectionSelfLocal().global(quad.point(j)),
+                     rp.first->second[it.numberInSelf()][i][j]);
+        }
+      }
+    } // end for
+
+  } // end if
+}
+
+
 template <class FunctionSpaceType>
 const FunctionSpaceType::Range
 CachingBaseFunctionSet<FunctionSpaceType>::
 extractGradientComp(const JacobianRange& jr, int idx) const {
   Range result;
   for (int i = 0; i < DimRange; ++i) {
-    result[i] = jr(idx, i);
+    // * check order of operators
+    result[i] = jr[idx][i];
   }
   return result;
 }

fem/common/cachingbase.hh

 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 // vi: set et ts=4 sw=2 sts=2:
-#ifndef __DUNE_CACHINGBASE_HH__
-#define __DUNE_CACHINGBASE_HH__
+#ifndef DUNE_CACHINGBASE_HH
+#define DUNE_CACHINGBASE_HH
 
 #include <map>
 
@@ -95,6 +95,12 @@ namespace Dune {
     const std::vector<JacobianRange>& gradients(int baseFunct,
                                                 const QuadratureType& quad) const;
 
+    //! Alternative acces to precomputed base function values on the faces
+    template <class QuadratureType>
+    const std::vector<Range>& faces(int faceIndex,
+                                    int baseFunct,
+                                    const QuadratureType& quad) const;
+
     //! get a reference of the base function baseFunct
     //! this is the same concept for all basis, we have a number of
     //! base functions, and store the pointers in a vector
@@ -110,9 +116,19 @@ namespace Dune {
     template <class QuadratureType>
     void registerQuadrature(const QuadratureType & quad);
 
+    //! Register face quadrature
+    template <class QuadratureType, class EntityType>
+    void registerQuadrature(const QuadratureType& quad,
+                            const EntityType& en);
+
   private:
     //- Local typedefs
-    typedef typename std::map<IdentifierType, std::vector<std::vector<Range> > > RangeMap;
+    typedef typename std::map<IdentifierType,
+        std::vector<std::vector<std::vector<Range> > > > FaceMap;
+    typedef typename FaceMap::iterator FaceMapIterator;
+    typedef typename FaceMap::const_iterator ConstFaceMapIterator;
+    typedef typename std::map<IdentifierType,
+        std::vector<std::vector<Range> > > RangeMap;
     typedef typename RangeMap::iterator RangeMapIterator;
     typedef typename RangeMap::const_iterator ConstRangeMapIterator;
     typedef typename std::map<IdentifierType,
@@ -136,6 +152,9 @@ namespace Dune {
     //! map with cached values for base function gradients
     JacobianMap grads_;
 
+    //! map with cached values for face values of base functions
+    FaceMap faces_;
+
   }; // end class CachingBaseFunctionSet