# Automatically generating boundary meshes for region intersections

I am trying to solve a PDE having different coefficients in different regions. I have been having trouble getting Mathematica to do it automatically, so I have gone down the route of creating my own meshes. I would like to have high quality boundaries between regions with different coefficients. In the documentation it looks like this can be done by specifying coordinates, but I would like to find a way to automate this using regions. The code below shows something like what I want to do. Note, my final problem will be in 3D and will have complicated geometry so building the boundary regions from coordinates manually would be very difficult. The last line of the code is invalid, but shows kind of what I would like to do.

<< NDSolveFEM
a = ImplicitRegion[x >= 0 && y >= 0 && x <= 1 && y <= 1, {x, y}];
b = ImplicitRegion[x >= 1/4 && y >= 1 && x <= 3/4 && y <= 2, {x, y}];
o = RegionUnion[a, b];
m = ToElementMesh[o];
m["Wireframe"]
isec = RegionIntersection[a, b];
bm = ToBoundaryMesh[o, "BoundaryElements" -> {isec}] As an alternative, here is a small function that merges boundary element meshes and preservers the possibly internal boundaries.

ClearAll[BoundaryMeshMerge]
BoundaryMeshMerge[bmIn1_] /; BoundaryElementMeshQ[bmIn1] := bmIn1;

BoundaryMeshMerge[bmIn1_, bmIn2_] /;
BoundaryElementMeshQ[bmIn1] && BoundaryElementMeshQ[bmIn2] :=
Block[
{bm1 = bmIn1, bm2 = bmIn2, c1, c2, nc1, newC, newBEle, eleType,
inci1, inci2},
c1 = bm1["Coordinates"];
c2 = bm2["Coordinates"];
nc1 = Length[c1];
newC = Join[c1, c2];
eleType =
inci1 = ElementIncidents[bm1["BoundaryElements"]];
inci2 = ElementIncidents[bm2["BoundaryElements"]] + nc1;
newBEle = MapThread[#1[#2] &, {eleType, Join[inci1, inci2]}];
If[BoundaryIntersectionQ[newC, newBEle],
{newC, newBEle} = FixBoundaryIntersection[newC, newBEle];
];
ToBoundaryMesh["Coordinates" -> newC,
"BoundaryElements" -> newBEle]
]

BoundaryMeshMerge[bmIn1_?BoundaryElementMeshQ,
bmInRest__?BoundaryElementMeshQ] :=
Fold[BoundaryMeshMerge, bmIn1, {bmInRest}]


Here is how it works: Generate boundary meshed for each sub-region.

<< NDSolveFEM
a = ImplicitRegion[x >= 0 && y >= 0 && x <= 1 && y <= 1, {x, y}];
b = ImplicitRegion[x >= 1/4 && y >= 1 && x <= 3/4 && y <= 2, {x, y}];
bm1 = ToBoundaryMesh[a, "MaxBoundaryCellMeasure" -> Infinity];
bm2 = ToBoundaryMesh[b, "MaxBoundaryCellMeasure" -> Infinity];


Merge them:

(bmM = BoundaryMeshMerge[bm1, bm2])["Wireframe"] Note, that the internal boundary is maintained. Next you connect the boundary mesh just generated to the symbolic region description via a NumericalRegion. This can help for to get a better accuracy of the mesh for curved boundaries (see below).

nr = ToNumericalRegion[RegionUnion[a, b]];
SetNumericalRegionElementMesh[nr, bmM];
ToElementMesh[nr]["Wireframe"] One thing to be a bit careful about is that the automatic region hole finding may get a bit confused but you can overwrite that. Another thing to be aware of is to regenerate the numerical region nr every time before you call ToElementMesh as a NumericalRegion caches some results which may lead to unexpected results when calling ToElementMesh.

Consider this example:

a = Disk[];
b = Disk[{1/2, 0}];
b1 = ToBoundaryMesh[a];
b2 = ToBoundaryMesh[b];
bmM = BoundaryMeshMerge[b1, b2];

Total[ToElementMesh[bmM]["MeshElementMeasure"], 2] -
Integrate[1, {x, y} \[Element] RegionUnion[a, b]]
-0.010391591157441482


Accuracy is not too good. Now we use a NumericalRegion and help the region hole finding algorithm:

nr = ToNumericalRegion[RegionUnion[a, b]];
SetNumericalRegionElementMesh[nr, bmM];
mesh = ToElementMesh[nr, "RegionHoles" -> None];
(* mesh["Wireframe"] *)


This gives a better result:

Total[mesh["MeshElementMeasure"], 2] -
Integrate[1, {x, y} \[Element] RegionUnion[a, b]]
-3.1916053924163634*^-6


Hope this helps.

• Wow, thank you very much. That will help a lot. – tensor Sep 26 '17 at 14:30

I don't know if it help, but worth the shot. Is something like this what you want?

<< NDSolveFEM
a = ImplicitRegion[x >= 0 && y >= 0 && x <= 1 && y <= 1, {x, y}];
b = ImplicitRegion[x >= 1/4 && y >= 1 && x <= 3/4 && y <= 2, {x, y}];
o = RegionUnion[a, b];
m = ToElementMesh[o,
MeshRefinementFunction ->
Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
If[y <= 1.05 && y >= 0.95 && x <= 0.75 && x >= 0.25  ,
area > 0.0002, area > 1]]]];
m["Wireframe"]
isec = RegionIntersection[a, b] • Thank you. That's clever. It isn't quite what I wanted, but I can probably make it work for my purposes, especially if I can find a Mathematica function to give me the minimum distance from a point to a region. I was hoping there was something direct that solves this it would seem to be really valuable for finite elements. – tensor Sep 24 '17 at 20:15
• @tensor Have you seen RegionDistance? – C. E. Sep 25 '17 at 8:08