added orthonormal bases for triangle, quadrilateral, tetrahedron,

pyramid, prism and hexahedron up to order of 8

[[Imported from SVN: r1388]]

fem/dgspace.hh

@@ -7,6 +7,7 @@
 
 #include <dune/fem/common/discretefunctionspace.hh>
 #include <dune/fem/dgspace/monomialbase.hh>
+#include <dune/fem/dgspace/orthonormalbase.hh>
 #include <dune/fem/dgspace/dgmapper.hh>
 #include <dune/fem/dofmanager.hh>
 
@@ -17,14 +18,16 @@ namespace Dune {
   //!  DiscreteFunctionSpace for discontinous functions
   //
   //**********************************************************************
-  template< class FunctionSpaceType, class GridType, int polOrd , class
-      DofManagerType = DofManager<GridType,DefaultGridIndexSet<GridType,LevelIndex> > >
+  template< class FunctionSpaceType, class GridType, int polOrd ,
+      template<class> class BaseFunctionSet,
+      class DofManagerType =
+        DofManager<GridType,DefaultGridIndexSet<GridType,LevelIndex> > >
   class DGDiscreteFunctionSpace
     : public DiscreteFunctionSpaceInterface
       <  FunctionSpaceType , GridType,
           DGDiscreteFunctionSpace < FunctionSpaceType , GridType,
-              polOrd, DofManagerType >,
-          MonomialBaseFunctionSet < FunctionSpaceType > >
+              polOrd, BaseFunctionSet, DofManagerType >,
+          BaseFunctionSet < FunctionSpaceType > >
   {
     enum { DGFSpaceId = 123456789 };
 
@@ -43,15 +46,15 @@ namespace Dune {
     typedef typename DofManagerType::MemObjectType MemObjectType;
 
     // BaseFunctionSet we are using
-    typedef MonomialBaseFunctionSet<FunctionSpaceType> BaseFunctionSetType;
+    typedef BaseFunctionSet<FunctionSpaceType> BaseFunctionSetType;
 
     typedef DGDiscreteFunctionSpace
-    < FunctionSpaceType , GridType , polOrd , DofManagerType >
+    < FunctionSpaceType , GridType , polOrd , BaseFunctionSet, DofManagerType >
     DGDiscreteFunctionSpaceType;
 
     typedef DiscreteFunctionSpaceInterface
     <FunctionSpaceType , GridType, DGDiscreteFunctionSpaceType,
-        MonomialBaseFunctionSet <FunctionSpaceType> >
+        BaseFunctionSet <FunctionSpaceType> >
     DiscreteFunctionSpaceType;
 
     /** \todo Please doc me! */
@@ -90,6 +93,7 @@ namespace Dune {
     const BaseFunctionSetType &
     getBaseFunctionSet ( EntityType &en ) const
     {
+      setType(base_,en.geometry().type());
       return base_;
     }
 
@@ -127,7 +131,14 @@ namespace Dune {
     };
 
   private:
-    BaseFunctionSetType base_;
+    void setType(MonomialBaseFunctionSet<FunctionSpaceType> & b,
+                 ElementType t) const {};
+    void setType(OrthonormalBaseFunctionSet<FunctionSpaceType> & b,
+                 ElementType t) const
+    {
+      b.changeType(t);
+    };
+    mutable BaseFunctionSetType base_;
     // mapper for function space
     DGMapperType mapper_;
   };

fem/dgspace/monomialbase.hh

 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 // vi: set et ts=4 sw=2 sts=2:
-#ifndef __DUNE__MONOMIALBASE_HH__
-#define __DUNE__MONOMIALBASE_HH__
+#ifndef DUNE__MONOMIALBASE_HH
+#define DUNE__MONOMIALBASE_HH
 
 #include <dune/common/matvec.hh>
 #include <dune/common/simplevector.hh>
@@ -129,4 +129,4 @@ namespace Dune {
 
 } // end namespace Dune
 
-#endif
+#endif // DUNE__MONOMIALBASE_HH

fem/dgspace/orthonormalbase.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <string>
+
+#include <dune/fem/dgspace/orthonormalbase.hh>
+
+template <class FunctionSpaceType>
+OrthonormalBaseFunctionSet<FunctionSpaceType >::
+OrthonormalBaseFunctionSet( FunctionSpaceType & fuspace,
+                            int polOrder )
+  :
+    BaseFunctionSetDefault
+    <FunctionSpaceType,OrthonormalBaseFunctionSet <FunctionSpaceType > >
+      ( fuspace ),
+    polOrder_(polOrder),
+    type_(Dune::line)
+{
+  assert(DimRange == 1);
+  switch (DimDomain)
+  {
+  case 2 :
+    numOfBaseFct_ = (polOrder + 2) * (polOrder + 1) / 2;
+    break;
+  case 3 :
+    numOfBaseFct_ = ((polOrder+1)*(polOrder+2)*(2*polOrder+3)/6 +
+                     (polOrder+1)*(polOrder+2)/2)/2;
+    break;
+  default :
+    DUNE_THROW(NotImplemented, "OrthonormalBaseFunctionSet only supports 2D and 3D Domain");
+  }
+}
+
+template <class FunctionSpaceType>
+void OrthonormalBaseFunctionSet<FunctionSpaceType >::
+real_evaluate( int baseFunct, const FieldVector<deriType, 0> &diffVariable,
+               const Domain & x,  Range & phi ) const
+{
+  switch(type_)
+  {
+  case triangle :
+    phi[0] = eval_triangle_2d (baseFunct, x); break;
+  case quadrilateral :
+    phi[0] = eval_quadrilateral_2d (baseFunct, x); break;
+  case tetrahedron :
+    phi[0] = eval_tetrahedron_3d (baseFunct, x); break;
+  case pyramid :
+    phi[0] = eval_pyramid_3d (baseFunct, x); break;
+  case prism :
+    phi[0] = eval_prism_3d (baseFunct, x); break;
+  case hexahedron :
+    phi[0] = eval_hexahedron_3d (baseFunct, x); break;
+  default :
+    DUNE_THROW(NotImplemented, "ElementType not suppoerted by OrthonormalBaseFunctionSet");
+  }
+}
+
+template <class FunctionSpaceType>
+void OrthonormalBaseFunctionSet<FunctionSpaceType >::
+jacobian ( int baseFunct, const Domain & x, JacobianRange & grad ) const
+{
+  switch(type_)
+  {
+  case triangle :
+    return grad_triangle_2d (baseFunct, x, grad);
+  case quadrilateral :
+    return grad_quadrilateral_2d (baseFunct, x, grad);
+  case tetrahedron :
+    return grad_tetrahedron_3d (baseFunct, x, grad);
+  case pyramid :
+    return grad_pyramid_3d (baseFunct, x, grad);
+  case prism :
+    return grad_prism_3d (baseFunct, x, grad);
+  case hexahedron :
+    return grad_hexahedron_3d (baseFunct, x, grad);
+  default :
+    DUNE_THROW(NotImplemented, "ElementType not suppoerted by OrthonormalBaseFunctionSet");
+  }
+}
+
+template <class FunctionSpaceType>
+void OrthonormalBaseFunctionSet<FunctionSpaceType >::
+real_evaluate( int baseFunct, const FieldVector<deriType, 1> &diffVariable,
+               const Domain & x,  Range & phi ) const
+{
+  JacobianRange grad;
+  jacobian(baseFunct, x, grad);
+  phi[0] = grad[diffVariable[0]];
+}
+
+template <class FunctionSpaceType> template <int diffOrd>
+void OrthonormalBaseFunctionSet<FunctionSpaceType >::
+evaluate( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
+          const Domain & x,  Range & phi ) const
+{
+  assert(baseFunct < numOfBaseFct_);
+  real_evaluate(baseFunct, diffVariable, x, phi);
+}
+
+template <class FunctionSpaceType> template <int diffOrd, class QuadratureType>
+void OrthonormalBaseFunctionSet<FunctionSpaceType >::
+evaluate( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
+          QuadratureType & quad, int quadPoint, Range & phi ) const
+{
+  Domain x = quad.point(quadPoint);
+  real_evaluate(baseFunct, diffVariable, x, phi);
+}

fem/dgspace/orthonormalbase.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE__ORTHONORMALBASE_HH
+#define DUNE__ORTHONORMALBASE_HH
+
+#include <dune/common/matvec.hh>
+#include <dune/common/simplevector.hh>
+#include <dune/fem/common/basefunctions.hh>
+
+namespace Dune {
+
+  /** @defgroup OrthonormalBaseFunctionSet The OrthonormalBaseFunctionSet
+     @ingroup BaseFunctionSet
+     The OrthonormalBaseFunctionSet is an special implementation of the interface
+     BaseFunctionSetInterface. This class offers orthonormal functions as a base.
+     Currently only functionspaces mapping from 2D and 3D real to real.
+
+     @{
+   */
+
+  //*************************************************************************
+  //
+  //  --OrthonormalBaseFunctionSet
+  //
+  //*************************************************************************
+
+  //! \todo Please doc me!
+  template<class FunctionSpaceType>
+  class OrthonormalBaseFunctionSet
+    : public BaseFunctionSetDefault
+      <FunctionSpaceType, OrthonormalBaseFunctionSet<FunctionSpaceType> >
+  {
+    typedef typename FunctionSpaceType::Domain Domain;
+    typedef typename FunctionSpaceType::DomainField DomainField;
+    typedef typename FunctionSpaceType::Range Range;
+    typedef typename FunctionSpaceType::JacobianRange JacobianRange;
+    enum { DimDomain = FunctionSpaceType::DimDomain };
+    enum { DimRange  = FunctionSpaceType::DimRange  };
+  public:
+
+    //! Constructor
+    OrthonormalBaseFunctionSet (FunctionSpaceType & fuspace, int polOrder);
+
+    //! change type
+    void changeType(ElementType & type)
+    {
+      type_ = type;
+    }
+
+    //! return the number of base functions for this BaseFunctionSet
+    int getNumberOfBaseFunctions () const
+    {
+      return numOfBaseFct_;
+    };
+
+    //! evaluate base function baseFunct with the given diffVariable and a
+    //! point x and range phi
+    template <int diffOrd>
+    inline
+    void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
+                    const Domain & x,  Range & phi ) const;
+
+    //! evaluate base fucntion baseFunct at a given quadrature point
+    //! the identifier of the quadrature is stored to check , whether the
+    //! qaudrature has changed an the values at the quadrature have to be
+    //! calulated again
+    template <int diffOrd, class QuadratureType>
+    inline
+    void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
+                    QuadratureType & quad, int quadPoint, Range & phi ) const;
+
+    //! evaluate the gradient
+    void jacobian ( int baseFunct, const Domain & x, JacobianRange & grad ) const;
+
+    //! \todo Please doc me!
+    void print (std::ostream& s, int baseFunct) const {
+      s << "print is missing\n";
+    };
+    //! \todo Please doc me!
+    void print (std::ostream& s) const {
+      s << "[";
+      for (int i = 0; i < numOfBaseFct_ - 1; i++) {
+        s << "["; print(s, i); s << "], ";
+      }
+      s << "["; print(s, numOfBaseFct_); s << "]]";
+    }
+
+  private:
+    //! polynomial order of the base
+    int polOrder_;
+
+    //! number of different basefunctions
+    int numOfBaseFct_;
+
+    //! Elementtype for which this basis is orthonormal
+    ElementType type_;
+
+    //! eval
+    inline void real_evaluate ( int baseFunct,
+                                const FieldVector<deriType, 0> &diffVariable,
+                                const Domain & x,
+                                Range & phi ) const;
+
+    //! grad
+    inline void real_evaluate ( int baseFunct,
+                                const FieldVector<deriType, 1> &diffVariable,
+                                const Domain & x,
+                                Range & phi ) const;
+
+    double eval_triangle_2d ( int i, const Domain & xi ) const;
+    double eval_quadrilateral_2d ( int i, const Domain & xi ) const;
+    double eval_tetrahedron_3d ( int i, const Domain & xi ) const;
+    double eval_pyramid_3d ( int i, const Domain & xi ) const;
+    double eval_prism_3d ( int i, const Domain & xi ) const;
+    double eval_hexahedron_3d ( int i, const Domain & xi ) const;
+
+    void grad_triangle_2d ( int i, const Domain & xi,
+                            JacobianRange & grad ) const;
+    void grad_quadrilateral_2d ( int i, const Domain & xi,
+                                 JacobianRange & grad ) const;
+    void grad_tetrahedron_3d ( int i, const Domain & xi,
+                               JacobianRange & grad ) const;
+    void grad_pyramid_3d ( int i, const Domain & xi,
+                           JacobianRange & grad ) const;
+    void grad_prism_3d ( int i, const Domain & xi,
+                         JacobianRange & grad ) const;
+    void grad_hexahedron_3d ( int i, const Domain & xi,
+                              JacobianRange & grad ) const;
+
+  }; // end class OrthonormalBaseFunctionSet
+
+  //! overloading the out stream for printing of monomial bases
+  template< class T >
+  inline std::ostream&
+  operator<< (std::ostream& s, OrthonormalBaseFunctionSet<T> & m)
+  {
+    m.print(s);
+    return s;
+  }
+
+  /** @} end documentation group */
+
+#include "orthonormalbase.cc"
+#include "orthonormalbase_impl.cc"
+
+} // end namespace Dune
+
+#endif // DUNE_ORTHONORMALBASE_HH