6
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ Funny, it seems that NIntegrate thinks Annulus is a Disk in this case. If you create mesh for only half on annulus AnnulusMesh[{0, 0}, {1,2}, {0,Pi}, {36,10}] then the value of NIntegrate is ok. $\endgroup$ – Pinti Dec 3 '19 at 15:54
  • $\begingroup$ @Pinti, this was the right clue. $\endgroup$ – user21 Dec 4 '19 at 8:51
  • 1
    $\begingroup$ @UlrichNeumann This bug is fixed in MeshTools 1.0.1 release. Thank you for using this package! $\endgroup$ – Pinti Dec 6 '19 at 15:16
  • $\begingroup$ @pinti Thanks, I'll update soon. $\endgroup$ – Ulrich Neumann Dec 7 '19 at 9:19
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks.That means , checking the mesh via mesh2D["Wireframe" ] is dangerous? $\endgroup$ – Ulrich Neumann Dec 4 '19 at 8:13
  • $\begingroup$ @UlrichNeumann, no. This is only an issue in NIntegrate because it does remesh the structure you give it. $\endgroup$ – user21 Dec 4 '19 at 8:49
  • $\begingroup$ mesh2D["Wireframe" ] before setting RegionHoles shows a plot with hole. That's what I tried to state in my last comment! $\endgroup$ – Ulrich Neumann Dec 4 '19 at 8:58
  • $\begingroup$ @UlrichNeumann, nothing bad (dangerous) is going to happen when you look at the wire frame. Maybe you mean that the wire frame is misleading? $\endgroup$ – user21 Dec 4 '19 at 9:24
  • 1
    $\begingroup$ @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. $\endgroup$ – Pinti Dec 6 '19 at 15:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for?Browse other questions tagged or ask your own question.