FEM: Nicer Element Shape for Spherical Region

Question

I'm trying to generate a mesh for later use in the Finite Element Method of the DSolve command. It is basically a parallelepiped with a spherical indentation. I'm trying to generate the mesh as follows:

Needs["NDSolve`FEM`"]
region = ImplicitRegion[!(Norm[{x, y, z}] < 1), {{x, -5, 5}, {y, -5, 5}, {z, 0, 5}}];
mesh = ToElementMesh[region];
Show[mesh["Wireframe"], PlotRange -> {{-2, 2}, {-2, 2}, {0, 1}}]

However, as you an see on the picture below this produces a really irregular and inaccurate (as I think) mesh for the spherical part. Is there a way to improve it and to generate more and homogeneously distributed elements on the spherical part?

Answer 1

The mesh seems to be fine and you can see that it is by doing:

region = ImplicitRegion[! (Norm[{x, y, z}] < 1), {{x, -5, 5}, {y, -5, 5}, {z, 0, 5}}];

m = DiscretizeRegion[region, {{-2, 2}, {-2, 2}, {0, 1}}]

To view as wireframe you can do:

Needs["NDSolve`FEM`"]

mesh = ToElementMesh[m] // Quiet;

Then:

Show[mesh["Wireframe"]]

If you want to keep your original region, you can vary MaxCellMeasure in DiscretizeRegion (of course, this will increase computation time and the number of mesh elements)

m2 = DiscretizeRegion[region, MaxCellMeasure -> 0.001]

Show[m2, PlotRange -> {{-2, 2}, {-2, 2}, {0, 1}}]

Edit

It seems there's a way to do what you want and hopefully keep performance hit to a minimum using MeshRefinementFunction:

m3 = DiscretizeRegion[region, MeshRefinementFunction -> Function[{vertices, vol}, 
     Block[{x, y, z}, {x, y, z} = Mean[vertices]; 
       If[-1 < x < 1 && -1 < y < 1 && 0 < z < 1, vol > 0.001, vol > 0.1]]]]

Then:

mesh3 = ToElementMesh[m3] // Quiet;

Show[mesh3["Wireframe"], PlotRange -> {{-2, 2}, {-2, 2}, {0, 1}}]

Looks reasonable.

Answer 2

Here are a few additions to @RunnyKine suggestions. If you are ever in doubt about the quality of a mesh (an ElementMesh to be exact) you can query the mesh.

Needs["NDSolve`FEM`"]
region = ImplicitRegion[! (Norm[{x, y, z}] < 1), {{x, -5, 5}, {y, -5, 
     5}, {z, 0, 5}}];
mesh = ToElementMesh[region];

Min[mesh["Quality"]]
0.004439742441262357`

So the minimum mesh quality does not look too good. The overall mesh quality distribution seems sort of OK.

Histogram[mesh["Quality"]]

We can do better by increasing the sample points of the underlying RegionPlot (which is the only method available in 3D) with:

mesh = ToElementMesh[region, 
   "BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 25}];
Show[mesh["Wireframe"], PlotRange -> {{-2, 2}, {-2, 2}, {0, 1}}]

Looking at the quality again:

Histogram[mesh["Quality"]]

More information can be found in the documentation for ToElementMesh and in the mesh generation tutorial.