BoundaryMeshRegion Error: Boundary surface is not closed

Question

I am trying to create a BoundaryMeshRegion comprised of many cubes with random locations, in a relatively small region.

I have had success implementing my technique for a smaller number of cubes (say, 30), but have run into an error which sometimes pops up after repeated runs of the same number. With an increase in the number of cubes, this error frequently pops up, of the form (ignore the blah blah blah...that is just to eliminate a ton of numbers):

BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges Line[{{{blah},{blah}, {blah}}}] only come from a single face.

I have tried to fix this using the suggestion in this thread with very limited success: Understanding BoundaryMeshRegion error

As for the exact method in which these cubes are created (inside of a Do loop and appending the results out one by one, etc., etc.), I need it to be done this way, as this is the technique which must be used as a part in a separate, much more complicated simulation.

Is there anything that can be done to eliminate this problem, such that any number of cubes can be used (whilst maintaining the spirit of this block-building technique)? Best of luck!

P.S., I highly recommend running the same number of cubes repeatedly. I thought I had fixed it, but after another few simulations, the error popped up again.

Here is the code.

Clear["Global`*"];

cubeNumber = 40;
cubeTable = {};
pos0 = Table[{RandomReal[{0.5, 0.7}], RandomReal[{0.5, 0.7}], 
   RandomReal[{0, 0.1}]}, cubeNumber]

Do[
  cubes = 
   Graphics3D[
    Cube[{pos0[[k]][[1]], pos0[[k]][[2]], pos0[[k]][[3]]}, 0.05]];
  
  AppendTo[cubeTable, cubes]; 
  
  , {k, 1, cubeNumber}];

cubicMesh = 
 BoundaryDiscretizeRegion[
  RegionUnion[
   Flatten[Table[
     DiscretizeGraphics[cubeTable[[k]]], {k, Length[cubeTable]}]]]]

Answer 1

At first we use the OP's codes but remove Graphics3D and Part.

Clear["Global`*"];
cubeNumber = 40;
cubeTable = {};
pos0 = Table[{RandomReal[{0.5, 0.7}], RandomReal[{0.5, 0.7}], 
   RandomReal[{0, 0.1}]}, cubeNumber];
Do[cubes = Cube[pos0[[k]], 0.05];
  AppendTo[cubeTable, cubes];, {k, 1, cubeNumber}];

After that we again recommend to use NDSolveFEM Or OpenCascadeLink to make the BoundaryMeshRegion effective.(Compare with BoundaryDiscretizeRegion[regs])

Needs["NDSolve`FEM`"];
regs = cubeTable // RegionUnion;
bm = ToBoundaryMesh[regs];
bmr = BoundaryMeshRegion[bm]
bmr // Volume
RegionBoundary[bmr] // Area
bm["Wireframe"]

regs = cubeTable // RegionUnion;
Needs["OpenCascadeLink`"];
(* Needs["NDSolve`FEM`"]; *)
shape = OpenCascadeShape[regs];
bm = OpenCascadeShapeSurfaceMeshToBoundaryMesh[shape];
bmr = BoundaryMeshRegion[bm]
bmr // Volume
bm["Wireframe"]

0.00348886.

Answer 2

I know that you said that your approach is dictated by other limitations, but there are still quite a few things you could improve in your code. Here is refactored example:

SeedRandom[1];

cubeNumber = 40;

pos0 = 
  RandomVariate[
    UniformDistribution[{{0.5, 0.7}, {0.5, 0.7}, {0, 0.1}}],
    cubeNumber
  ];

cubeTable = Cube[#, 0.05]& /@ pos0;

Region[RegionUnion @@ cubeTable]

At least you could consider:

1. Generating the positions in one call to RandomVariate rather than a table.
2. Mapping Cube over the list of positions rather than using a table. Note also that your Cube[{pos0[[k]][[1]], pos0[[k]][[2]], pos0[[k]][[3]]}, 0.05] is completely equivalent to Cube[pos0[[k]], 0.05].
3. Consider avoiding the discretization of graphics and perhaps use regions instead.