Cube Tetrahedron Picking

Author

Eric W. Weisstein
March 30, 2004

This notebook downloaded from http://mathworld.wolfram.com/notebooks/ComputationalGeometry/CubeTetrahedronPicking.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/CubeTetrahedronPicking.html.

©2005 Wolfram Research, Inc. except for portions noted otherwise

Code

Options[IntegrateWithAssumptions]=
{
Simplification->Simplify,
Assumptions->True
};

IntegrateWithAssumptions[integrand_,dx__,opts___?OptionQ]:=Module[
{
ineqs=#2<#1<#3&@@@{dx},
conds,
simp=Simplification/.{opts}/.Options[IntegrateWithAssumptions],
assum=Assumptions/.{opts}/.Options[IntegrateWithAssumptions]
},
conds=And@@@NestList[Most,Most[ineqs],Length[ineqs]-1];
Fold[
Simplify[Integrate[#1,#2[[1]],Assumptions->assum&&#2[[2]]],assum&&#2[[2]]]&,
integrand,Transpose[{Reverse[{dx}],conds}]]
]

IntegrateWithAssumptions2[integrand_,dx__,opts___?OptionQ]:=Module[
{allbutlast,
int,
x,x0,x1
},
{x,x0,x1}=First[{dx}];
allbutlast=IntegrateWithAssumptions[integrand,Sequence@@Rest[{dx}],opts,Assumptions->x0<x<x1];
int=Integrate[allbutlast,{dx}[[1,1]]];
Simplify[Limit[int,x->x1,Direction->1]-Limit[int,x->x0,Direction->-1]]
]

Figure

t=Tetrahedron[Table[Random[],{4},{3}]]

t={Polygon[{{0.7889609391893055`,0.32272043783455256`,0.12747961419004444`},{0.6714964275851176`,0.31917719893703433`,0.9101851914643386`},{0.13755919918922033`,0.39724563103322286`,0.48812768116071337`}}],Polygon[{{0.7889609391893055`,0.32272043783455256`,0.12747961419004444`},{0.6714964275851176`,0.31917719893703433`,0.9101851914643386`},{0.637677796472786`,0.70362487416762`,0.3451062331385839`}}],Polygon[{{0.7889609391893055`,0.32272043783455256`,0.12747961419004444`},{0.13755919918922033`,0.39724563103322286`,0.48812768116071337`},{0.637677796472786`,0.70362487416762`,0.3451062331385839`}}],Polygon[{{0.6714964275851176`,0.31917719893703433`,0.9101851914643386`},{0.13755919918922033`,0.39724563103322286`,0.48812768116071337`},{0.637677796472786`,0.70362487416762`,0.3451062331385839`}}]};

Show[Graphics3D[{WireFrame[Cuboid[{0,0,0}]][[1]],t},Boxed->False]]

LiveGraphics3D

WriteLiveForm["cubetet.m",Show[Graphics3D[{WireFrame[Cuboid[{0,0,0}]][[1]],Orient[t]/.Polygon[x_]:>Polygon[c+#&/@x]},Boxed->False]]]

Numerical Integration

Setup  array and do test integrations

a=Array[x,{4,3}];

{#,0,1}&/@Flatten[Array[x,{4,3}]]

MonteCarlo

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->MonteCarlo]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->MonteCarlo,MaxPoints->10^6,PrecisionGoal->10]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->MonteCarlo,MaxPoints->10^7,PrecisionGoal->10]//Timing

v6

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->MonteCarlo,MaxPoints->10^7,PrecisionGoal->10]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->MonteCarlo,MaxPoints->10^7,PrecisionGoal->10]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->MonteCarlo,MaxPoints->10^8,PrecisionGoal->10]//Timing

QuasiMonteCarlo

NIntegrate[Abs[TetrahedronVolume[a]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->2]//Timing

NIntegrate[Abs[TetrahedronVolume[a]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^7,PrecisionGoal->3]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->10]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->10]//Timing

QuasiMonteCarlo rolls over and gives wrong results (bug!) for (10^7)points

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^7,PrecisionGoal->10]//Timing

...this is wrong...

v6

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->10]//Timing

NIntegrate[Evaluate[Abs[TetrahedronVolume[a]]],Evaluate[Sequence@@({#,0,1}&/@Flatten[a])],Method->QuasiMonteCarlo,MaxPoints->10^7,PrecisionGoal->10]//Timing

Brute Force Exact Integration

a=Array[x,{4,3}];

{#,0,1}&/@Flatten[Array[x,{4,3}]]

TetrahedronVolume[a]

Rather wrong

Integrate[TetrahedronVolume[a],Sequence@@({#,0,1}&/@Flatten[a])]

Simplify`PiecewiseIntegrate;
PiecewiseDump`$ShowIntegrals=True;
i=Integrate[Abs[TetrahedronVolume[a]],Sequence@@({#,0,1}&/@Flatten[a])];//Timing

...uses > 1.5 GB...

...ran out of memory after 90 minutes...

Exact Integration Using Geometric Probability

By Zinani (2003)

N[3977/216000-π^2/2160,20]

Intermediate Quantities

Area[l_]:=Plus@@(Det/@Partition[l,2,1,1])/2

Quadrilateral

p={{0,b,1},{a,0,1},{0,c,0},{a c/b,0,0}};

p2={{0,0},{1,0},{c/b,1},{0,1}};

Show[Graphics[Polygon[p2]]/.{c->.2,b->1},AspectRatio->Automatic]

Together[Area[p2]]

a[4]=With[{
V=((a(b^2+b c+c^2))/(6b)),
A=((b+c)/(2b)),
j=a^2 b,
e=((3b^4+12b^3c+14b^2c^2+12b c^3+3c^4)/(36(b+c)^4))
},
4Integrate[(V^2+(1-V)^2)j A^4 e,{b,0,1},{a,b,1},{c,0,b/a}]
]//Timing