[dune-pdelab] [dune-pdelab-commit] dune-pdelab r1558 - trunk/dune/pdelab/localoperator

Christian Engwer christian.engwer at uni-muenster.de
Mon Jul 25 09:41:53 CEST 2011


Hi Felix,

> you almost got me confused, so I wrote an example program illustrating
> the points I made before (works with trunk). It is attached as well as
> its output. It l2 projects a smooth function into a monomial fe space on
> a grid with a single general hexahedra and performs the integrations
> (for various orders of fe space and gaussian rules). Then it compares to
> the corresponding simplex grid. To the interpretation:
>
> Assume you integrate a function with standard finite element basis of
> polynomial order P on a 3d grid with general hexahedra. The local
> shapefunctions are defined on the reference element, where they have
> order P. Hence the function itself has also order P. However, what you
> actually integrate on each cell is the function times the integration
> element (determinant of the jacobian of the transformation). As
> mentioned before, in this case, this function has polynomial order 2.
> Hence you need a gaussian quadrature which works for the order P+2. If
> you have a simplex grid, then the integration element has always
> polynomial order zero. Hence you require different quadrature orders for
> simplices and hexahedrals. Keeping in mind, that the gaussian quadrature
> rule for order 0 and 1, 2 and 3, 4 and 5 are the same, the code results
> confirm the explanation above.

OK, you are right. I forgot about the integration element. I just
thouight about the rhs which would require a highr order integration
due to the transformation onto the reference element.

> Actually in our usual udg approach with monomial Pk basis and marching
> cubes sub triangulation, the transformation from the entity part to the
> bounding box does NOT increase the quadrature order. As this
> transformation is in the worst case a Q1 function it will map e.g. the
> monom x^2 to a function with highest term x^2 * y^2 * z^2 . However, the
> gaussian rule integrating x^k will also integrate x^k * y^k * z^k (this
> is, I think, actually due to the tensor construction of the gaussian
> quadrature rules). However, when using a Qk basis we have to be careful,
> as then the mapping will map e.g. the monom x*y to a function with a
> term like x^2 * y^2.

That's exactly the situation I had in mind. It is the same as the one
with the rhs, although in case of the rhs.

Christian

> Am Freitag, den 22.07.2011, 22:11 +0200 schrieb Christian Engwer:
> > > now I don't exactly know what you mean by "normal". If this operator
> > > would be applied on a grid with general hexahedrals, then the
> > > integration order will have to be higher than on a simplex grid. E.g. in
> > > 3d, the Jacobian Determinant of the transformation has second polynomial
> > > order (first in 2d) and the JacobianITs themselves have first order
> > > (first in 2d too). This increases the order of the actual function
> > > integrated on the ref-element by gaussian quadrature. Actually in case
> > > of UDG, the quadrature order can be slightly reduced, as the jacobianIT
> > > matrix used for the transformation of the fe gradients is constant on
> > > each element due to the affine mapping of the bounding boxes.
> > 
> > no, actually not. You define the shapefnkts and the quadraturepoints
> > on the reference element, thus you don't have to consider the
> > transformation. This doesn't hold anymore for source terms, but this
> > is already a problem if you just have structured grid, but source
> > terms of higher order.
> > 
> > In case of udg you need the additional order, as you only integrate a
> > subset of the reference element of your shape functions.
> > 
> > Christian
> > 
> > > If "normal" for you means simplex of affine, then you are right of
> > > course, but it thats ok for everyone, I'd rather have the more general
> > > formulation.
> > > 
> > > Felix
> > > 
> > > Am Freitag, den 22.07.2011, 17:57 +0200 schrieb Christian Engwer:
> > > > Hi Felix,
> > > > 
> > > > I believe that (for all normal discretizations) the integration order
> > > > does not have to consider the geometry type. This is due to the
> > > > tensor-like construction of the gauss-quadrature. 
> > > > 
> > > > In the case of UDG this is different.
> > > > 
> > > > Christian
> > > > 
> > > > On Fri, Jul 22, 2011 at 05:06:32PM +0200, felix at conan.iwr.uni-heidelberg.de wrote:
> > > > > Author: felix
> > > > > Date: Fri Jul 22 17:06:32 2011
> > > > > New Revision: 1558
> > > > > URL: http://svn.dune-project.org/websvn/listing.php?repname=dune-pdelab&path=/&rev=1558&sc=1
> > > > > 
> > > > > Log:
> > > > > Really choose the right quadrature order automatically (assuming mu
> > > > > and rho to be constant). This should be sufficient for all element
> > > > > types.
> > > > > 
> > > > > Modified:
> > > > >    trunk/dune/pdelab/localoperator/stokesdg.hh
> > > > > 
> > > > > Modified: trunk/dune/pdelab/localoperator/stokesdg.hh
> > > > > ==============================================================================
> > > > > --- trunk/dune/pdelab/localoperator/stokesdg.hh	Fri Jul 22 12:50:50 2011	(r1557)
> > > > > +++ trunk/dune/pdelab/localoperator/stokesdg.hh	Fri Jul 22 17:06:32 2011	(r1558)
> > > > > @@ -110,9 +110,12 @@
> > > > >                  typedef typename LFSV::Traits::SizeType size_type;
> > > > >  
> > > > >                  // select quadrature rule
> > > > > -                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = v_order + superintegration_order;
> > > > >                  Dune::GeometryType gt = eg.geometry().type();
> > > > > +                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > +                const int det_jac_order = gt.isSimplex() ?  0 : (dim-1);
> > > > > +                // quad order is velocity order + det_jac order + superintegration
> > > > > +                const int qorder = v_order + det_jac_order + superintegration_order;
> > > > > +
> > > > >                  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
> > > > >                  
> > > > >                  // loop over quadrature points
> > > > > @@ -190,9 +193,11 @@
> > > > >                  typedef typename LFSV::Traits::SizeType size_type;
> > > > >  
> > > > >                  // select quadrature rule
> > > > > +                Dune::GeometryType gtface = ig.geometry().type();
> > > > >                  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = 2*v_order - 1 + superintegration_order;
> > > > > -                Dune::GeometryType gtface = ig.geometryInInside().type();
> > > > > +                const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
> > > > > +                const int jac_order = gtface.isSimplex() ? 0 : 1;
> > > > > +                const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
> > > > >  
> > > > >                  const RF penalty_factor = prm.getFaceIP(ig);
> > > > > @@ -241,7 +246,7 @@
> > > > >                              }
> > > > >                          }
> > > > >                          //================================================//
> > > > > -                        // \mu \int \sigma / |\gamma|^\beta v u_0
> > > > > +                        // \int \sigma / |\gamma|^\beta v u_0
> > > > >                          //================================================//
> > > > >                          factor = penalty_factor * weight;
> > > > >                          for (unsigned int i=0;i<vsize;++i) 
> > > > > @@ -321,10 +326,11 @@
> > > > >                  typedef typename LFSV::Traits::SizeType size_type;
> > > > >  
> > > > >                  // select quadrature rule
> > > > > -                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = 2*v_order - 2 + superintegration_order;
> > > > > -
> > > > >                  Dune::GeometryType gt = eg.geometry().type();
> > > > > +                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > +                const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
> > > > > +                const int jac_order = gt.isSimplex() ? 0 : 1;
> > > > > +                const int qorder = 2*v_order - 2 + 2*jac_order + det_jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
> > > > >                  
> > > > >                  // loop over quadrature points
> > > > > @@ -342,7 +348,8 @@
> > > > >                      BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
> > > > >                                              eg.geometry(), local, grad_phi_v);
> > > > >  
> > > > > -                    const RF weight = it->weight() * eg.geometry().integrationElement(it->position());
> > > > > +                    const RF detj = eg.geometry().integrationElement(it->position());
> > > > > +                    const RF weight = it->weight() * detj;
> > > > >                      
> > > > >                      //================================================//
> > > > >                      // \int (mu*grad_u*grad_v)
> > > > > @@ -431,9 +438,10 @@
> > > > >                  typedef typename LFSV::Traits::SizeType size_type;
> > > > >  
> > > > >                  // select quadrature rule
> > > > > +                Dune::GeometryType gtface = ig.geometry().type();
> > > > >                  const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
> > > > > -                const int qorder = 2*v_order - 1 + superintegration_order;
> > > > > -                Dune::GeometryType gtface = ig.geometryInInside().type();
> > > > > +                const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
> > > > > +                const int qorder = 2*v_order + det_jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
> > > > >  
> > > > >                  const RF penalty_factor = prm.getFaceIP(ig);
> > > > > @@ -672,8 +680,9 @@
> > > > >  
> > > > >                  // select quadrature rule
> > > > >                  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = 2*v_order - 1 + superintegration_order;
> > > > > -                Dune::GeometryType gtface = ig.geometryInInside().type();
> > > > > +                Dune::GeometryType gtface = ig.geometry().type();
> > > > > +                const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
> > > > > +                const int qorder = 2*v_order + det_jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
> > > > >  
> > > > >                  // evaluate boundary condition type
> > > > > @@ -832,9 +841,11 @@
> > > > >                  typedef typename LFSV::Traits::SizeType size_type;
> > > > >  
> > > > >                  // select quadrature rule
> > > > > -                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = 3*v_order - 1 + superintegration_order;
> > > > >                  Dune::GeometryType gt = eg.geometry().type();
> > > > > +                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > +                const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
> > > > > +                const int jac_order = gt.isSimplex() ? 0 : 1;
> > > > > +                const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
> > > > >                  
> > > > >                  // loop over quadrature points
> > > > > @@ -935,8 +946,10 @@
> > > > >  
> > > > >                  // select quadrature rule
> > > > >                  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = 3*v_order - 1 + superintegration_order;
> > > > >                  Dune::GeometryType gt = eg.geometry().type();
> > > > > +                const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
> > > > > +                const int jac_order = gt.isSimplex() ? 0 : 1;
> > > > > +                const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
> > > > >                  
> > > > >                  // loop over quadrature points
> > > > > @@ -1056,10 +1069,10 @@
> > > > >                  typedef typename LFSV::Traits::SizeType size_type;
> > > > >  
> > > > >                  // select quadrature rule
> > > > > -                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > -                const int qorder = 2*v_order + superintegration_order;
> > > > > -
> > > > >                  Dune::GeometryType gt = eg.geometry().type();
> > > > > +                const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
> > > > > +                const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
> > > > > +                const int qorder = 2*v_order + det_jac_order + superintegration_order;
> > > > >                  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
> > > > >                  
> > > > >                  // loop over quadrature points
> > > > > 
> > > > > _______________________________________________
> > > > > dune-pdelab-commit mailing list
> > > > > dune-pdelab-commit at dune-project.org
> > > > > http://lists.dune-project.org/mailman/listinfo/dune-pdelab-commit
> > > > > 
> > > > 
> > > > _______________________________________________
> > > > dune-pdelab mailing list
> > > > dune-pdelab at dune-project.org
> > > > http://lists.dune-project.org/mailman/listinfo/dune-pdelab
> > > 
> > > -- 
> > > Felix Heimann
> > > Universität Heidelberg
> > > Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
> > > Arbeitsgruppe Paralleles Rechnen
> > > IWR 368, Raum 422
> > > Tel: 06221 / 54 8881
> > > 
> > > 
> 
> -- 
> Felix Heimann
> Universität Heidelberg
> Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
> Arbeitsgruppe Paralleles Rechnen
> IWR 368, Raum 422
> Tel: 06221 / 54 8881
> 

> // -*- tab-width: 4; indent-tabs-mode: nil -*-
> /** \file 
>     \brief Solve Poisson problem on various grids (sequential)
> */
> #ifdef HAVE_CONFIG_H
> #include "config.h"     
> #endif
> #include<iostream>
> #include<vector>
> #include<map>
> #include<string>
> #include<dune/common/mpihelper.hh>
> #include<dune/common/exceptions.hh>
> #include<dune/common/fvector.hh>
> #include<dune/common/float_cmp.hh>
> #include<dune/common/static_assert.hh>
> #include<dune/grid/yaspgrid.hh>
> #include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
> #include<dune/istl/bvector.hh>
> #include<dune/istl/operators.hh>
> #include<dune/istl/io.hh>
> 
> #if HAVE_UG 
> #include<dune/grid/uggrid.hh>
> #endif
> 
> #include<dune/pdelab/backend/istlvectorbackend.hh>
> #include<dune/pdelab/finiteelementmap/monomfem.hh>
> #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
> #include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
> #include<dune/pdelab/gridfunctionspace/interpolate.hh>
> #include<dune/pdelab/constraints/constraints.hh>
> #include<dune/pdelab/common/function.hh>
> #include<dune/pdelab/common/functionutilities.hh>
> #include<dune/pdelab/common/vtkexport.hh>
> #include <dune/grid/common/gridfactory.hh>
> 
> // declaration of a basic unit cube that uses the GridFactory
> template< int dim >
> class BasicGeneralCube;
> 
> // unit cube in 3 dimensions with two variants: tetraheda and hexahedra
> template<>
> class BasicGeneralCube< 3 >
> {
> protected:
>   template< class Grid >
>   static void insertVertices ( Dune::GridFactory< Grid > &factory )
>   {
>     Dune::FieldVector< double, 3 > pos;
> 
>     pos[0] = 0;  pos[1] = 0;  pos[2] = 0;    factory.insertVertex(pos);
>     pos[0] = 1;  pos[1] = 0;  pos[2] = 0;    factory.insertVertex(pos);
>     pos[0] = 0;  pos[1] = 1.5;  pos[2] = 0;    factory.insertVertex(pos);
>     pos[0] = 2;  pos[1] = 1.8;  pos[2] = 0;    factory.insertVertex(pos);
>     pos[0] = 0;  pos[1] = 0;  pos[2] = 1;    factory.insertVertex(pos);
>     pos[0] = 3;  pos[1] = 0;  pos[2] = 2;    factory.insertVertex(pos);
>     pos[0] = 0;  pos[1] = 3;  pos[2] = 2.5;    factory.insertVertex(pos);
>     pos[0] = 4;  pos[1] = 1.3;  pos[2] = 3;    factory.insertVertex(pos);
>   }
> 
>   template< class Grid >
>   static void insertSimplices ( Dune::GridFactory< Grid > &factory )
>   {
>     const Dune::GeometryType type( Dune::GeometryType::simplex, 3 );
>     std::vector< unsigned int > cornerIDs( 4 );
> 
>     cornerIDs[0] = 0;  cornerIDs[1] = 1;  cornerIDs[2] = 2;  cornerIDs[3] = 4;
>     factory.insertElement( type, cornerIDs );
> 
>     cornerIDs[0] = 1;  cornerIDs[1] = 3;  cornerIDs[2] = 2;  cornerIDs[3] = 7;
>     factory.insertElement( type, cornerIDs );
> 
>     cornerIDs[0] = 1;  cornerIDs[1] = 7;  cornerIDs[2] = 2;  cornerIDs[3] = 4;
>     factory.insertElement( type, cornerIDs );
> 
>     cornerIDs[0] = 1;  cornerIDs[1] = 7;  cornerIDs[2] = 4;  cornerIDs[3] = 5;
>     factory.insertElement( type, cornerIDs );
> 
>     cornerIDs[0] = 4;  cornerIDs[1] = 7;  cornerIDs[2] = 2;  cornerIDs[3] = 6;
>     factory.insertElement( type, cornerIDs );
>   }
> 
>   template< class Grid >
>   static void insertCubes ( Dune::GridFactory< Grid > &factory )
>   {
>     const Dune::GeometryType type( Dune::GeometryType::cube, 3 );
>     std::vector< unsigned int > cornerIDs( 8 );
>     for( int i = 0; i < 8; ++i )
>       cornerIDs[ i ] = i;
>     factory.insertElement( type, cornerIDs );
>   }
> };
> 
> template< int dim, int variant >
> class UGGeneralCube : public BasicGeneralCube< dim >
> {
> public:
>   typedef Dune::UGGrid<dim> GridType;
> 
> private:  
>   GridType* grid_;
> 
> public:
>   UGGeneralCube ()
>   {
>     Dune::GridFactory< GridType > factory;
>     BasicGeneralCube< dim >::insertVertices( factory );
>     if( variant == 1 )
>       BasicGeneralCube< dim >::insertCubes( factory );
>     else if( variant == 2 )
>       BasicGeneralCube< dim >::insertSimplices( factory );
>     else
>       DUNE_THROW( Dune::NotImplemented, "Variant " 
>                   << variant << " of UG unit cube not implemented." );
>     grid_ = factory.createGrid();
>   }
> 
>   ~UGGeneralCube ()
>   {
>     delete grid_;
>   }
> 
>   GridType &grid ()
>   {
>     return *grid_;
>   }
> };
> //===============================================================
> // Problem setup and solution 
> //===============================================================
> 
> // function for Dirichlet boundary conditions and initialization
> template<typename GV, typename RF>
> class G
>   : public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
>                                                   G<GV,RF> >
> {
> public:
>   typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
>   typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,G<GV,RF> > BaseT;
> 
>   G (const GV& gv) : BaseT(gv) {}
>   inline void evaluateGlobal (const typename Traits::DomainType& x, 
> 							  typename Traits::RangeType& y) const
>   {
>     typename Traits::DomainType center;
>     for (int i=0; i<GV::dimension; i++) center[i] = 0.5;
>     center -= x;
> 	y = exp(-center.two_norm2());
>   }
> };
> 
> 
> // generate a P1 function and output it
> template<typename GV, int k> 
> void integrate (const GV& gv, const int quadorder, std::string filename, Dune::GeometryType gtype)
> {
>   // constants and types
>   typedef typename GV::Grid::ctype DF;
>   const int dim = GV::dimension;
>   typedef double R;
> 
>   // make finite element map
>   typedef Dune::PDELab::MonomLocalFiniteElementMap<DF,R,dim,k> FEM;
>   FEM fem(gtype);
> 
>   // make grid function space
>   typedef Dune::PDELab::ISTLVectorBackend<1> VBE;
>   typedef Dune::PDELab::NoConstraints CON;
>   typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE> GFS; 
>   GFS gfs(gv,fem);
> 
>   // make coefficent Vector and initialize it from a function
>   typedef typename Dune::PDELab::BackendVectorSelector<GFS,R>::Type V;
>   V x0(gfs);
>   x0 = 0.0;
>   // for(unsigned int i=0; i<gfs.size(); ++i)
>   //   VBE::access(x0,i) = (double) i;
> 
>   typedef G<GV,R> GType;
>   GType g(gv);
>   Dune::PDELab::interpolate(g,gfs,x0);
> 
>   // make discrete function object
>   typedef Dune::PDELab::DiscreteGridFunction<GFS,V> DGF;
>   DGF dgf(gfs,x0);
> 
>   typename DGF::Traits::RangeType integral(0);
>   Dune::PDELab::integrateGridFunction(dgf,integral,quadorder);
>   std::cout << "Integration with order " << quadorder << " : " << integral << std::endl;
> 
>   // output grid function with VTKWriter
>   Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTKOptions::conforming);
>   vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(dgf,"solution"));
>   vtkwriter.write(filename,Dune::VTKOptions::ascii);
> }
> 
> 
> //===============================================================
> // Main program with grid setup
> //===============================================================
> 
> template<typename GV>
> void test_integrations(const GV& gv, Dune::GeometryType gtype)
> {
>   // solve problem
>   const int maxq = 6;
>   std::cout << "FE order = 1" << std::endl;
>   for(int q=0; q<maxq; ++q)
>     integrate<GV,1>(gv,q,"integration_domain",gtype);
> 
>   std::cout << "FE order = 2" << std::endl;
>   for(int q=0; q<maxq; ++q)
>     integrate<GV,2>(gv,q,"integration_domain",gtype);
> 
>   std::cout << "FE order = 3" << std::endl;
>   for(int q=0; q<maxq; ++q)
>     integrate<GV,3>(gv,q,"integration_domain",gtype);
> 
> }
> 
> int main(int argc, char** argv)
> {
>   try{
>     //Maybe initialize Mpi
>     Dune::MPIHelper::instance(argc, argv);
> 
>     // UG Pk 2D test
> #if HAVE_UG
> 
>     std::cout << "Integration with cubes:" << std::endl;
>     {
>       // make grid 
>       const int dim = 3;
>       const int variant = 1;
>       typedef UGGeneralCube<dim,variant> Cube;
>       typedef Cube::GridType Grid;
>       Cube cube;
>       Grid & grid = cube.grid();
> 
>       // get view
>       typedef Grid::LeafGridView GV;
>       const GV& gv=grid.leafView(); 
>  
>       Dune::GeometryType bbox_type; bbox_type.makeCube(dim);
>       test_integrations(gv,bbox_type);
>     }
> 
>     std::cout << "Integration with simplices:" << std::endl;
>     {
>       // make grid 
>       const int dim = 3;
>       const int variant = 2;
>       typedef UGGeneralCube<dim,variant> Cube;
>       typedef Cube::GridType Grid;
>       Cube cube;
>       Grid & grid = cube.grid();
> 
>       // get view
>       typedef Grid::LeafGridView GV;
>       const GV& gv=grid.leafView(); 
>  
>       Dune::GeometryType bbox_type; bbox_type.makeSimplex(dim);
>       test_integrations(gv,bbox_type);
>     }
> 
> #endif
> 
> 	// test passed
> 	return 0;
>   }
>   catch (Dune::Exception &e){
>     std::cerr << "Dune reported error: " << e << std::endl;
> 	return 1;
>   }
>   catch (...){
>     std::cerr << "Unknown exception thrown!" << std::endl;
> 	return 1;
>   }
> }

> Integration with cubes:
> DimX=1, DimY=1, DimZ=1
> DimX=1, DimY=1, DimZ=1
> FE order = 1
> Integration with order 0 : 3.04346
> Integration with order 1 : 3.04346
> Integration with order 2 : 2.12961
> Integration with order 3 : 2.12961
> Integration with order 4 : 2.12961
> Integration with order 5 : 2.12961
> FE order = 2
> Integration with order 0 : 3.22151
> Integration with order 1 : 3.22151
> Integration with order 2 : 2.02376
> Integration with order 3 : 2.02376
> Integration with order 4 : 2.0234
> Integration with order 5 : 2.0234
> FE order = 3
> Integration with order 0 : 3.47061
> Integration with order 1 : 3.47061
> Integration with order 2 : 1.90052
> Integration with order 3 : 1.90052
> Integration with order 4 : 2.04595
> Integration with order 5 : 2.04595
> 
> 
> Integration with simplices:
> DimX=1, DimY=1, DimZ=1
> DimX=1, DimY=1, DimZ=1
> FE order = 1
> Integration with order 0 : 2.27982
> Integration with order 1 : 2.27982
> Integration with order 2 : 2.27982
> Integration with order 3 : 2.27982
> Integration with order 4 : 2.27982
> Integration with order 5 : 2.27982
> FE order = 2
> Integration with order 0 : 1.43913
> Integration with order 1 : 1.43913
> Integration with order 2 : 1.94343
> Integration with order 3 : 1.94343
> Integration with order 4 : 1.94343
> Integration with order 5 : 1.94343
> FE order = 3
> Integration with order 0 : 1.66484
> Integration with order 1 : 1.66484
> Integration with order 2 : 1.94556
> Integration with order 3 : 1.96563
> Integration with order 4 : 1.96563
> Integration with order 5 : 1.96563





More information about the dune-pdelab mailing list