[Dune] [Dune-Commit] dune-grid r4301 - trunk/grid/test
Andreas Dedner
andreas.dedner at mathematik.uni-freiburg.de
Tue Jul 29 15:33:38 CEST 2008
> > * Integral over unit outer normals is not necessarily
> zero for curved surfaces
> > (take, e.g, the upper half sphere). Therefore we
> issue only a warning.
> >
> The integral is zero for closed surfaces, be they curved
> or not. Since
> the boundary of elements is closed I don't see why we
> should only have
> a warning here.
How about the half sphere above the xy plane, e.g. we
only positive z values. All Normals will point downwards
I think, e.g. have negative z-values (in fact I would
assume
that the normal is (0,0,-1) everywhere), i.e., the
integral will not be zero. For dimworld=dimgrid it should
be zero but otherwise...
And even if dimworld=dimgrid an integration error has to
be taken into account which on coarse grids might be quite
a bit more than 1e-8.
>
> > * Added a check whether normals are orthogonal to the
> face. As this check uses
> > the jacobianInverseTransposed-method on the
> "integrationGlocal" geometry, it
> > is disabled by default. To enable it, simply add
> > #define ENABLE_NORMAL_CHECK 1
> > to your check program.
> >
> As far as I know the method jacobianInverseTransposed()
> on
> intersectionGlobal()
> is mandatory. Why make this test optional?
>
You are right, but at the moment most grids do not
implement this method correctly. Until a few months
ago the return value for the jacobianInverseTranspose
was wrong for dimgrid!=dimworld so most grids could
not even implement the correct jit since for the
intersectionGlobal geometry we always have
dimworld!=dimgeometry
Best
Andreas
> greetings,
> Oliver
>
> >
> > Modified: trunk/grid/test/checkintersectionit.cc
> >
>
===================================================================
> > --- trunk/grid/test/checkintersectionit.cc 2008-07-25
> 13:55:19 UTC (rev 4300)
> > +++ trunk/grid/test/checkintersectionit.cc 2008-07-25
> 16:23:08 UTC (rev 4301)
> > @@ -9,9 +9,19 @@
> > \brief Tests for the IntersectionIterator
> > */
> >
> > +struct CheckIntersectionIteratorErrorState
> > +{
> > + unsigned int sumNormalsNonZero;
> > +
> > + CheckIntersectionIteratorErrorState ()
> > + : sumNormalsNonZero( 0 )
> > + {}
> > +};
> > +
> > +
> > +
> > /** \brief Helper routine: Test a general geometry
> > */
> > -
> > template <class GeometryImp>
> > void checkGeometry(const GeometryImp& geometry)
> > {
> > @@ -48,7 +58,11 @@
> >
> > // Check whether point is within the
> intersection
> > if (!geometry.checkInside(testPoint))
> > - DUNE_THROW(GridError, "Test point is not
> within geometry!");
> > + {
> > + std :: cerr << "Test point (" << testPoint
> << ") not within geometry."
> > + << std :: endl;
> > + //DUNE_THROW(GridError, "Test point is not
> within geometry!");
> > + }
> >
> > // Transform to global coordinates
> > FieldVector<ctype, dimworld> global =
> geometry.global(testPoint);
> > @@ -77,9 +91,10 @@
> >
> > /** \brief Test the IntersectionIterator
> > */
> > -template <class GridViewType>
> > +template <class GridViewType, class ErrorState >
> > void checkIntersectionIterator(const GridViewType&
> view,
> > - const typename
> GridViewType::template Codim<0>::Iterator& eIt)
> > + const typename
> GridViewType::template Codim<0>::Iterator& eIt,
> > + ErrorState &errorState
> )
> > {
> > using namespace Dune;
> >
> > @@ -90,9 +105,12 @@
> > const bool checkOutside = (grid.name() !=
> "AlbertaGrid");
> > const typename GridViewType::IndexSet& indexSet =
> view.indexSet();
> >
> > - const int dimworld = GridType::dimensionworld;
> > + typedef typename GridType :: ctype ctype;
> > +
> > + const int dim = GridType :: dimension;
> > + const int dimworld = GridType :: dimensionworld;
> >
> > - FieldVector<double,dimworld> sumNormal(0.0);
> > + FieldVector< ctype,dimworld > sumNormal( ctype( 0 )
> );
> >
> > //
> /////////////////////////////////////////////////////////
> > // Check the types defined by the iterator
> > @@ -128,9 +146,6 @@
> > const QuadratureRule<double, interDim>& quad
> > = QuadratureRules<double,
> interDim>::rule(iIt->intersectionSelfLocal().type(),
> interDim);
> >
> > - for (size_t i=0; i<quad.size(); i++)
> > - sumNormal.axpy(quad[i].weight(),
> iIt->integrationOuterNormal(quad[i].position()));
> > -
> > typedef typename IntersectionIterator::Entity
> EntityType;
> > typedef typename EntityType::EntityPointer
> EntityPointer;
> >
> > @@ -231,14 +246,39 @@
> > // (Ab)use a quadrature rule as a set of test
> points
> > for (size_t i=0; i<quad.size(); i++)
> > {
> > - // check integrationOuterNormal
> > - double det =
>
intersectionGlobal.integrationElement(quad[i].position());
> > - det -=
>
iIt->integrationOuterNormal(quad[i].position()).two_norm();
> > - if( std::abs( det ) > 1e-8 )
> > + const FieldVector< ctype, dim-1 > &pt =
> quad[ i ].position();
> > + FieldVector< ctype, dimworld > normal =
> iIt->integrationOuterNormal( pt );
> > +
> > + sumNormal.axpy( quad[ i ].weight(), normal
> );
> > +
> > + ctype det =
> intersectionGlobal.integrationElement( pt );
> > + if( std :: abs( det - normal.two_norm() ) >
> 1e-8 )
> > {
> > - DUNE_THROW(GridError, "integrationElement
> and length of integrationOuterNormal do no match!");
> > + std :: cerr << "integrationOuterNormal
> yields wrong length: "
> > + << "|n| = " <<
> normal.two_norm() << ", "
> > + << "integrationElement = " <<
> det << std :: endl;
> > }
> > -
> > +
> > +#if ENABLE_NORMAL_CHECK
> > + // Check whether the normal is orthogonal to
> the intersection, i.e.,
> > + // whether (J^-1 * n) = 0. Here J is the
> jacobian of the intersection
> > + // geometry (intersectionGlobal) and n is a
> normal.
> > + //
> > + // This check is disabled by default, since
> not all grid support the
> > + // method jacobianInverseTransposed on a
> face geometry.
> > + const FieldMatrix< ctype, dimworld, dim-1 >
> &jit
> > + =
> intersectionGlobal.jacobianInverseTransposed( pt );
> > + FieldVector< ctype, dim-1 > x( ctype( 0 ) );
> > + jit.umtv( normal, x );
> > + if( x.infinity_norm() > 1e-8 )
> > + {
> > + std :: cerr << "integrationOuterNormal is
> not orthogonal to intersection: "
> > + << "n = " << normal << std ::
> endl;
> > + }
> > +#else
> > +//#warning "Normals are not checked."
> > +#endif
> > +
> > FieldVector<double,dimworld> globalPos =
> intersectionGlobal.global(quad[i].position());
> > FieldVector<double,dimworld> localPos =
>
eIt->geometry().global(intersectionSelfLocal.global(quad[i].position()));
> >
> > @@ -287,12 +327,16 @@
> >
>
DUNE_THROW(GridError,"Entity::hasBoundaryIntersections
> implemented wrong! \n");
> > }
> >
> > - //
>
////////////////////////////////////////////////////////////////////////
> > - // Check whether the integral over the outer
> normal really is zero
> > - //
>
////////////////////////////////////////////////////////////////////////
> > - if( sumNormal.two_norm() > 1e-8 )
> > - DUNE_THROW(GridError,"Sum of
> integrationOuterNormals is " << sumNormal.two_norm() << "
> but it should be Zero!");
> > -
> > + // Chech whether integral over the outer normal is
> zero
> > + //
> > + // Note: This is wrong on curved surfaces (take,
> e.g., the upper half sphere).
> > + // Therefore we only issue a warning.
> > + if( sumNormal.two_norm() > 1e-8 )
> > + {
> > + std :: cout << "Warning: Integral over outer
> normals is nonzero: "
> > + << sumNormal << std :: endl;
> > + ++errorState.sumNormalsNonZero;
> > + }
> > }
> >
> > /** \brief Test both IntersectionIterators
> > @@ -307,9 +351,17 @@
> > ElementIterator eIt = view.template begin<0>();
> > ElementIterator eEndIt = view.template end<0>();
> >
> > + CheckIntersectionIteratorErrorState errorState;
> > for (; eIt!=eEndIt; ++eIt)
> > - checkIntersectionIterator(view, eIt);
> > + checkIntersectionIterator( view, eIt, errorState
> );
> >
> > + if( errorState.sumNormalsNonZero > 0 )
> > + {
> > + std :: cerr << "Warning: Integral over outer
> normals is not always zero."
> > + << std :: endl;
> > + std :: cerr << " This behaviour may be
> correct for entities with"
> > + << " nonzero curature." << std ::
> endl;;
> > + }
> > }
> >
> > template <class GridType>
> >
> > Modified: trunk/grid/test/test-alberta.cc
> >
>
===================================================================
> > --- trunk/grid/test/test-alberta.cc 2008-07-25 13:55:19
> UTC (rev 4300)
> > +++ trunk/grid/test/test-alberta.cc 2008-07-25 16:23:08
> UTC (rev 4301)
> > @@ -1,6 +1,7 @@
> > -
> > #include <config.h>
> >
> > +#define ENABLE_NORMAL_CHECK 1
> > +
> > /*
> >
> > Instantiate Alberta-Grid and feed it to the generic
> gridcheck()
> >
> >
> > _______________________________________________
> > Dune-Commit mailing list
> > Dune-Commit at dune-project.org
> >
>
http://lists.dune-project.org/mailman/listinfo/dune-commit
> >
>
>
> --
>
************************************************************************
> * Oliver Sander ** email:
> sander at mi.fu-berlin.de *
> * Freie Universität Berlin ** phone: + 49 (30) 838
> 75217 *
> * Institut für Mathematik II ** URL :
> page.mi.fu-berlin.de/~sander *
> * Arnimallee 6 **
> -------------------------------------*
> * 14195 Berlin, Germany ** Member of MATHEON
> (www.matheon.de) *
>
************************************************************************
>
>
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune
More information about the Dune
mailing list