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.
Both computers have Windows 10 with Mathematica v 11.3. How I can make Region plot right image?
Answer
This is a bug and you should report it to support@wolfram.com. 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]];
NDSolve`FEM`ToElementMesh[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]
Comments
Post a Comment