[Dune] gridcheck.cc : Integral over outer normals is not always zero.

[Aleksejs Fomins]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Dear Dune,

I have performed preliminary tests for this problem. The curvilinear geometry go up to 5th order interpolation. Experiments show that integrating the normal can require up to order 9 quadrature, for the cases explored, but perhaps even higher for cases not considered.
I will file a bug report tomorrow, requesting to edit the gridcheck.

Speaking of which, another question. Has any work been done on adaptive quadrature? Since curvilinear geometry frequently encounters the need to integrate the square root of jacobian, which is not a polynomial, it is hard to guess how high of an interpolation order will be necessary to integrate it to the desired accuracy. So far we attempt to consequently run all quadrature orders until a good enough order has been reached. But it would be better if the quadrature points could be reused in a hierarchic manner.

Greetings,
Aleksejs





On 19/03/15 09:01, Christian Engwer wrote:
>     const Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre;
>     const Dune::QuadratureRule< ctype, Intersection::mydimension > &quadrature
>       = Dune::QuadratureRules< ctype, Intersection::mydimension >::rule( intersection.type(), 3, qt );
>     for( std::size_t i = 0; i < quadrature.size(); ++i )
>       sumNormal.axpy( quadrature[ i ].weight(), intersection.integrationOuterNormal( quadrature[ i ].position() ) );
> 
> The test currently requests a 3rd order quadrature rule. I guess in your case this is not sufficcient, so you could simply increase the order and check again.
> 
> Christian
> 
> On Thu, Mar 19, 2015 at 08:46:24AM +0100, Aleksejs Fomins wrote:
>> Dear Martin,
>>
>> Thank you for your reply. I am using a volume grid. I misunderstood the warning. The confusion came from the fact that I use curvilinear elements, not a curvilinear 2D world, so I was puzzled as of how the curvature of the elements can affect the divergence theorem, and of course it can't.
>>
>> So far my normals seemed fine, I was able to calculate the integral in Gauss theorem.
>>
>> I will replicate this exact normal integral test, and if it passes, I will report a bug.
>>
>> Thanks,
>> Aleksejs
>>
>>
>> On 18/03/15 20:26, Martin Nolte wrote:
>>> Hi Aleksejs,
>>>
>>> the divergence theorem only holds for volumes, not for surfaces. On a surface, the surface integral over all normals equals the integral over the mean curvature. And this is exactly what the message states: In case of a surface grid, the integral might be non-zero if the mean curvature of is nonzero.
>>>
>>> If you encounter this message for a volume grid (i.e., a grid with dimgrid = dimworld), then something went wrong. In this case, the divergence theorem holds and, as you stated, the surface integral over all normals should be zero. Apart from a bug in your normal implementation, this might also result from an insufficient quadrature order. Would you consider filing a bug report in the latter case?
>>>
>>> Best,
>>>
>>> Martin
>>>
>>> PS: I am merely the author of the warning, not of the test. Originally, the test simply failed for surface grids with element geometries of nonzero mean curvature (which is plain wrong).
>>>
>>> On 03/18/2015 04:23 PM, Aleksejs Fomins wrote:
>>> Dear Dune,
>>>
>>> When running the gridcheck.cc I notice I encounter the following warning:
>>>
>>> -- Checking Intersection Iterator
>>> Warning: Integral over outer normals is not always zero.
>>>           This behaviour may be correct for entities with nonzero curvature.
>>> Warning: Integral over outer normals is not always zero.
>>>           This behaviour may be correct for entities with nonzero curvature.
>>>
>>> Could the person who wrote this test please explain what this means.
>>>
>>> I assume that this refers to the integral Int(vec{n} dS) over the surface of an element.
>>> If this is indeed the case, then this integral should be zero by divergence theorem.
>>> In particular, could you explain the case in which the integral would not be zero
>>>
>>> Greetings,
>>> Aleksejs
>>>>
>>>> _______________________________________________
>>>> 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
>>
> 
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.12 (GNU/Linux)

iQIcBAEBAgAGBQJVCwWXAAoJEDkNM7UOEwMZRpsP+wVofhrIn5OEDnJN8NMHmh+P
cIw3qCxrNTSyLhXLR1lee/6CV+leEieI9rdTY/aKD2aSpIFgo/2X1UGvYqj9UW7P
z/p2zEU1axlqFVputTOYyIDynV3Nuw5YtQuIKMqFnNvnXaVjftqTGQ0YbNWE2AJ7
g7l6tovNdLGaxpFrYoajiu3taYo7PWIRbRD0a0sYpfXmZhQHriL6GoWd/Q2ovtPL
XyOJpy2ND3U/Ax2Pgjsp7F3sYVeizWmm0LevByAy+WNywXRi1pIqRs+WH/zgjrYj
yLx8cgOn7eN5/qQTDYVy0SlxcmHXALUWqPDm2xwGbHsXM1cxXI7vWk0yFBpfNWAF
oZSgx5L0TXeZkWGbnkIOiKbN7xpVasng9kXqxahfdB2r4YDNX2ZanUoH2GfNRUPb
uw4Bw+IF2DFD4pk+PbOzn9MAP5AUcAWwaZ5NdHzqyfE+G/cPRKggPJUolxzFPNNg
pMCkBimJL8CXGc5zd0289f/ECC8XslV2FdW7Q/39WBBU74psTmV+GlmuAIvH3MvY
YtHPGnMEVQ6wEGgGl8GIkRnTDv0LTIcL2umuylnHJzXb8pt2PJ06KsLx2F90wWjQ
A6zx8qbdH8wTI3hAFvvJXvaHv4xTJx+SIRQ3slDpcR8qkK1H+zzK0r2ArhnswAQz
SB2xP8vPjLiqG3EP7ce8
=lUDM
-----END PGP SIGNATURE-----