3D FEM with holes

Question

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

EDITED TO ADD:

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];
Show[mesh["Wireframe"]]

SECOND VERSION USING Cylinder and Ball

ToElementMesh[RegionUnion[Cylinder[{{0, 0, 0}, {0, 0, 5}}], Ball[{0, 0, 2.5}, 1/8]], "BoundaryMeshGenerator" -> {"Continuation"}, "RegionHoles" -> {0, 0, 2.5}]["Wireframe"]

Answer 1

Update

This now behaves much better:

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

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];
mesh["Wireframe"]

mesh["MeshOrder"]
2

nr["Predicates"]
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

nr["PredicateVariables"]
{x$237, x$238, x$239}

Old

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];
Show[bmesh["Wireframe"]]

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

ToElementMesh[bmesh]

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

Answer 2

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.