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

Christian Engwer christian.engwer at uni-muenster.de
Fri Jul 22 17:57:52 CEST 2011


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
> 




More information about the dune-pdelab mailing list