# RegionIntersection returns unevaluated for 3D regions

I'm trying to find the volume an ellipsoid and a prism share. I tried using RegionIntersection, but it returns unevaluated. Input:

sample3D=Ellipsoid[{0, 0, 0}, {572, 96, 572}];
beam3D=Prism[{{2864.79, 0, -2810}, {0, 10.0001, -2810}, {-2864.79,
0, -2810}, {2864.79, 0, 2810}, {0, 10.0001, 2810}, {-2864.79, 0,
2810}}];
Volume@RegionIntersection[beam3D, sample3D]


Returns:

Volume[BooleanRegion[#1 && #2 &, {Prism[{{1145.93, 0, -2810}, {0,
10.0004, -2810}, {-1145.93, 0, -2810}, {1145.93, 0, 2810}, {0,
10.0004, 2810}, {-1145.93, 0, 2810}}],
Ellipsoid[{0, 0, 0}, {572, 96, 572}]}]]


Both the ellipsoid and prism are of the same dimension and are recognized by RegionQ as regions.

A two-dimensional example of a triangle and an ellipse works fine. Input:

sample2D=Disk[{0, 0}, {572, 10}];
beam2D=Triangle[{{2864.79, 0}, {0, 10.0001}, {-2864.79, 0}}];
Area@RegionIntersection[sample2D, beam2D]


Returns:

8927.06


I already read somewhere, that Mathematica versions below 11.3 had problems with Region functions, but I'm using version 12.0.

Edit:

The original problem can be worked around by applying BoundaryDiscretizeRegion to both regions. Input:

sample3D$$edit1 = BoundaryDiscretizeRegion@Ellipsoid[{0, 0, 0}, {572, 96, 572}]; beam3D$$edit1 = BoundaryDiscretizeRegion@
Prism[{{2864.79, 0, -2810}, {0, 10.0001, -2810}, {-2864.79,
0, -2810}, {2864.79, 0, 2810}, {0, 10.0001, 2810}, {-2864.79, 0,
2810}}];
Volume@RegionIntersection[beam3D$$edit1, sample3D$$edit1]


Returns:

9.37725*10^6


However, this only works for some points for Prism and not for others. Other points generate a bunch of error messages. Input:

sample3D$$edit2 = BoundaryDiscretizeRegion@ Ellipsoid[{0, 0, 0}, {572, 96, 572}]; beam3D$$edit2 = BoundaryDiscretizeRegion@
Prism[{{21824., 11.9101, -2810}, {6823.95, 38.09, -2810}, {-15000,
0, -2810}, {21824., 11.9101, 2810}, {6823.95, 38.09, 2810}, {-15000,
0, 2810}}];
Volume@RegionIntersection[beam3D$$edit2, sample3D$$edit2]


Returns error messages:

BoundaryDiscretizeRegion::regpnd: A non-degenerate region is expected at position 1 of BoundaryDiscretizeRegion[Prism[{{21824.,11.9101,-2810},{6823.95,38.09,-2810},{-15000,0,-2810},{21824.,11.9101,2810},{6823.95,38.09,2810},{-15000,0,2810}}]].
RegionIntersection::reg: BoundaryDiscretizeRegion[Prism[{{21824.,11.9101,-2810},{6823.95,38.09,-2810},{-15000,0,-2810},{21824.,11.9101,2810},{6823.95,38.09,2810},{-15000,0,2810}}]] is not a correctly specified region.
Volume::reg: RegionIntersection[BoundaryDiscretizeRegion[Prism[{{21824.,11.9101,-2810},{6823.95,38.09,-2810},{-15000,0,-2810},{21824.,11.9101,2810},{6823.95,38.09,2810},{-15000,0,2810}}]],] is not a correctly specified region.


Inspired by this answer I added Rationalize to the Prisms points and this works for most of my applications, albeit rather slowly. Input:

sample3D$$edit3 := BoundaryDiscretizeRegion@Ellipsoid[{0, 0, 0}, {572, 96, 572}]; beam3D$$edit3 :=
BoundaryDiscretizeRegion@
Prism[Rationalize[{{150739.44950992631,
23.691003428484592, -2810}, {135739.44950992634,
26.308997333058976, -2810}, {-15000,
0, -2810}, {150739.44950992631, 23.691003428484592,
2810}, {135739.44950992634, 26.308997333058976,
2810}, {-15000, 0, 2810}}, 0]];
Volume@RegionIntersection[beam3D$$edit3, sample3D$$edit3] // AbsoluteTiming


Returns:

{11.139, 486462.}


Eventhough this works, it is too slow for my ultimate goal and still returns error messages for some points:

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

• In Version 13.0 I get a result with: Volume[MeshRegion[NDSolveFEMToElementMesh[RegionIntersection[beam3D, sample3D]]]] Jan 24, 2022 at 8:04
• Strange, I have also version 13.0 on Windows 10 and get the same output like rowsi Jan 24, 2022 at 11:06
• DiscretizeRegion[beam3D] also hangs on 13.0 which isn't great even if you just wanted an approximate result. Wolfram needs better QA for basic cases like this :( Jan 24, 2022 at 11:09
• @flinty, don't misunderstand me this should work and I have reported that internally, but it's certainly not a basic case. The number of characters needed to show an issue does not say anything about the complexity to solve the issue; especially in the wolfram language. Jan 24, 2022 at 11:43
• @DanielHuber, does the version I posted does not work for you in V13? Jan 24, 2022 at 11:44

I managed to fix my problem with a workaround. It consists of two changes to my original problem:

First I split the Prism from my original problem into two prisms whose intersection formed the original prism. Second, I switched from BoundaryDiscretizeRegion to BoundaryDiscretizeGraphics, which works a lot faster. I also added MaxCellMeasure->0.05 to my Ellipsoid. Both higher and lower values for MaxCellMeasure give the error message:

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


But 0.05 seems to be a "sweetspot" and all calculations for my application work with this value. The now working piece looks like this:

incPrism :=
BoundaryDiscretizeGraphics[
Prism[Rationalize[Join[incPrismPtsNegative, incPrismPtsPositive],
0]]];
outPrism :=
BoundaryDiscretizeGraphics[
Prism[Rationalize[Join[outPrismPtsNegative, outPrismPtsPositive],
0]]];
sample3D :=
BoundaryDiscretizeGraphics[
Ellipsoid[{0, 0, 0}, {SampleWidth/2, SampleHeight, SampleWidth/2}],
MaxCellMeasure -> 0.05]
ParallelTable[
Volume@RegionIntersection[outPrism, sample3D], {alpha,
0.01 \[Degree], 0.2 \[Degree], 0.01 \[Degree]}]
ParallelTable[
Volume@RegionIntersection[incPrism, sample3D], {alpha,
0.01 \[Degree], 0.2 \[Degree], 0.01 \[Degree]}]


With SampleWidth and SampleHeight being non-negative numbers for the ellipsoids radii, incPrismPtsNegative, incPrismPtsPositive, outPrismPtsNegative and outPrismPtsPositive being the corresponding points for incPrism and outPrism, which are dependant on the angle alpha, which is used in ParallelTable. Returns:

{2.63661*10^6, 5.26535*10^6, 7.88319*10^6, 1.04872*10^7, 1.30742*10^7,
1.56414*10^7, 1.81854*10^7, 2.07027*10^7, 2.31899*10^7,
2.56421*10^7, 2.80515*10^7, 3.04142*10^7, 3.27272*10^7, 3.49878*10^7,
3.71931*10^7, 3.93402*10^7, 4.14256*10^7, 4.34455*10^7,
4.53931*10^7, 4.55803*10^7}
{4.55883*10^7, 4.55883*10^7, 4.55882*10^7, 4.5588*10^7, 4.55878*10^7,
4.55876*10^7, 4.55873*10^7, 4.5587*10^7, 4.55867*10^7, 4.55863*10^7,
4.55859*10^7, 4.5585*10^7, 4.55849*10^7, 4.55844*10^7, 4.55838*10^7,
4.55832*10^7, 4.55825*10^7, 4.55818*10^7, 4.55811*10^7, 4.55803*10^7}


With further tinkering I found that MaxCellMeasure->0.05 is not the be-all and end-all to my problem. However, playing around with the value for MaxCellMeasure gives satisfactory results most of the time. For my application I still count this as a success.