# Sometimes Region can't build right region

During my work with NDSolveFEM I found strange behavior of function Region.

a = 600;
reg1 = Rectangle[{0, 0}, {6 a, a}];
hole1 = Disk[{3 a, 0.75 a}, a/8];
hole2 = Disk[{3 a, 0.25 a}, a/8];
hole3 = Disk[{2.5 a, 0.5 a}, a/8];
resreg = RegionDifference[reg1, RegionUnion[hole1, hole2, hole3]]


So, it is really simple region with 3 holes near the centre of rectangle 6a * a.

resreg is:

BooleanRegion[#1 && ! (#2 || #3 || #4) &, {Rectangle[{0, 0},
{3600, 600}], Disk[{1800, 450.}, 75], Disk[{1800, 150.}, 75], Disk[{1500., 300.}, 75]}]


If I try to plot it by Region on my first computer, I will receive the simple Quad region 1 * 1:

So, this region is unmeshable and I can't understand the reason of this behavior.

If I try it on my notebook. I will receive the right plot and meshable region.

Mesh of region:

Both computers have Windows 10 with Mathematica v 11.3. How I can make Region plot right image?

This is a bug and you should report it to [email protected]. Here is why - simply because DiscretizeRegion has the same problem and does not give a message about it:

a = 600;
reg1 = Rectangle[{0, 0}, {6 a, a}];
hole1 = Disk[{3 a, 0.75 a}, a/8];
hole2 = Disk[{3 a, 0.25 a}, a/8];
hole3 = Disk[{2.5 a, 0.5 a}, a/8];
resreg = RegionDifference[reg1, RegionUnion[hole1, hole2, hole3]];


The finite element functions have the same problem, but they give a message about what is going on. Start with a fresh kernel:

Quit


And then:

a = 600;
reg1 = Rectangle[{0, 0}, {6 a, a}];
hole1 = Disk[{3 a, 0.75 a}, a/8];
hole2 = Disk[{3 a, 0.25 a}, a/8];
hole3 = Disk[{2.5 a, 0.5 a}, a/8];
resreg = RegionDifference[reg1, RegionUnion[hole1, hole2, hole3]];
NDSolveFEMToElementMesh[resreg]


This gives a message:

So now we know that the RegionBounds computation either timed out or was not able to compute the bounds. However, when you call

RegionBounds[resreg]


You will get the bounds. Once the bounds are computed they are cached. That's why we needed the Quit above to make sure that we start from the same state.

As a workaround you can compute the RegionBounds before calling Region, or, if you happen to already have a FEM mesh just use:

 Region[MeshRegion[myFEMMesh]]


or

MeshRegion[myFEMMesh]


If Region cannot do it alone, you try BoundaryDiscretizeRegion with the option MaxCellMeasure to control the resolution of the discretization. In this case, it helps.

length = 0.5;
R = BoundaryDiscretizeRegion[resreg, MaxCellMeasure -> {1 -> length}]


You can obtain a mesh with DiscretizeRegion:

S = DiscretizeRegion[R]


There is a workaround:

a=600;q = ImplicitRegion[x >= 0 && x <= 6*a && y >= 0 &&
y <= a && (x - 3 a)^2 + (y - 3/4*a)^2 >=
a^2/64 && (x - 3 a)^2 + (y - 1/4*a)^2 >=
a^2/64 && (x - 5/2*a)^2 + (y - 1/2*a)^2 >= a^2/64, {x, y}];Region[q]