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:

rng = 10.;
solnRegn = 
  ImplicitRegion[z >= 0, {{x, -rng, rng}, {y, -rng, rng}, {z, 0, 
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, 
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.

  • 2
    $\begingroup$ 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. $\endgroup$
    – Tim Laska
    Apr 18, 2020 at 17:23
  • $\begingroup$ @Tim Laska. This looks interesting. I will have a play with it. Thanks for your help! $\endgroup$ Apr 18, 2020 at 18:51
  • 2
    $\begingroup$ WARNING: RegionUnion just crashed my kernel with these meshes. $\endgroup$
    – Szabolcs
    Apr 18, 2020 at 19:33
  • $\begingroup$ I have trouble understanding what your final objective is. Cloud you elaborate a bit on that? $\endgroup$
    – user21
    Apr 20, 2020 at 8:23
  • $\begingroup$ @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). $\endgroup$ Apr 20, 2020 at 8:46

3 Answers 3


As mentioned in the comment, you could try this:

cube = OpenCascadeShape[Cuboid[{-rng, -rng, 0}, {rng, rng, rng}]]
(bmeshcube = 
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]]

OpenCascade Solution


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

rng = 10;
mesh = ToElementMesh[Cuboid[{-rng, -rng, 0}, {rng, rng, rng}](*,

enter image description here

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

  • $\begingroup$ 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? $\endgroup$ Apr 20, 2020 at 10:11
  • $\begingroup$ @PaulHarrison, you should clarify this in your question. And no, you will need a different approach. Probably based on the OpenCascadeLink. $\endgroup$
    – user21
    Apr 20, 2020 at 11:16
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, 
mr1 = MeshRegion[ToElementMesh[edge1, "MeshOrder" -> 1]]
mesh = ToElementMesh[
reg = RegionUnion[mr0, mr1];
boolreg = BooleanRegion[Or, {mr0, mr1}]

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.


  • 2
    $\begingroup$ 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. $\endgroup$ Apr 20, 2020 at 9:02

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.