Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question

2 Answers 2

up vote 13 down vote accepted

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.

share|improve this answer
    
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. –  Stan Sep 2 at 13:07
    
@Stan, please see my update. –  RunnyKine Sep 2 at 13:18
    
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. –  Stan Sep 2 at 13:32
    
@Stan, see my new edit. Hopefully this is what you want. –  RunnyKine Sep 2 at 13:47
    
@Stan The Method -> "Continuation" option apparently helps with corners and cusps but it doesn't do anything in this case. –  RunnyKine Sep 2 at 13:54

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.

share|improve this answer
1  
I was hoping you would chime in. +1 –  RunnyKine Sep 2 at 14:13
    
@RunnyKine, I was about not too :-) your answer is pretty good :-) –  user21 Sep 2 at 14:18
    
@user21: this is a very nice addition! Thanks! –  Stan Sep 2 at 14:28
    
I'm glad you did :) I like the analysis. –  RunnyKine Sep 2 at 14:32

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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