11
$\begingroup$

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?

enter image description here

Improve this question
$\endgroup$

2 Answers 2

Reset to default
14
$\begingroup$

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}}]

Mathematica graphics

To view as wireframe you can do:

Needs["NDSolve`FEM`"] mesh = ToElementMesh[m] // Quiet;

Then:

Show[mesh["Wireframe"]]

Mathematica graphics

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]

Mathematica graphics

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}}]

Mathematica graphics

Looks reasonable.

Improve this answer
$\endgroup$
7
  • $\begingroup$ I'm not sure if it is a display problem or a discretization problem. You have chosen a smaller region, and the sphere part looks nice. But whenever my sphere is small in comparison to the big parallelepiped, e.g. sphere of radius 1 in a parallelepiped of dimensions (-5,5)x(-5,5)x(0,5), then it starts to look ugly as in my picture. $\endgroup$ Commented Sep 2, 2014 at 13:07
  • $\begingroup$ @Stan, please see my update. $\endgroup$ Commented Sep 2, 2014 at 13:18
  • $\begingroup$ yes this helps, but with the impact on the performance that you mentioned. I guess there is no easy way to generate a fine-grained mesh around the sphere and a coarse-grained one farther away? Something like to pre-generate a fine-grained mesh on the surface of the sphere separately, and then use it to generate the rest. Or to generate the parallelepiped and the sphere separately, and then to "subtract" the sphere from the block. Not sure if this makes sense. $\endgroup$ Commented Sep 2, 2014 at 13:32
  • $\begingroup$ @Stan, see my new edit. Hopefully this is what you want. $\endgroup$ Commented Sep 2, 2014 at 13:47
  • $\begingroup$ @Stan The Method -> "Continuation" option apparently helps with corners and cusps but it doesn't do anything in this case. $\endgroup$ Commented Sep 2, 2014 at 13:54
13
$\begingroup$

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"]]

enter image description here

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}}]

enter image description here

Looking at the quality again:

Histogram[mesh["Quality"]]

enter image description here

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

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ I was hoping you would chime in. +1 $\endgroup$ Commented Sep 2, 2014 at 14:13
  • 1
    $\begingroup$ @RunnyKine, I was about not too :-) your answer is pretty good :-) $\endgroup$ Commented Sep 2, 2014 at 14:18
  • $\begingroup$ @user21: this is a very nice addition! Thanks! $\endgroup$ Commented Sep 2, 2014 at 14:28
  • $\begingroup$ I'm glad you did :) I like the analysis. $\endgroup$ Commented Sep 2, 2014 at 14:32

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.