「SDOI2016」平凡的骰子 - 计算几何+三维重心+二面角 | Bill Yang's Blog

「SDOI2016」平凡的骰子 - 计算几何+三维重心+二面角

题目大意

    这是一枚平凡的骰子。它是一个均质凸多面体,表面有$n$个端点,有$f$个面,每一面是一个凸多边形,且任意两面不共面。将这枚骰子抛向空中,骰子落地的时候不会发生二次弹跳(这是一种非常理想的情况)。你希望知道最终每一面着地的概率。每一面着地的概率可以用如下的方法计算:我们假设$O$为骰子的重心,并以$O$为球心,做半径为$1$的单位球面(记为$S$)。我们知道$S$的表面积即单位球的表面积,为$4\pi$,这里$\pi$为圆周率。对于骰子的某一面$C$来说,球面$S$上存在一块区域$T$满足:当下落时若骰子所受重力方向与$S$的交点落在$T$中,则$C$就是最终着地的一面。那么$C$着地的概率为区域$T$的面积除以$4\pi$。为了能更好地辅助计算球面上一块区域的面积,我们给出单位球面$S$上三角形的面积计算公式。考虑单位球面$S$上的三个两两相交的大圆,交点依次为$A$,$B$和$C$。则曲面三角形$ABC$的面积为$\text{Area}(ABC)=\alpha+\beta+\gamma-\pi$,其中$\alpha$,$\beta$和$\gamma$分别对应了三个二面角的大小。如下图所示。

    我们保证:每一面着地的时候,重心的垂心都恰好在这一面内。也就是说不会出现摆不稳的情况。


题目分析

这是一道结论题。
需要知道以下几个东西才能做。

  1. 如何求凸多面体重心?
    答:将凸多面体剖分成若干个四面体,每个四面体计算体积,所有剖分出的向量带权(体积)的复合平均即为重心。
  2. 如何四面体剖分?
    答:对于每个面进行三角形剖分,与剖分点连线即形成四面体。
  3. 如何求剖分出的四面体体积?
    答:剖分出的三个向量作混合积除以$6$既是体积。
    混合积指任意两个做叉积后于另外一个点积。
    叉积得到底面面积,点积得到与高的乘积,除以$6$即为体积。
  4. 如何求二面角?
    答:通过法向量的夹角计算,法向量可以通过叉积得到。
  5. 如何求曲面三角形面积?
    答:$\text{Area}(ABC)=\alpha+\beta+\gamma-\pi$,题目已描述。
  6. 如何求曲面多边形面积?
    答:可以剖分成曲面三角形,也可以对曲面多边形的内角求和减去$(\text{点数}-2)\pi$(根据高斯波捏公式可推导)。

然后就可以做了。
然而这之前以上六个我都不会。


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;
inline const int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(!isdigit(x)) {
if(x=='-')bj=-1;
x=getchar();
}
while(isdigit(x)) {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}
const int maxn=105;
const double pi=acos(-1);
struct Point {
double x,y,z;
Point(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
Point operator + (const Point& b) {
return Point(x+b.x,y+b.y,z+b.z);
}
void operator += (const Point& b) {
*this=*this+b;
}
Point operator - (const Point& b) {
return Point(x-b.x,y-b.y,z-b.z);
}
void operator -= (const Point& b) {
*this=*this-b;
}
Point operator * (double c) {
return Point(x*c,y*c,z*c);
}
void operator *= (double c) {
*this=*this*c;
}
Point operator / (double c) {
return Point(x/c,y/c,z/c);
}
void operator /= (double c) {
*this=*this/c;
}
double length() {
return sqrt(x*x+y*y+z*z);
}
} p[maxn];
typedef Point Vector;
double Dot(Vector a,Vector b) {
return a.x*b.x+a.y*b.y+a.z*b.z;
}
Vector Cross(Vector a,Vector b) {
return Vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
double Triple(Vector a,Vector b,Vector c) {
return Dot(a,Cross(b,c));
}
double Cal(Vector a,Vector b,Vector c) {
Vector u=Cross(a,b),v=Cross(a,c);
return acos(Dot(u,v)/u.length()/v.length());
}
int n,f,F[maxn][maxn];
int main() {
n=Get_Int();
f=Get_Int();
for(int i=1; i<=n; i++) {
double x,y,z;
scanf("%lf%lf%lf",&x,&y,&z);
p[i]=Point(x,y,z);
}
for(int i=1; i<=f; i++) {
F[i][0]=Get_Int();
for(int j=1; j<=F[i][0]; j++)F[i][j]=Get_Int();
}
double vol=0;
Vector ctr;
for(int i=1; i<=f; i++)
for(int j=2; j<F[i][0]; j++) {
Vector ci=p[F[i][1]]+p[F[i][j]]+p[F[i][j+1]];
double vi=Triple(p[F[i][1]],p[F[i][j]],p[F[i][j+1]]);
vol+=vi;
ctr+=ci*vi;
}
ctr/=vol*4;
for(int i=1; i<=n; i++)p[i]-=ctr;
for(int i=1; i<=f; i++) {
double ans=0;
for(int j=2; j<F[i][0]; j++)ans+=Cal(p[F[i][j]],p[F[i][j-1]],p[F[i][j+1]]);
ans+=Cal(p[F[i][F[i][0]]],p[F[i][F[i][0]-1]],p[F[i][1]])+Cal(p[F[i][1]],p[F[i][F[i][0]]],p[F[i][2]]);
printf("%0.7lf\n",(ans-(F[i][0]-2)*pi)/(4*pi));
}
return 0;
}
本篇文章有用吗?有用还不快点击!
0%