[Dune] dune-grid r4301 - trunk/grid/test
Oliver Sander
sander at mi.fu-berlin.de
Tue Jul 29 17:32:38 CEST 2008
Hi Andreas!
> 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.
>
>
I have settled this in a private thread with Martin. You
guys are right, the integral is only zero in flat spaces.
>
> 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
>
But tests exist to tell me what works and what doesn't.
Nobody looks at warnings, and somebody not looking closely
will think that all tests pass and hence everything works.
Only to find out later that it is not so. I prefer to be honest
about what doesn't work and rather try to implement the
missing features.
:-)
Oliver
P.S.: Yes I know, some missing features are supposed to be
implemented by me...
> 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
>>
>
>
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune
>
--
************************************************************************
* 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