[dune-pdelab] [dune-pdelab-commit] dune-pdelab r1558 - trunk/dune/pdelab/localoperator
Felix Heimann
felix.heimann at iwr.uni-heidelberg.de
Fri Jul 22 18:15:29 CEST 2011
Hi Christian,
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.
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