BoundaryDiscretizeRegion given an ImplicitRegion never returns

I am trying to generate random points on a "simple" 3D surface for testing purposes. The surface is a plane with etched or embossed rectangular structures.

I am using ImplicitRegion to model it :

GenerateCubes[y_, d_] :=
Cuboid[{#1, y, Min[0, #2]}, {#1 + 10, y + 10, Max[0, #2]}] &,
{Range[10, 100 - 1, 20], d*2^Range[5]}]

reg = ImplicitRegion[
(RegionMember[Cuboid[{0, 0, -6.5}, {110, 90, 0.}], {x, y, z}] ||
Or @@ (RegionMember[#, {x, y, z}] & /@ GenerateCubes[10, 0.2]) ||
Or @@ (RegionMember[#, {x, y, z}] & /@ GenerateCubes[30, 0.1])) &&
And @@ (Not[RegionMember[#, {x, y, z}]] & /@ GenerateCubes[50, -0.1]) &&
And @@ (Not[RegionMember[#, {x, y, z}]] & /@ GenerateCubes[70, -0.2]),
{{x, 0, 110}, {y, 0, 90}, {z, -6.5, 6.5}}];


This is what the object looks like :

RegionPlot3D[reg,
PlotPoints -> 50, PlotRange -> {{0, 110}, {0, 90}, {-6.5, 6.5}}, Axes -> True]


Now when I try to discretize it (in order to randomly sample points in the region), the notebook freezes:

BoundaryDiscretizeRegion[reg]


What's happening? Is the region too complicated to be discretized?

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jul 16 '15 at 1:25
• You're right, I copied an old code, sorry about the mistake ! I edited the question. Jul 16 '15 at 1:55

The region is a union of cuboids, which in FEM-talk is a HexahedronElement. We can just take all the coordinate boundaries, divide space into cuboids according to the coordinates, and select the ones whose mid point lies in the region defined by the OP's expression, which I store in exp.

Needs["NDSolveFEM"];

exp = (RegionMember[Cuboid[{0, 0, -6.5}, {110, 90, 0.}], {x, y, z}] ||
Or @@ (RegionMember[#, {x, y, z}] & /@
GenerateCubes[10, 0.2]) ||
Or @@ (RegionMember[#, {x, y, z}] & /@
GenerateCubes[30, 0.1])) &&
And @@ (Not[RegionMember[#, {x, y, z}]] & /@
GenerateCubes[50, -0.1]) &&
And @@ (Not[RegionMember[#, {x, y, z}]] & /@
GenerateCubes[70, -0.2]) // Rationalize;

coords = Union @@
Cases[exp, a_ <= # <= b_ :> {a, b}, Infinity] & /@ {x, y, z};

cuboids = Tuples[#][[{1, 5, 7, 3, 2, 6, 8, 4}]] & /@
Select[
Tuples[Partition[#, 2, 1] & /@ coords],
exp /. Thread[{x, y, z} -> Mean /@ #] &
];
ecoords = Union @@ cuboids;
nf = Nearest[ecoords -> Automatic];
emesh = ToElementMesh[
"Coordinates" -> ecoords,
"MeshElements" -> {HexahedronElement[Map[First@*nf, cuboids, {2}]]}
];

mr = MeshRegion[emesh];
RegionPlot3D@ mr


One advantage to direct construction is that the edges of the region are sharp, exact, not approximated.

• I must say, Michael, you have become quite a supporter of the FEM stuff; I like and apprechiate that ;-) Jul 16 '15 at 12:22
• Thanks. A couple of reasons are good functionality and good documentation that supports learning it. Also a few of your posts on the site have been helpful. :) Jul 16 '15 at 14:00
• I was looking for a solution that would not involve statically designing the region, more like the original question or the solution proposed by bbgodfrey , but it seems to be impossible at the moment. It is especially important as some other objects I work with might not be so easy to generate using FEM. I don't care for approximate regions using DiscretizeRegion as long as I have a general solution independent of the objects I am using, but I don't think it is currently possible to discretize regions with creases (negative portions), am I right ? Jul 16 '15 at 15:08
• @user1766172 If exp changes and the corresponding region is still the addition and subtraction of cuboids, then the above code should compute the corresponding mesh. I don't see how it's any more static than your original code. Now, if you wish to have other shapes, spheres, tetrahedra or rotated cuboids, then my code won't work on that. I think the problem with discretizing is the same as the problem with Simplify[exp] or even just LogicalExpand[exp]. The mix of many &&, ||, <= seems to cause a combinatorial explosion of cases. Jul 16 '15 at 15:20
• @MichaelE2 Yes I was speaking of other shapes, especially spheres and ellipsoids. But you summed it up perfectly, at the end of the day the region will be implicitly represented by a combination of too many conditions. I think I'll have to go for a numeric scheme if I want a generic solution rather than an analytic one using discretization. Thanks. Jul 16 '15 at 15:49

BoundaryDiscretizeRegion can handle a variant of your ImplicitRegion just fine.

reg = ImplicitRegion[
(RegionMember[Cuboid[{0, 0, -6.5}, {110, 90, 0.}], {x, y, z}] ||
Or @@ (RegionMember[#, {x, y, z}] & /@ GenerateCubes[10, 0.2]) ||
Or @@ (RegionMember[#, {x, y, z}] & /@ GenerateCubes[30, 0.1])) ||
Or @@ (RegionMember[#, {x, y, z}] & /@ GenerateCubes[50, 0.3]) ||
Or @@ (RegionMember[#, {x, y, z}] & /@ GenerateCubes[70, 0.4]),
{{x, 0, 110}, {y, 0, 90}, {z, -6.5, 6.5}}];


BoundaryDiscretizeRegion[reg]


On the other hand, BoundaryDiscretizeRegion for the code in the question runs a long time. Evidently, BoundaryDiscretizeRegion has difficulty coping with the subtractive ("etched") portion of the ImplicitPlot. I shall consider whether a work-around exists.

By the way, use the option MaxCellMeasure to improve the resolution of the BoundaryDiscretizeRegion

Solution

An alternative, more compact way to produce the plot in the question is

reg = BooleanRegion[Xor, Flatten[
{Cuboid[{0, 0, -6.5}, {110, 90, 0.}], GenerateCubes[10, 0.2],
GenerateCubes[30, 0.1], GenerateCubes[50, -0.1], GenerateCubes[70, -0.2]}]];


A RegionPlot3D of this expression is, of course, identical to that in the question. Since this plot is constructed from a surface mesh, we should be able to extract it. In fact, BoundaryDiscretizeGraphics is designed for that purpose. Unfortunately, it responds that it is "not implemented for" this mesh. However, it is easy to construct a do-it-yourself equivalent. With the plot set equal to p,

Quiet[Cases[p, GraphicsGroup[{z_}] -> z, Infinity] /. i_Integer -> p[[1, 1, i]]];


is a mesh representation of the surface. It is not so pretty, because the zoning is not regular, but it is a reasonably accurate mesh. Greater accuracy, if needed, can be obtained by increasing PlotPoints in the original plot.

Graphics3D[%]


Better Solution

Define the region with taller Cuboids, simply to see better the final result.

r = BooleanRegion[Xor, Flatten[
{Cuboid[{0, 0, -65}, {110, 90, 0}], GenerateCubes[10, 02], GenerateCubes[30, 01],
GenerateCubes[50, -01], GenerateCubes[70, -02]}]];


Next, plot the Mesh only for this object.

p = RegionPlot3D[r, PlotPoints -> 50, PlotRange -> {{0, 110}, {0, 90}, {-65, 65}},
Mesh -> 33, PlotStyle -> None]


I presume that only the top surface is desired. If so, it is easy to delete the sides and bottom.

DeleteCases[p // Normal, ({_, _, -65.} | {0., _, _} | {_, 90., _} | {110., _, _} |
{_, 0., _}), Infinity];


Finally, convert this into a region.

Quiet@DiscretizeGraphics[%]


This seems like a fairly general approach. (The specific numbers introduced into the third step are just the boundaries of the box.) Increase Mesh -> 33 to a larger value for greater resolution.

• You're right, this version is discretized correctly. However, I definitely need to take into account the subtractive portions.Do you mean to say that if I wait long enough, it should also work with the original implicit region (with subtractive region) ? Jul 16 '15 at 2:31
• Yes it works well for this exact object, but as I mentioned, I was looking for a more generic solution. BoundaryDiscretizeRegion was IMO the best way to achieve this generically, but it seems to only work for some specific cases... Is it normal behaviour, meaning that in future releases it might handle the original case much better ? Thanks Jul 16 '15 at 15:14
• @user1766172 My final solution is, I believe, quite general. Nonetheless, I wonder whether a mesh is even what you want. Instead, you could obtain a List of Polygons that make up a 3D plot of whatever you like, and distribute your random numbers first to List elements weighted by their Areas and then to positions within the selected Polygons. Jul 16 '15 at 17:20
• The solution is good too, I can only accept one of the answers though. I wanted to use a mesh -- in case the solution was generic enough -- since I was planning on using actual CAD objects imported as STL or PLY files. But it seems it won't be this easy if the objects are complex, as I won't be able to use any of the solutions proposed here generically for custom 3D objects. Jul 16 '15 at 19:02