RegionUnion Fails for Overlapping 3D Regions

Question

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["NDSolve`FEM`"]
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:

Any help gratefully received. Thanks.

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 '20 at 17:23
@Tim Laska. This looks interesting. I will have a play with it. Thanks for your help! — Paul Harrison, Apr 18 '20 at 18:51
WARNING: RegionUnion just crashed my kernel with these meshes. — Szabolcs, Apr 18 '20 at 19:33
I have trouble understanding what your final objective is. Cloud you elaborate a bit on that? — user21, Apr 20 '20 at 08: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 '20 at 08:46
@user21 Ps. The reason the bevelled corners were a problem is that (as I understand it) BCs, where unspecified, are inserted as Neumann BCs. My BCs are (obviously) only specified at the main faces of the box (z==0, x=-rng etc.). The little bevelled edges effectively form 12 additional faces, which by default get their own Neumann BCs, and distort the fields in the corners. I was concerned about that and thought the best way to deal with it was to follow the (your) prescription at the SE post I linked to my original question above, which was also about similarly cut-off corners (albeit in 2D). — Paul Harrison, Apr 20 '20 at 08:58
@PaulHarrison, I still do not understand the problem. Why do you not just use Cubiod or ImplicitRegion[True,...] to represent your region? — user21, Apr 20 '20 at 09:37
@user21 Thanks for the tips. I can try these. Cuboid, I had considered, but was unaware initially that if I specify explicitly a box-shaped volume, that the mesh creation would soften the corners, but not do the same to a Cuboid (which is a box-shaped region). Why would it not do the same to Cuboid? I also do not understand what putting True as an argument to ImplicitRegion does. I guess I can find it in the documentation, somewhere. — Paul Harrison, Apr 20 '20 at 09:48
@PaulHarrison, the fact that ImplicitRegion[z>=0,...] bevels the box that you specify is not wanted but a problem with the generality of the algorithm creating the mesh. Cuboid gives additional information (that the region is a cuboid, something with ImplicitRegion[z>=0,..] does not give) and a better algorithm can be used for the Cuboid case. — user21, Apr 20 '20 at 09:52
@user21 Thanks, now I understand! [At least part of "the problem" is that the user is not an expert :-) ]. Perhaps this will also help with my wider problems. Will try it now. — Paul Harrison, Apr 20 '20 at 09:55

score 4 · Answer 1 · answered Apr 18 '20 at 18:40

As mentioned in the comment, you could try this:

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

score 0 · Answer 2 · answered Apr 20 '20 at 09:57

0

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

Needs["NDSolve`FEM`"]
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.

answered Apr 20 '20 at 09:57

user21

39,710
8
110
167

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 '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 '20 at 11:16

score -1 · Answer 3 · answered Apr 18 '20 at 19:43

Needs["NDSolve`FEM`"]
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 '20 at 09:02

RegionUnion Fails for Overlapping 3D Regions

3 Answers3