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.

Answer 1

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]]

Answer 2

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.

Answer 3

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.