# RegionUnion Fails for Overlapping 3D Regions

Bug introduced in 12.0 or earlier and persisting through 12.1.0

Warning: The below code may crash your kernel!

I am trying to make a (simple) 3D mesh of a box for a FEM solution (I am using MMA 12.1 on a Mac). By default, the mesh "bevels" the edges of my box. So, I try to follow this, at least as closely as possible: ElementMesh from ImplicitRegion cuts corners of region The idea is to mesh the box edges separately (once each) and merge all of them into a single box with much sharper corners. But the first merge (RegionUnion) fails. Here is my simplified code which demonstrates the issue:

Needs["NDSolveFEM"]
rng = 10.;
solnRegn =
ImplicitRegion[z >= 0, {{x, -rng, rng}, {y, -rng, rng}, {z, 0,
rng}}];
mr0 = MeshRegion[ToElementMesh[solnRegn, "MeshOrder" -> 1]];
mesh = ToElementMesh[mr0];
Print[Magnify[mesh["Wireframe"], 1.5]];
Clear[mesh]; edge1 =
ImplicitRegion[z > x + 2 rng - 1, {{x, -rng, -rng + 1}, {y, -rng, rng}, {z, rng - 1,
rng}}];
mr1 = MeshRegion[ToElementMesh[edge1, "MeshOrder" -> 1]];
mesh = ToElementMesh[mr1]; Print[
Magnify[mesh["Wireframe"], 1.5]]; Clear[mesh];
reg = RegionUnion[mr0, mr1];
mesh = ToElementMesh[reg];
Print[Magnify[mesh["Wireframe"], 1.5]];


Here is my output:

• I think that you should consider using the new OpenCascadeLink. OpenCascade is an open source cad package and it does a good job preserving sharp edges on boolean operations. I showed a couple of usage examples in my answer here. – Tim Laska Apr 18 at 17:23
• @Tim Laska. This looks interesting. I will have a play with it. Thanks for your help! – Paul Harrison Apr 18 at 18:51
• WARNING: RegionUnion just crashed my kernel with these meshes. – Szabolcs Apr 18 at 19:33
• I have trouble understanding what your final objective is. Cloud you elaborate a bit on that? – user21 Apr 20 at 8:23
• @user21 I am solving a Poisson-type equation (Laplacian and load term) in 3D in the box with Dirichlet and (implicit) Neumann BCs. Hope that helps? In fact, I have found a work-around now, since I find that I can refine the corners of the box using a smaller MaxBoundaryCellMeasure, while keeping my MaxCellMeasure quite large. Probably dumb that I did not try that sooner. The workaround does not, of course, solve the bug which Szabolcs has posted for this case. (I am having another problem now, which I will post as a separate question). – Paul Harrison Apr 20 at 8:46

As mentioned in the comment, you could try this:

Needs["OpenCascadeLink"]
Needs["NDSolveFEM"]
cube = OpenCascadeShape[Cuboid[{-rng, -rng, 0}, {rng, rng, rng}]]
(bmeshcube =
pp = Polygon[{{-rng + 1, -rng, rng}, {-rng, -rng, rng}, {-rng, -rng,
rng - 1}}];
axis = {{-rng, -rng, rng}, {-rng, rng, rng}};
Show[Graphics3D[{{Red, pp}, {Blue, Thick, Arrow[axis]}}],
bmeshprism["Wireframe"], Boxed -> False]
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]


I think the issue at hand might be a misunderstanding and that something like:

Needs["NDSolveFEM"]
rng = 10;
mesh = ToElementMesh[Cuboid[{-rng, -rng, 0}, {rng, rng, rng}](*,
"MeshElementType"\[Rule]TetrahedronElement*)];
mesh["Wireframe"]


will make the issue of the unwanted beveled edges go away.

• Thanks. My next step will be to replace the bottom of the box with a function ie. to require z>=f[x,y]. Is this still going to work with that? – Paul Harrison Apr 20 at 10:11
• @PaulHarrison, you should clarify this in your question. And no, you will need a different approach. Probably based on the OpenCascadeLink. – user21 Apr 20 at 11:16
Needs["NDSolveFEM"]
rng = 10.;
solnRegn =
ImplicitRegion[-rng <= x <= rng || -rng <= y <= rng ||
0 <= z <= rng, {x, y, z}];
mr0 = MeshRegion[ToElementMesh[solnRegn, "MeshOrder" -> 1]]
edge1 = ImplicitRegion[{z > x + 2 rng - 1 && z < rng &&
x > -rng}, {{x, -rng, -rng + 1}, {y, -rng, rng}, {z, rng - 1,
rng}}];
mr1 = MeshRegion[ToElementMesh[edge1, "MeshOrder" -> 1]]
mesh = ToElementMesh[
mr1];(*Print[Magnify[mesh["Wireframe"],1.5]];*)
(*Clear[mesh];*)
reg = RegionUnion[mr0, mr1];
boolreg = BooleanRegion[Or, {mr0, mr1}]
Region[boolreg]
ToElementMesh[boolreg]


The question includes code. That code is not a Mathematica style very much. Print, Magnify and alike are not of interest to the tasks.

mesh = ToElementMesh[mr0]
mesh1 = ToElementMesh[reg]
`

are obsolete too.

I did reconstruct some information into a simple functioning input for Mathematica.

• Thanks for your feedback, but with all due respect, I think that stating that a given idiom is "obsolete" is not really helping. The commands are still part of the language, are they not? It is a problem when they don't work as advertised. Who is dictating what is an obsolete idiom and what is not? All languages support modern and older idioms (after all, we learn a language for life). You also said you were able to reconstruct a simple functioning input, but you did not post it. You could improve your answer by so doing. – Paul Harrison Apr 20 at 9:02