3

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.

Thanks a bunch in advance!

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.
rowsi
  • 395
  • 1
  • 10
  • 4
    In Version 13.0 I get a result with: Volume[MeshRegion[NDSolveFEMToElementMesh[RegionIntersection[beam3D, sample3D]]]] – user21 Jan 24 '22 at 08:04
  • Strange, I have also version 13.0 on Windows 10 and get the same output like rowsi – Daniel Huber Jan 24 '22 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 :( – flinty Jan 24 '22 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. – user21 Jan 24 '22 at 11:43
  • @DanielHuber, does the version I posted does not work for you in V13? – user21 Jan 24 '22 at 11:44
  • Thanks a lot for your comments so far. I have been able to get a result by applying BoundaryDiscretizeRegion to both regions. Eventhough it works for the points mentioned in my question, it still doesn't evaluate for others. I'll update my question accordingly. – rowsi Jan 24 '22 at 13:32
  • @user21 Sadly, I am not able to upgrade to 13.0 as there are some issues with access to the data on our network drive at my institution with this version. Local copies of my notebooks would work, but administrator credentials aren't available to me. I contacted our IT guys to resolve the issue. – rowsi Jan 24 '22 at 14:38
  • Hi rowsi, I recommend when you update your post to include different versions of such things as beam3D and related terms, you separate their definitions out by changing the names. For example, beam3D$edit1 and so-on. That way, future users who want to help you can better differentiate between what is working and what isn’t, and then reference these things in their answers. Can you, please, make this update? – CA Trevillian Jan 24 '22 at 15:38
  • @CATrevillian edited as per your suggestion. Thanks a bunch, gonna keep that in mind in the future. – rowsi Jan 24 '22 at 15:57
  • @user21 As I wrote, I do not get your result but the result rowsi gets. – Daniel Huber Jan 24 '22 at 16:26
  • With version 13.0.0 Volume[MeshRegion[NDSolve`FEM`ToElementMesh[RegionIntersection[beam3D, sample3D]]]] returns 9.33744*10^6 for the original input. But for sample3D$edit2 and beam3D$edit2 it doesn't work (both with and without BoundaryDiscretizeRegion). – Alexey Popkov Jan 25 '22 at 09:07

1 Answers1

1

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}

Thanks a lot for your comments everyone!

Edit

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.

rowsi
  • 395
  • 1
  • 10