[dune-pdelab] [dune-pdelab-commit] dune-pdelab r1558 - trunk/dune/pdelab/localoperator
Christian Engwer
christian.engwer at uni-muenster.de
Fri Jul 22 22:11:12 CEST 2011
> 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
>
>
More information about the dune-pdelab
mailing list