[dune-pdelab] [dune-pdelab-commit] dune-pdelab r1558 - trunk/dune/pdelab/localoperator
Felix Heimann
felixheimann at gmx.de
Sun Jul 24 21:14:01 CEST 2011
Hi Christian,
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.
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.
Best Regards,
Felix
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: integration.cc
Type: text/x-c++src
Size: 8427 bytes
Desc: not available
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20110724/53964ea7/attachment.cc>
-------------- next part --------------
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