-
Notifications
You must be signed in to change notification settings - Fork 6
/
volume_fractions.cpp
139 lines (126 loc) · 5.43 KB
/
volume_fractions.cpp
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#include "util.h"
#include "volume_fractions.h"
// Assumes phi0<0 and phi1>=0, phi2>=0, or vice versa.
// In particular, phi0 must not equal either of phi1 or phi2.
template<class T>
static T sorted_triangle_fraction(T phi0, T phi1, T phi2)
{
return phi0*phi0/(2*(phi0-phi1)*(phi0-phi2));
}
float area_fraction(float phi0, float phi1, float phi2)
{
if(phi0<0){
if(phi1<0){
if(phi2<0) return 0;
else return 1-sorted_triangle_fraction(phi2, phi0, phi1);
}else if(phi2<0) return 1-sorted_triangle_fraction(phi1, phi2, phi0);
else return sorted_triangle_fraction(phi0, phi1, phi2);
}else if(phi1<0){
if(phi2<0) return 1-sorted_triangle_fraction(phi0, phi1, phi2);
else return sorted_triangle_fraction(phi1, phi2, phi0);
}else if(phi2<0) return sorted_triangle_fraction(phi2, phi0, phi1);
else return 0;
}
double area_fraction(double phi0, double phi1, double phi2)
{
if(phi0<0){
if(phi1<0){
if(phi2<0) return 0;
else return 1-sorted_triangle_fraction(phi2, phi0, phi1);
}else if(phi2<0) return 1-sorted_triangle_fraction(phi1, phi2, phi0);
else return sorted_triangle_fraction(phi0, phi1, phi2);
}else if(phi1<0){
if(phi2<0) return 1-sorted_triangle_fraction(phi0, phi1, phi2);
else return sorted_triangle_fraction(phi1, phi2, phi0);
}else if(phi2<0) return sorted_triangle_fraction(phi2, phi0, phi1);
else return 0;
}
float area_fraction(float phi00, float phi10, float phi01, float phi11)
{
float phimid=(phi00+phi10+phi01+phi11)/4;
return (area_fraction(phi00, phi10, phimid)
+area_fraction(phi10, phi11, phimid)
+area_fraction(phi11, phi01, phimid)
+area_fraction(phi01, phi00, phimid))/4;
}
double area_fraction(double phi00, double phi10, double phi01, double phi11)
{
double phimid=(phi00+phi10+phi01+phi11)/4;
return (area_fraction(phi00, phi10, phimid)
+area_fraction(phi10, phi11, phimid)
+area_fraction(phi11, phi01, phimid)
+area_fraction(phi01, phi00, phimid))/4;
}
//============================================================================
// Assumes phi0<0 and phi1>=0, phi2>=0, and phi3>=0 or vice versa.
// In particular, phi0 must not equal any of phi1, phi2 or phi3.
template<class T>
static T sorted_tet_fraction(T phi0, T phi1, T phi2, T phi3)
{
return phi0*phi0*phi0/((phi0-phi1)*(phi0-phi2)*(phi0-phi3));
}
// Assumes phi0<0, phi1<0, and phi2>=0, and phi3>=0 or vice versa.
// In particular, phi0 and phi1 must not equal any of phi2 and phi3.
template<class T>
static T sorted_prism_fraction(T phi0, T phi1, T phi2, T phi3)
{
T a=phi0/(phi0-phi2),
b=phi0/(phi0-phi3),
c=phi1/(phi1-phi3),
d=phi1/(phi1-phi2);
return a*b*(1-d)+b*(1-c)*d+c*d;
}
float volume_fraction(float phi0, float phi1, float phi2, float phi3)
{
sort(phi0, phi1, phi2, phi3);
if(phi3<=0) return 1;
else if(phi2<=0) return 1-sorted_tet_fraction(phi3, phi2, phi1, phi0);
else if(phi1<=0) return sorted_prism_fraction(phi0, phi1, phi2, phi3);
else if(phi0<=0) return sorted_tet_fraction(phi0, phi1, phi2, phi3);
else return 0;
}
double volume_fraction(double phi0, double phi1, double phi2, double phi3)
{
sort(phi0, phi1, phi2, phi3);
if(phi3<=0) return 1;
else if(phi2<=0) return 1-sorted_tet_fraction(phi3, phi2, phi1, phi0);
else if(phi1<=0) return sorted_prism_fraction(phi0, phi1, phi2, phi3);
else if(phi0<=0) return sorted_tet_fraction(phi0, phi1, phi2, phi3);
else return 0;
}
float volume_fraction(float phi000, float phi100,
float phi010, float phi110,
float phi001, float phi101,
float phi011, float phi111)
{
// This is the average of the two possible decompositions of the cube into
// five tetrahedra.
return ( volume_fraction(phi000, phi001, phi101, phi011)
+volume_fraction(phi000, phi101, phi100, phi110)
+volume_fraction(phi000, phi010, phi011, phi110)
+volume_fraction(phi101, phi011, phi111, phi110)
+2*volume_fraction(phi000, phi011, phi101, phi110)
+volume_fraction(phi100, phi101, phi001, phi111)
+volume_fraction(phi100, phi001, phi000, phi010)
+volume_fraction(phi100, phi110, phi111, phi010)
+volume_fraction(phi001, phi111, phi011, phi010)
+2*volume_fraction(phi100, phi111, phi001, phi010))/12;
}
double volume_fraction(double phi000, double phi100,
double phi010, double phi110,
double phi001, double phi101,
double phi011, double phi111)
{
// This is the average of the two possible decompositions of the cube into
// five tetrahedra.
return ( volume_fraction(phi000, phi001, phi101, phi011)
+volume_fraction(phi000, phi101, phi100, phi110)
+volume_fraction(phi000, phi010, phi011, phi110)
+volume_fraction(phi101, phi011, phi111, phi110)
+2*volume_fraction(phi000, phi011, phi101, phi110)
+volume_fraction(phi100, phi101, phi001, phi111)
+volume_fraction(phi100, phi001, phi000, phi010)
+volume_fraction(phi100, phi110, phi111, phi010)
+volume_fraction(phi001, phi111, phi011, phi010)
+2*volume_fraction(phi100, phi111, phi001, phi010))/12;
}