6

Fixed in Version 12.1.1

In a very simple example I try to describe an annulus using MeshTools package

mesh2D = AnnulusMesh[{0, 0}, {    1 ,  2 }, {-Pi, Pi}, {36, 10}];

The area given by

{mesh2D["MeshElementMeasure"][[1]] // Total, Pi (2^2 - 1^2)} // N
(*{9.377, 9.42478}*)

differs significantly from evaluation of NIntegrate inside the element mesh

NIntegrate[1, Element[{x, y}, mesh2D]]
(*12.5027*)

What's wrong here? Thanks!

user21
  • 39,710
  • 8
  • 110
  • 167
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55

1 Answers1

8

Fixed in version 12.1.1

Needs["MeshTools`"]
mesh2D = AnnulusMesh[{0, 0}, {1, 2}, {-Pi, Pi}, {36, 10}];
NIntegrate[1, Element[{x, y}, mesh2D]]
9.37700159401424`

Old Answer:

AnnulusMesh does not set the region hole of the mesh region. Then, in NIntegrate the mesh when the mesh is re-meshed that region hole is fully meshed.

Needs["MeshTools`"]
mesh2D = AnnulusMesh[{0, 0}, {1, 2}, {-Pi, Pi}, {36, 10}];
mesh2D["RegionHoles"]
Automatic

SetRegionHoles[mesh2D, {{0, 0}}]
{{0., 0.}}


NIntegrate[1, Element[{x, y}, mesh2D]]
9.37700159401424`

Note that an annulus meshed with ToElementMesh automatically sets that region hole property.

ToElementMesh[Annulus[{0, 0}, {1, 2}]]["RegionHoles"]
{{2.5326962749261384`*^-16, 2.7929047963226594`*^-16}}

I think the best way forward is to add a region hole to AnnulusMesh. I'll have a look how time consuming it would be for NIntegrate to auto search for region holes if the mesh["RegionHoles"] is Automatic; but that may be prohibitive.

In other words this is happening:

ToElementMesh[mesh2D]["Wireframe"]

enter image description here

But you want this to happen:

ToElementMesh[mesh2D, "RegionHoles" -> {{0, 0}}]["Wireframe"]

enter image description here

The reason NIntegrate converts quad and hex meshes to triangle and tet meshes is that the main mechanism for NIntegrate is to do adaptive refinement which is only available for triangle and tet meshes. So for a quad or hex meshes an additional cost of conversion comes into play because we want this to work and give good results:

Needs["MeshTools`"]
mesh2D = AnnulusMesh[{0, 0}, {1, 2}, {-Pi, Pi}, {36, 10}];
SetRegionHoles[mesh2D, {{0, 0}}];
nr = ToNumericalRegion[Annulus[{0, 0}, {1, 2}]];
SetNumericalRegionElementMesh[nr, mesh2D];
\[Pi] (2^2 - 1^1) - FEMNIntegrate[1, {x, y}, nr]
-6.6228668185175366`*^-6

Note, how the quality is much better than for the original annulus mesh. Probably the design of AnnulusMesh could be improved by allowing

AnnulusMesh[Annulus[{0,0},{1,2}],{-Pi,Pi},{36,10}]

Because then that same symbolic description used for the creation of an AnnulusMesh could be used to create the numerical region.

user21
  • 39,710
  • 8
  • 110
  • 167
  • Thanks.That means , checking the mesh via mesh2D["Wireframe" ] is dangerous? – Ulrich Neumann Dec 04 '19 at 08:13
  • @UlrichNeumann, no. This is only an issue in NIntegrate because it does remesh the structure you give it. – user21 Dec 04 '19 at 08:49
  • mesh2D["Wireframe" ] before setting RegionHoles shows a plot with hole. That's what I tried to state in my last comment! – Ulrich Neumann Dec 04 '19 at 08:58
  • @UlrichNeumann, nothing bad (dangerous) is going to happen when you look at the wire frame. Maybe you mean that the wire frame is misleading? – user21 Dec 04 '19 at 09:24
  • Yes, misleading, it shows an annulus! – Ulrich Neumann Dec 04 '19 at 09:52
  • @UlrichNeumann, the visualization just plots the elements present but the remeshing needs information if the inner boundary is an internal boundary or if it is a region hole. And that went missing. – user21 Dec 04 '19 at 10:00
  • @user21 Thanks for this clarification. I didn't know about this effect of "RegionHoles" on NIntegrate. I will fix this in MeshTools package immediately. But why does NIntegrate need to do internal remeshing if I pass it a nice ElementMesh anyway? – Pinti Dec 04 '19 at 13:13
  • @Pinti, see update. – user21 Dec 05 '19 at 07:30
  • @Pinit, I just put in a fix. If a mesh has "RegionHoles" Automatic then the region hole computation is triggered. This would have avoided the issue seen here; but still I believe it is better to add "RegionHoles" explicitly to your mesh. – user21 Dec 05 '19 at 14:55
  • 1
    @user21 Thank you for all explanations. I have already updated MeshTools package with bugfix. In a near future I am going to add "MeshOrder"->2 option for AnnulusMesh so its area can be computed with higher accuracy. – Pinti Dec 06 '19 at 15:19