6

Bug introduced in 10.0 and fixed in 10.2

Could any one please confirm the following bug in mathematica 10.0.2 ?

If I define this ball

Ω = ImplicitRegion[0 <= x^2 + y^2 + z^2 <= 1, {x, y, z}];

and try to compute its edge:

  Needs["NDSolve`FEM`"];
  mesh = ToElementMesh[DiscretizeRegion[Ω]]

This crashes the kernel?

whereas this does not:

    mesh = ToElementMesh[DiscretizeRegion[Ball[]]]
    mesh["Wireframe"]

Mathematica graphics

This is related to this question

user21
  • 39,710
  • 8
  • 110
  • 167
chris
  • 22,860
  • 5
  • 60
  • 149
  • 1
    I'm on v. 10.1, so I can't confirm your behavior on 10.0.2, but I do obtain an error (TetGenTetrahedralize::reterr: Tetrahedralize returned an error, 3. >>), although it does not crash my kernel. – MarcoB Jun 09 '15 at 18:08
  • @MarcoB thanks! – chris Jun 09 '15 at 18:09
  • No problem! I also noticed that, even though ToElementMesh completes evaluation with the Ball[] example, the ElementMesh object that it returns is still wonky: if I try and extract any of its properties, it returns: ElementMesh::noinfo: Input expression ElementMesh[{{-0.999955,0.999955},{-1.,1.},{-1.,1.}},{TetrahedronElement[<8343>]}] contains insufficient information to interpret the result. >>. Is that the same for you? – MarcoB Jun 09 '15 at 18:11
  • I'm on 10.1 too and though I get an error, it still returns the ElementMesh which can be displayed using mesh["Wireframe"] – RunnyKine Jun 09 '15 at 18:13
  • @MarcoB no; no 10.0.2 I get mesh["Properties"] {BoundaryConnectivity,BoundaryElements,Bounds,ContinuousElementConnectivity,Coordinates,ElementConnectivity,EmbeddingDimension,HangingNodes,MeshElements,MeshOrder,PointElements,Properties,Quality,RegionHoles,VertexBoundaryConnectivity,VertexElementConnectivity} – chris Jun 09 '15 at 18:13
  • 2
    You are actually correct. I was trying to extract properties from the ElementMesh object returned by ToElementMesh. Apparently one can't do that without assigning it to a variable first. I didn't know that. Now that I do, I've gone back to your own original Omega region: even though I still obtain that error I mentioned, I can extract all properties from it just fine. I guess this means that whatever was wrong in v. 10.0.2 has been at least partially fixed in v. 10.1 . Sorry for the confusion! – MarcoB Jun 09 '15 at 18:21
  • @MarcoB ToElementMesh[DiscretizeRegion[Ball[]]]["Wireframe"] works for me, for any other property, without assigning the mesh to a variable. -- The difference in Through[{Mean, StandardDeviation}[#]] & /@ mesh["Quality"] is quite a bit worse for the ImplicitRegion, but that's probably to be expected. – Michael E2 Jun 09 '15 at 20:00
  • @MarcoB could you please try {ev, if, mesh} = helmholzSolve3D[\[CapitalOmega], 4, MaxCellMeasure -> 0.25] with this code http://mathematica.stackexchange.com/a/85552/1089 on 10.1 ? – chris Jun 09 '15 at 20:34
  • I did and it failed due to running out of memory. I posted the results as an improper answer below. I'll remove when no more useful. – MarcoB Jun 09 '15 at 20:50

2 Answers2

8

There are a few things to be said here. First, it should not crash and a future version will behave better in this scenario. What I am not quite sure abuot is why the call ToElementMesh[DiscretizeRegion[...]] in the first place. Note that both

mesh = ToElementMesh[\[CapitalOmega]];

and 

DiscretizeRegion[\[CapitalOmega]];

work fine. In your example you first make a mesh region and then remesh to make an ElementMesh. (Since DiscretizeRegion will return a first odrer mesh ToEM[DR[...]] will (should) return a first order mesh. You can get that diretly with ToElementMesh[region, "MeshOrder"->1]. Additionally, the quality of the mesh will not be as good as it could be since ToEM has no way of knowing what the exact boundary is in ToEM[DR[...]])

Now for the real problem: When you do a direct conversion from the mesh region to an element mesh via you get a message:

mr = DiscretizeRegion[\[CapitalOmega]];
mr["MakeRepresentation"["ElementMesh"]]

ElementMesh::femimq: "The element mesh has insufficient quality of -0.0124712. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements."

So the mesh region has bad elements. This is fixed in the development version. It might be good to have this message also in the ToElementMesh[DiscretizeRegion[...]] case such that one knows that the returned element mesh is wonky.

Another note, ploting a "Wireframe" should work in many cases - even if the mesh quality is less then zero. The reason for that is that then one can still look at the mesh and possibly find the bad elements. While this may not be possible in this case, it's useful for smaller meshes.

user21
  • 39,710
  • 8
  • 110
  • 167
  • Very insightful. I understand that MMA typically "fails late", but indeed it would be very useful to have Mathematica complain earlier in this case. Thank you anyway! – MarcoB Jun 10 '15 at 06:18
1

This is in response to Chris's request in comments, but not an answer to his question.

I tried running the code you asked:

{ev, if, mesh} = helmholzSolve3D[\[CapitalOmega], 4, MaxCellMeasure -> 0.25]

with the definition of Omega that you have in the current question, and the definition of helmholzSolve3D from https://mathematica.stackexchange.com/a/85552/1089.

It runs for a few seconds, then unfortunately it fails with the following:

Mathematica graphics

I am on version 10.1.0, on Windows 7 (64-bit). My machine has 8GB of total RAM. Let me know if you want to try something different!

Community
  • 1
MarcoB
  • 67,153
  • 18
  • 91
  • 189
  • I tried with version 10.1 which I installed; it runs; but then crashed my computer by requesting too much memory (16Gb). It seems ToElementMesh produces too many simplices. – chris Jun 10 '15 at 09:04