I'm trying to create a finite element mesh with some spherical holes.

To my mind the behavior is erratic but generally fails. I'm using MMA 10.0.2. In the following 3 different sets of holes are defined by coords3Da etc. By setting coords3D to coords3Da (or coords3Db or coords3Dc) and running the code you will see that it generates a mesh for b and c but fails for a. I haven't been able to generate a mesh for any 3-hole version. I also tried a second version in which I used MMA's defined functions "Cylinder" and "Ball." That attempt fails to generate a mesh with a single hole. (It does generate a mesh if the RegionHoles option is removed.)

Is there any way to tweak one of these to get it to generate meshes for a few tens of holes? I should mention that this code was morphed from the code here: How to create subregions for the NDSolve FEM Solver


I've since succeeded in generating meshes with as many as six holes but it's all hit and miss. I'm after a robust way to generate meshes.

<< NDSolve`FEM`;

cylinder[x_, y_, z_, h_, r_] := x^2 + y^2 <= r^2 && 0 <= z <= h;
holes3D[{x0_, y0_, z0_}, r_] := ((x - x0)^2 + (y - y0)^2 + (z - z0)^2 >= r^2);

coords3Da = {{-0.5, 0, .5}, {0, 0.5, 3.5}};
coords3Db = {{0, 0, 2.5}, {1/2, 0, 1.5}};
coords3Dc = {{0, 0, 2.5}}

coords3D = coords3Db;

sd3D = And @@ (holes3D[#, 3/25] & /@ coords3D);
regions = And[cylinder[x, y, z, 5, 1], sd3D];

Ω2cyl = ImplicitRegion[regions, {x, y, z}]
mesh = ToElementMesh[Ω2cyl, "BoundaryMeshGenerator" ->{"Continuation"},"RegionHoles" -> coords3D];

Mathematica graphics


ToElementMesh[RegionUnion[Cylinder[{{0, 0, 0}, {0, 0, 5}}], Ball[{0, 0, 2.5}, 1/8]], "BoundaryMeshGenerator" -> {"Continuation"}, "RegionHoles" -> {0, 0, 2.5}]["Wireframe"]
  • $\begingroup$ For your second updated example, I do not think you want a RegionUnion here; at least that's not the same as the examples above. $\endgroup$
    – user21
    Commented Nov 7, 2016 at 13:52

2 Answers 2



This now behaves much better:

mesh = ToElementMesh[\[CapitalOmega]2cyl, 
   "BoundaryMeshGenerator" -> {"BoundaryDiscretizeRegion", 
     Method -> {"MarchingCubes", PlotPoints -> 33}}, 
   "RegionHoles" -> coords3D];

enter image description here

Update 2

In Version 11.1 we have way to generate high fidelity second order meshes from individual boundary meshes as documented in the section on Numerical Regions. This has the advantage that regions where we have a large base and small components can be meshed well. Here is an example:

We create the primitives and the region:

primitives = {Cylinder[{{0, 0, 0}, {0, 0, 5}}], 
   Ball[#, 1/8] & /@ coords3D};
rd = Fold[RegionDifference, primitives[[1]], primitives[[2]]];

With ToNumericalRegion we create a NumericalRegion and boundary element meshes for each of the primitives.

nr = ToNumericalRegion[rd];
bms = ToBoundaryMesh[#] & /@ Flatten[primitives];

Note, that in principal you have all the freedom you want to generate each of these boundary meshes.

Get the coordinates for each boundary mesh and determine the offset that we need for the boundary mesh elements to point to the correct coordinates:

coords = #["Coordinates"] & /@ bms;
coordOffsets = Most[FoldList[Plus, 0, Length /@ coords]];

Get the boundary mesh elements and generate a new combined boundary mesh.

beles = #["BoundaryElements"] & /@ bms;
bmesh = ToBoundaryMesh["Coordinates" -> Join @@ coords, 
   "BoundaryElements" -> 
    Flatten[MapIndexed[(index = #2[[1]]; 
        MapThread[#1[#2] &, {Head /@ beles[[index]], 
          coordOffsets[[index]] + 
           ElementIncidents[beles[[index]]]}]) &, beles]], 
   "RegionHoles" -> coords3D];

Now, we have a boundary mesh representation of the region rd. That's not really new. What you can do now, however, is helpful - you can associate the new boundary mesh with the numerical region.

SetNumericalRegionElementMesh[nr, bmesh]

When you now call ToElementMesh on the numerical region, a second order mesh can be generated. Also this second order mesh will have it's mid-side nodes moved to the correct position because the numerical region has a symbolic representation of the region.

mesh = ToElementMesh[nr];

enter image description here


0 <= x$239/5 <= 1 && x$237^2 + x$238^2 <= 1 && x$237^2 + x$238^2 + (-2.5 + x$239)^2 > 
  1/64 && (-1/2 + x$237)^2 + x$238^2 + (-1.5 + x$239)^2 > 1/64

{x$237, x$238, x$239}


This is probaly not an ideal solution, but perhaps a start. In 3D the only boundary mesh generator is (V10.1) the "RegionPlot" one. One thing that can be done is specify the "SamplePoints" in the x,y and z direction independently:

bmesh = ToBoundaryMesh[\[CapitalOmega]2cyl, 
   "BoundaryMeshGenerator" -> {"RegionPlot", 
     "SamplePoints" -> {17, 17, 64}}, "RegionHoles" -> coords3D];

I'd then generate the mesh from the boundary mesh with


Or, probably better, pass the option to ToElementMesh directly.


I highly suggest you skip the ImplicitRegion definition and to build all the boudary mesh elements "by hand" so to do not let the ToElementMesh function build the boundary mesh elements for you. Here's my sketched workflow:

1) create a single spherical boundary mesh and extrapolate the coordinates "coords" and the incidents "incs" of all the boundary mesh elements

{coords, incs} = 

ToBoundaryMesh[Sphere[]][#] & /@ {"Coordinates", 
 "BoundaryElements"} // {#[[1]], #[[2, 1, 1]]} &;

2) Create the list of all the coordinates by traslating coords to the points you want the holes to be. Lets say your holes are centered in the points given by ctrs then

Table[# + ctrs[[i]] & /@ coords, {i, 1, Length@ctrs}] // Join @@ # &;

3) Create the list of all the incidents; just "append" incs n times to itself, where n is the number of your holes, like this

Table[# + n Max@# &@incs, {n, 0, Length@ctrs - 1}] // Join @@ # &;

5) now you can built a ToBoundaryMesh object made of TriangleElement with the two previous lists and you must append the coordinates and incidents of the mesh elements relative to the external cylinder. Then pass it to ToElementMesh specifying the RegionHoles coordinates. I suggest to set the MaxCellBoudary to infinity when ToElementMesh because this won't affect the boundary mesh elements of the holes you manually defined

This method works flawlessly for thousands of holes.

  • $\begingroup$ (+1) While this is a way to go, the current method will only produce first order accurate meshes. One would need to add second order nodes, with e.g. MeshOrderAlteration[mesh, 2]. Then the last thing to be done is move the mid side notes onto the boundary. $\endgroup$
    – user21
    Commented Nov 7, 2016 at 11:04
  • $\begingroup$ Could you add an example for ctrs? $\endgroup$
    – user21
    Commented Nov 7, 2016 at 11:05
  • $\begingroup$ ctrs are the list you called coords3Da, coords3Db and so on $\endgroup$
    – Fortsaint
    Commented Nov 8, 2016 at 2:04
  • $\begingroup$ Have a look at the update 2 of my answer. That should helps this type of approach. $\endgroup$
    – user21
    Commented Mar 28, 2017 at 14:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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