Ich schlage hiermit vor, eine Sammlung von häufig gebrauchten Operatoren

in Dune anzulegen.  Den Anfang machen der Laplaceoperator und die
Massenmatrix.  Es ist der Code aus example/heatequation/fem.cc.  Dort
werden sie in kürze verschwinden.

[[Imported from SVN: r599]]

fem/operator/laplace.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_LAPLACE_CC__
+#define __DUNE_LAPLACE_CC__
+
+#include <dune/fem/feoperator.hh>
+#include <dune/fem/feop/spmatrix.hh>
+
+#include <dune/fem/fastquad.hh>
+
+
+namespace Dune
+{
+
+  /** \brief The Laplace operator
+   */
+  template <class DiscFunctionType, class TensorType>
+  class LaplaceFEOp :
+    public FiniteElementOperator<DiscFunctionType,
+        SparseRowMatrix<double>,
+        LaplaceFEOp<DiscFunctionType,TensorType> > {
+
+    typedef typename DiscFunctionType::FunctionSpaceType FunctionSpaceType;
+    typedef typename FunctionSpaceType::GridType GridType;
+    enum { dim = GridType::dimension };
+    typedef typename GridType::template Traits<0>::CoordType CoordType;
+    typedef typename FunctionSpaceType::JacobianRange JacobianRange;
+    typedef typename FunctionSpaceType::RangeField RangeFieldType;
+
+    typedef typename FiniteElementOperator<DiscFunctionType,
+        SparseRowMatrix<double>,
+        LaplaceFEOp<DiscFunctionType, TensorType> >::OpMode OpMode;
+
+    mutable JacobianRange grad;
+    mutable JacobianRange othGrad;
+    mutable JacobianRange mygrad[4];
+
+  public:
+    FastQuad < typename FunctionSpaceType::RangeField, typename
+        FunctionSpaceType::Domain , 1 > quad;
+
+    DiscFunctionType *stiffFunktion_;
+    TensorType *stiffTensor_;
+
+    LaplaceFEOp( const typename DiscFunctionType::FunctionSpace &f, OpMode opMode ) :    //= ASSEMBLED ) :
+                                                                                         FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,LaplaceFEOp<DiscFunctionType,TensorType> >( f, opMode ) ,
+                                                                                         quad ( *(f.getGrid().template lbegin<0> (0))), stiffFunktion_(NULL), stiffTensor_(NULL)
+    {}
+
+    LaplaceFEOp( const DiscFunctionType &stiff, const typename DiscFunctionType::FunctionSpace &f, OpMode opMode ) :    //= ASSEMBLED ) :
+                                                                                                                        FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,LaplaceFEOp<DiscFunctionType,TensorType> >( f, opMode ) ,
+                                                                                                                        quad ( *(f.getGrid().template lbegin<0> (0))), stiffFunktion_(&stiff), stiffTensor_(NULL)
+    {}
+
+    LaplaceFEOp( TensorType &stiff, const typename DiscFunctionType::FunctionSpace &f, OpMode opMode ) :    //= ASSEMBLED ) :
+                                                                                                            FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,LaplaceFEOp<DiscFunctionType,TensorType> >( f, opMode ) ,
+                                                                                                            quad ( *(f.getGrid().template lbegin<0> (0))), stiffFunktion_(NULL), stiffTensor_(&stiff)
+    {}
+
+    const SparseRowMatrix<double>* getMatrix() const {
+      assert(this->matrix_);
+      return this->matrix_;
+    }
+
+    SparseRowMatrix<double>* newEmptyMatrix( ) const
+    {
+      return new SparseRowMatrix<double>( this->functionSpace_.size ( this->functionSpace_.getGrid().maxlevel() ) ,
+                                          this->functionSpace_.size ( this->functionSpace_.getGrid().maxlevel() ) ,
+                                          15 * (dim-1) , 0.0 );
+    }
+
+    void prepareGlobal ( const DiscFunctionType &arg, DiscFunctionType &dest )
+    {
+      this->arg_  = &arg.argument();
+      this->dest_ = &dest.destination();
+      assert(this->arg_ != NULL); assert(this->dest_ != NULL);
+      dest.clear();
+    }
+
+    template <class EntityType>
+    double getLocalMatrixEntry( EntityType &entity, const int i, const int j ) const
+    {
+      enum { dim = GridType::dimension };
+      typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
+
+      const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
+
+      // calc Jacobian inverse before volume is evaluated
+      const Mat<dim,dim,double>& inv = entity.geometry().Jacobian_inverse(quad.point(0));
+
+      const double vol = entity.geometry().integration_element(quad.point(0));
+
+      double val = 0.;
+      for ( int pt=0; pt < quad.nop(); pt++ )
+      {
+        baseSet.jacobian(i,quad,pt,grad);
+
+        // multiply with transpose of jacobian inverse
+        grad(0) = inv.mult_t (grad(0));
+
+        if( i != j )
+        {
+          baseSet.jacobian(j,quad,pt,othGrad);
+
+          // multiply with transpose of jacobian inverse
+          othGrad(0) = inv.mult_t(othGrad(0));
+
+          val += ( grad(0) * othGrad(0) ) * quad.weight( pt );
+        }
+        else
+        {
+          val += ( grad(0) * grad(0) ) * quad.weight( pt );
+        }
+      }
+      val *= vol;
+      return val;
+    }
+
+
+    template < class EntityType, class MatrixType>
+    void getLocalMatrix( EntityType &entity, const int matSize, MatrixType& mat) const {
+      enum { dim = GridType::dimension };
+      typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
+
+      const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
+
+      // calc Jacobian inverse before volume is evaluated
+      const Mat<dim,dim,double>& inv = entity.geometry().Jacobian_inverse(quad.point(0));
+
+      const double vol = entity.geometry().integration_element(quad.point(0));
+      int i,j;
+
+      for(i=0; i<matSize; i++)
+        for (j=0; j<=i; j++ )
+          mat(i,j)=0.0;
+
+      for ( int pt=0; pt < quad.nop(); pt++ )
+      {
+        for(i=0; i<matSize; i++)
+        {
+          baseSet.jacobian(i,quad,pt,mygrad[i]);
+
+          // multiply with transpose of jacobian inverse
+          mygrad[i](0) = inv.mult_t (mygrad[i](0));
+        }
+
+        typename FunctionSpaceType::Range ret;
+
+        if(stiffTensor_) {
+          stiffTensor_->evaluate(entity.geometry().global(quad.point(pt)),ret);
+          ret(0) *= quad.weight( pt );
+          for(i=0; i<matSize; i++)
+            for (j=0; j<=i; j++ )
+              mat(i,j) += ( mygrad[i](0) * mygrad[j](0) ) * ret(0);
+        }
+        else{
+          for(i=0; i<matSize; i++)
+            for (j=0; j<=i; j++ )
+              mat(i,j) += ( mygrad[i](0) * mygrad[j](0) ) * quad.weight( pt );
+        }
+
+
+
+      }
+
+      for(i=0; i<matSize; i++)
+        for (j=0; j<=i; j++ )
+          mat(i,j) *= vol;
+
+      for(i=0; i<matSize; i++)
+        for (j=matSize; j>i; j--)
+          mat(i,j) = mat(j,i);
+
+      return;
+    }
+
+  protected:
+  };   // end class
+
+} // end namespace
+
+#endif

fem/operator/massmatrix.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_MASSMATRIX_HH__
+#define __DUNE_MASSMATRIX_HH__
+
+namespace Dune {
+
+  /** \brief The mass matrix
+   */
+  template <class DiscFunctionType>
+  class MassMatrixFEOp :
+    public FiniteElementOperator<DiscFunctionType,
+        SparseRowMatrix<double>,
+        MassMatrixFEOp<DiscFunctionType> > {
+    typedef typename DiscFunctionType::FunctionSpaceType FunctionSpaceType;
+    typedef typename FiniteElementOperator<DiscFunctionType,
+        SparseRowMatrix<double>,
+        MassMatrixFEOp<DiscFunctionType> >::OpMode OpMode;
+
+    typedef typename FunctionSpaceType::RangeField RangeFieldType;
+    typedef typename FunctionSpaceType::GridType GridType;
+    typedef typename FunctionSpaceType::Range RangeType;
+
+  public:
+    FastQuad < typename FunctionSpaceType::RangeField, typename
+        FunctionSpaceType::Domain , 1 > quad;
+
+    MassMatrixFEOp( const typename DiscFunctionType::FunctionSpace &f, OpMode opMode ) :   //= ON_THE_FLY ) :
+                                                                                           FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,MassMatrixFEOp<DiscFunctionType> >( f, opMode ) ,
+                                                                                           quad ( *(f.getGrid().template lbegin<0> (0))) {}
+
+    SparseRowMatrix<double>* newEmptyMatrix( ) const {
+      return new SparseRowMatrix<double>( this->functionSpace_.size ( this->functionSpace_.getGrid().maxlevel() ) ,
+                                          this->functionSpace_.size ( this->functionSpace_.getGrid().maxlevel() ) ,
+                                          10, 0.0 );
+    }
+
+    template <class EntityType>
+    double getLocalMatrixEntry( EntityType &entity, const int i, const int j ) const
+    {
+      typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
+      double val = 0.;
+
+      Vec<GridType::dimension> tmp(1.0);
+
+      const BaseFunctionSetType & baseSet
+        = this->functionSpace_.getBaseFunctionSet( entity );
+
+      RangeType v1 (0.0);
+      RangeType v2 (0.0);
+
+      const double vol =
+        entity.geometry().integration_element( tmp );
+      for ( int pt=0; pt < quad.nop(); pt++ )
+      {
+        baseSet.eval(i,quad,pt,v1);
+        baseSet.eval(j,quad,pt,v2);
+        val += ( v1 * v2 ) * quad.weight( pt );
+      }
+      val *= vol;
+      return val;
+    }
+
+    template < class EntityType, class MatrixType>
+    void getLocalMatrix( EntityType &entity, const int matSize, MatrixType& mat) const
+    {
+      int i,j;
+      typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
+      const BaseFunctionSetType & baseSet
+        = this->functionSpace_.getBaseFunctionSet( entity );
+
+      static Vec<GridType::dimension> tmp(1.0);
+      const double vol = entity.geometry().integration_element(tmp);
+
+      static RangeType v[4];
+
+      for(i=0; i<matSize; i++)
+        for (j=0; j<=i; j++ )
+          mat(i,j)=0.0;
+
+      for ( int pt=0; pt < quad.nop(); pt++ )
+      {
+        for(i=0; i<matSize; i++)
+          baseSet.eval(i,quad,pt,v[i]);
+
+        for(i=0; i<matSize; i++)
+          for (j=0; j<=i; j++ )
+            mat(i,j) += ( v[i] * v[j] ) * quad.weight( pt );
+      }
+
+      for(i=0; i<matSize; i++)
+        for (j=0; j<=i; j++ )
+          mat(i,j) *= vol;
+
+      for(i=0; i<matSize; i++)
+        for (j=matSize; j>i; j--)
+          mat(i,j) = mat(j,i);
+
+
+      return;
+    }
+  protected:
+    DiscFunctionType *_tmp;
+  };
+
+} // namespace Dune
+
+#endif