we still need this things.

[[Imported from SVN: r2804]]

quadrature/fixedorder/quadtetratri.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_QUADTETRATRI_HH__
+#define __DUNE_QUADTETRATRI_HH__
+
+// the UG Quadratures
+#include "ugquadratures.hh"
+
+namespace Dune {
+
+  //******************************************************************
+  //
+  //!  Memorization of the number of quadrature points
+  //
+  //******************************************************************
+
+  static const double referenceVol_triangle    = 0.5;
+  static const double referenceVol_tetrahedron = 1.0/6.0;
+
+  //! specialization triangles
+  template <class Domain, class RangeField, int polOrd>
+  struct QuadraturePoints<Domain,RangeField,triangle,polOrd>
+  {
+    enum { dim = 2 };
+    enum { numberOfCorners = dim+1 };
+    typedef UG_Quadratures::QUADRATURE QUADRATURE;
+    static int numberOfQuadPoints ();
+    static int order ();
+    static Domain getPoint (int i);
+    static RangeField getWeight (int i);
+  };
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField, int polOrd>
+  int QuadraturePoints<Domain,RangeField,triangle,polOrd>::
+  numberOfQuadPoints()
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    return quad->nip;
+  }
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField, int polOrd>
+  int QuadraturePoints<Domain,RangeField,triangle,polOrd>::
+  order()
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    return quad->order;
+  }
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField,int polOrd>
+  RangeField QuadraturePoints<Domain,RangeField,triangle,polOrd>::
+  getWeight(int i)
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    RangeField ref(referenceVol_triangle);
+    RangeField
+    w(ref * static_cast<RangeField> (quad->weight[i]));
+    return w;
+  }
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField,int polOrd>
+  Domain QuadraturePoints<Domain,RangeField,triangle,polOrd>::
+  getPoint(int i)
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    Domain tmp;
+    for(int j=0; j<dim; j++)
+      tmp[j] = quad->local[i][j];
+    return tmp;
+  }
+
+  //*************************************************************
+  //! specialization tetrahedrons
+  template <class Domain, class RangeField, int polOrd>
+  struct QuadraturePoints<Domain,RangeField,tetrahedron,polOrd>
+  {
+    enum { dim = 3 }; // tetrahedrons
+    enum { numberOfCorners = dim+1 };
+    typedef UG_Quadratures::QUADRATURE QUADRATURE;
+    static int numberOfQuadPoints ();
+    static int order ();
+    static Domain getPoint (int i);
+    static RangeField getWeight (int i);
+  };
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField, int polOrd>
+  int QuadraturePoints<Domain,RangeField,tetrahedron,polOrd>::
+  numberOfQuadPoints()
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    return quad->nip;
+  }
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField, int polOrd>
+  int QuadraturePoints<Domain,RangeField,tetrahedron,polOrd>::
+  order()
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    return quad->order;
+  }
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField, int polOrd>
+  RangeField QuadraturePoints<Domain,RangeField,tetrahedron,polOrd>::
+  getWeight(int i)
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    RangeField ref(referenceVol_tetrahedron);
+    RangeField
+    w(ref * static_cast<RangeField> (quad->weight[i]));
+    return w;
+  }
+
+  //! the weight is the volume of the reference element
+  template <class Domain, class RangeField,int polOrd>
+  Domain QuadraturePoints<Domain,RangeField,tetrahedron,polOrd>::
+  getPoint(int i)
+  {
+    QUADRATURE * quad = UG_Quadratures::GetQuadratureRule(dim,numberOfCorners,polOrd);
+    Domain tmp;
+    for(int j=0; j<dim; j++)
+      tmp[j] = quad->local[i][j];
+    return tmp;
+  }
+
+} // end namespace Dune
+
+#endif

quadrature/fixedorder/ugquadratures.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+/** \file
+ * \brief Contains quadrature formulas taken form UG
+ */
+#ifndef DUNE_UG_QUADRATURES_HH
+#define DUNE_UG_QUADRATURES_HH
+
+namespace Dune {
+
+  namespace UG_Quadratures {
+
+    enum { dimension = 3 };
+
+    typedef double DOUBLE_VECTOR[dimension];
+    typedef double DOUBLE_VECTOR_3D[3];
+
+    /****************************************************************************/
+    /*                                                                          */
+    /* data structures exported by the corresponding source file                */
+    /*                                                                          */
+    /****************************************************************************/
+
+    /** \brief %Quadrature formula borrowed from UG */
+    struct QUADRATURE {
+
+      /** \brief Number of integration points     */
+      int nip;
+
+      /** \brief Order of quadrature rule         */
+      int order;
+
+      /** \brief %Array[nip] for local coordinates */
+      const DOUBLE_VECTOR_3D *local;
+
+      /** \brief %Array[nip] for weights           */
+      const double *weight;
+    };
+
+    struct GAUSS_POINT {
+      DOUBLE_VECTOR local;
+      DOUBLE_VECTOR global;
+      double weight;
+      DOUBLE_VECTOR Jinv[dimension];
+    };
+
+    /****************************************************************************/
+    /*                                                                          */
+    /* function declarations                                                    */
+    /*                                                                          */
+    /****************************************************************************/
+
+    /****************************************************************************/
+    /** \brief Provide a quadrature formula
+
+       \param dim - dimension of the quadrature domain
+       \param n - number of corners of the element
+       \param order - order of approximation
+
+       This function returns a pointer to a quadrature formula.
+
+       \section Example
+       \verbatim
+       QUADRATURE *quadrature;
+
+       if ((quadrature = GetQuadrature(2,3,2)) == NULL)
+              return(1);
+       sum = 0.0;
+       for (l=0; l<Q_NIP(quadrature); l++)
+       {
+        LOCAL_TO_GLOBAL(n,x,Q_LOCAL(quadrature,l),global);
+        (*SourceFunc)(global,val);
+        sum += val * Q_WEIGHT(quadrature,l);
+       }
+       sum = val * AreaOfTriangle;
+       \endverbatim
+
+       \return <ul>
+       <li> pointer to quadrature </li>
+       <li> NULL if the quadrature formula cannot be found </li>
+       </ul>
+     */
+    inline QUADRATURE *GetQuadrature(int dim, int n, int order);
+
+
+    /****************************************************************************/
+    /** \brief Provide a quadrature formula
+
+       \param dim - dimension of the formula
+       \param n - number of corners of the element
+       \param order - order of approximation
+
+       This function returns a pointer to a quadrature formula.
+
+       It is different from GetQuadrature() in the case that the quadrature
+       formula of order 'order' is not exactly available.  'GetQuadrature'
+       just returns the highest order formula, GetQuadratureRule() returns
+       the formula of the smallest degree that integrates exactly until 'order'.
+
+       \section Example
+       \verbatim
+       QUADRATURE *quadrature;
+
+       if ((quadrature = GetQuadratureRule(2,3,2)) == NULL)
+              return(1);
+       sum = 0.0;
+       for (l=0; l<Q_NIP(quadrature); l++)
+       {
+        LOCAL_TO_GLOBAL(n,x,Q_LOCAL(quadrature,l),global);
+        (*SourceFunc)(global,val);
+        sum += val * Q_WEIGHT(quadrature,l);
+       }
+       sum = val * AreaOfTriangle;
+       \endverbatim
+
+       \return <ul>
+       <li> pointer to quadrature </li>
+       <li> NULL if the quadrature formula cannot be found </li>
+       </ul>
+     */
+    /****************************************************************************/
+    inline QUADRATURE *GetQuadratureRule(int dim, int n, int order);
+
+    /****************************************************************************/
+    /** \brief Provide a quadrature formula with axial-symmetric integration points
+
+       \param dim - dimension of the formular
+       \param n - number of corners of the element
+       \param order - order of approximation
+
+       This function returns a pointer to a quadrature formula.  The
+       quadrature rule uses symmetric integration points in 1D and 2D,
+       i.e. the integration points are symmetric w.r.t. the point 0.5 in
+       the 1D case and symmetric to the x- and y-axis in 2D.
+
+       It is different from GetQuadrature in the case that the quadrature
+       formula of order 'order' is not exactly available.  'GetQuadrature'
+       just returns the highest order formula, GetQuadratureRule returns
+       the formula of the smallest degree that integrates exactly until 'order'.
+
+       \section Example
+       \verbatim
+       QUADRATURE *quadrature;
+
+       if ((quadrature = GetSymmetricQuadratureRule(2,3,2)) == NULL)
+              return(1);
+       sum = 0.0;
+       for (l=0; l<Q_NIP(quadrature); l++)
+       {
+        LOCAL_TO_GLOBAL(n,x,Q_LOCAL(quadrature,l),global);
+        (*SourceFunc)(global,val);
+        sum += val * Q_WEIGHT(quadrature,l);
+       }
+       sum = val * AreaOfTriangle;
+       \endverbatim
+
+       \return <ul>
+       <li> pointer to quadrature </li>
+       <li> NULL if the quadrature formula cannot be found </li>
+       </ul>
+     */
+    /****************************************************************************/
+
+    inline QUADRATURE *GetSymmetricQuadratureRule(int dim, int n, int order);
+    inline int         GaussPoints(int dim, int n, int order, DOUBLE_VECTOR *x, GAUSS_POINT *gp);
+
+#include "ugquadratures.cc"
+
+  } // end namespace UGQuadratures
+
+} // end namespace Dune
+
+#endif