[Dune] [Dune-Commit] dune-grid r4301 - trunk/grid/test
Oliver Sander
sander at mi.fu-berlin.de
Mon Jul 28 09:38:59 CEST 2008
Hi!
I have a few comments regarding this commit:
mnolte at dune-project.org schrieb:
> Author: mnolte
> Date: 2008-07-25 18:23:08 +0200 (Fri, 25 Jul 2008)
> New Revision: 4301
>
> Modified:
> trunk/grid/test/checkintersectionit.cc
> trunk/grid/test/test-alberta.cc
> Log:
> * 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.
> * 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?
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) *
************************************************************************
More information about the Dune
mailing list