Is there a better way to produce a shape with exclusions in RegionPlot3D? [duplicate]

Question

I'm attempting to create a certain shape. It's essentially a rectangular prism with an array of cylinders cut out of it, and then some half shells of cylinders along the sides. My attempt is basically this:

1. Create a list of these cylinders, each defined by ImplicitRegion, then combine them into a single region with RegionUnion.
2. Define the rectangular prism, also with ImplicitRegion
3. Define the total region with ImplicitRegion by including the points that are a part of the prism, but not the cylinder array
4. Plot this resulting region with RegionPlot3D

Here's the code for that, just plotting the cylinder array to begin:

 hCub = 1.5;
W = 5;
L = 8;
CYLrad = W/20;
shellthickness = CYLrad/6;
deltah = hCub*0;
pitch = 1;
xnum = 5;
ynum = 6;
xinnercylnum = xnum - 2;
xpitch = W/(xnum - 1);
ypitch = xpitch*Sqrt[3.]/2;
Clear[x, y, z];
firstcylpos = {xpitch, ypitch, 0};
cyltable = 
  RegionUnion@Flatten[#, 1] &@
   Table[Table[
     ImplicitRegion[
      0 <= (x - (j - 1)*xpitch)^2 + (y - i*ypitch)^2 <= CYLrad^2 && 
       0 <= (z - 0) <= hCub + deltah(*&&
      x\[LessEqual]cylcenter[[1]]*), {x, y, z}], {i, 0, ynum, 1}], {j,
      xnum}];
innercyltable = 
  RegionUnion@Flatten[#, 1] &@
   Table[Table[
     ImplicitRegion[
      0 <= (x - j*xpitch)^2 + (y - i*ypitch)^2 <= CYLrad^2 && 
       0 <= (z - 0) <= hCub + deltah(*&&
      x\[LessEqual]cylcenter[[1]]*), {x, y, z}], {i, ynum}], {j, 
     xinnercylnum}];
yzface1cylshelltable = 
  RegionUnion@
   Table[ImplicitRegion[(CYLrad - 
          shellthickness)^2 <= (x - W)^2 + (y - i*ypitch)^2 <= 
       CYLrad^2 && 0 <= (z - 0) <= hCub + deltah && x <= W, {x, y, 
      z}], {i, ynum}];
yzface2cylshelltable = 
  RegionUnion@
   Table[ImplicitRegion[(CYLrad - 
          shellthickness)^2 <= (x - 0)^2 + (y - i*ypitch)^2 <= 
       CYLrad^2 && 0 <= (z - 0) <= hCub + deltah && x >= 0, {x, y, 
      z}], {i, ynum}];
xzface2cylshelltable = 
  RegionUnion@
   Table[ImplicitRegion[(CYLrad - 
          shellthickness)^2 <= (x - i*xpitch)^2 + (y - 0)^2 <= 
       CYLrad^2 && 0 <= (z - 0) <= hCub + deltah && y >= 0, {x, y, 
      z}], {i, xinnercylnum}];
allcylshells = 
  RegionUnion[yzface1cylshelltable, yzface2cylshelltable, 
   xzface2cylshelltable];
cub = ImplicitRegion[
   0 <= x <= W && 0 <= y <= L && 0 <= z <= hCub, {x, y, z}];
cubnocyls = RegionDifference[cub, cyltable];
Timing@RegionPlot3D[{Evaluate@cubnocyls, Evaluate@innercyltable, 
   Evaluate@allcylshells}, 
  PlotStyle -> {{Black, Lighting -> "Neutral", 
     Opacity@1}, {Lighting -> "Neutral", 
     White}, {Lighting -> "Neutral", White}}, PlotPoints -> 20, 
  Boxed -> False, ImageSize -> 800]

Doing this takes about 12 second and produces:

This isn't an especially complicated structure (relatively), and it's at a much lower resolution than I'll eventually need to do (PlotPoints->20). I cranked it up to PlotPoints->100, it took ~11 minutes, and it still looked pretty bad.

I'm wondering if I'm just doing this in a very naive way, and there's a more efficient way to do it. Is there?

Answer 1

We can work in 2D and use RegionProduct to extrude into 3D. I've modified your code slightly:

hCub = 1.5;
W = 5;
L = 8;
CYLrad = W/20;
shellthickness = CYLrad/6;
deltah = hCub*0;
pitch = 1;
xnum = 5;
ynum = 6;
xinnercylnum = xnum - 2;
xpitch = W/(xnum - 1);
ypitch = xpitch*Sqrt[3.]/2;
firstcylpos = {xpitch, ypitch, 0};

mr = DiscretizeRegion[ImplicitRegion[
  0 <= x^2 + y^2 <= CYLrad^2, {x, y}]];
cyltable = 
  RegionUnion@Flatten[#, 1] &@
   Table[Table[TransformedRegion[mr, 
   TranslationTransform[{(j - 1)*xpitch, i*ypitch}]], {i, 0, ynum, 1}], {j,xnum}];
innercyltable = 
  RegionUnion@Flatten[#, 1] &@
   Table[Table[
     TransformedRegion[mr, 
   TranslationTransform[{j*xpitch, i*ypitch}]], {i, ynum}], {j, 
     xinnercylnum}];

mr = DiscretizeRegion[ImplicitRegion[(CYLrad - 
          shellthickness)^2 <= (x - W)^2 + y^2 <= 
       CYLrad^2 && x <= W, {x, y}]];
    yzface1cylshelltable = 
  RegionUnion@
   Table[TransformedRegion[mr, 
   TranslationTransform[{0, i*ypitch}]], {i, ynum}];

mr = DiscretizeRegion[ImplicitRegion[(CYLrad - 
          shellthickness)^2 <= x^2 + y^2 <= 
       CYLrad^2 && x >= 0, {x, y}]];
yzface2cylshelltable = 
  RegionUnion@
   Table[TransformedRegion[mr, 
   TranslationTransform[{0, i*ypitch}]], {i, ynum}];

mr = DiscretizeRegion[ImplicitRegion[(CYLrad - 
          shellthickness)^2 <= x^2 + y^2 <= 
       CYLrad^2 && y >= 0, {x, y}]];
xzface2cylshelltable = 
  RegionUnion@
   Table[TransformedRegion[mr, 
   TranslationTransform[{i*xpitch, 0}]], {i, xinnercylnum}];

allcylshells = 
  RegionUnion[yzface1cylshelltable, yzface2cylshelltable, 
   xzface2cylshelltable];

cub = DiscretizeGraphics[Rectangle[{0, 0}, {W, L}]];
cubnocyls = RegionDifference[cub, cyltable];

extruded = RegionProduct[#, Line[{{0}, {hCub + deltah}}]]& /@ 
  {cubnocyls, innercyltable, allcylshells};

styles = {{Black, Lighting -> "Neutral", 
     Opacity@1}, {Lighting -> "Neutral", 
     White}, {Lighting -> "Neutral", White}};

Show[MapThread[MeshRegion[#1, BaseStyle -> #2] &, {extruded, styles}]]