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
| coord interface_centroid(Point point,scalar f){
coord n = mycs (point,f);
double alpha = plane_alpha (f[], n);
coord centroid = {0.,0.,0.};
int nrcross=0.;
if (n.x!=0){// it could cross an x-edge
for (double zz=-0.5;zz<=0.5;zz+=1.){
for (double yy=-0.5;yy<=0.5;yy+=1.){
double tempx = (alpha-n.y*yy-n.z*zz)/n.x;
if (fabs(tempx)<=0.5){//Crosses the 'edge'
nrcross++;
centroid.x+=tempx;
centroid.y+=yy;
centroid.z+=zz;
}
}
}
}
if (n.y!=0){// it could cross an y-edge
for (double zz=-0.5;zz<=0.5;zz+=1.){
for (double xx=-0.5;xx<=0.5;xx+=1.){
double tempy = (alpha-n.x*xx-n.z*zz)/n.y;
if (fabs(tempy)<=0.5){//Crosses the 'edge'
nrcross++;
centroid.x+=xx;
centroid.y+=tempy;
centroid.z+=zz;
}
}
}
}
if (n.z!=0){// it could cross an z-edge
for (double xx=-0.5;xx<=0.5;xx+=1.){
for (double yy=-0.5;yy<=0.5;yy+=1.){
double tempz = (alpha-n.y*yy-n.x*xx)/n.z;
if (fabs(tempz)<=0.5){//Crosses the edge
nrcross++;
centroid.x+=xx;
centroid.y+=yy;
centroid.z+=tempz;
}
}
}
}
foreach_dimension() //normalize
centroid.x/=(double)nrcross;
foreach_dimension() //convert to grid units
centroid.x=x+(centroid.x*Delta);
return centroid;
}
double cell_interface_area(Point point,scalar c){ //Copied from fractions.h
if (c[] > 0 && c[] < 1.) {
coord n = mycs (point, c), p;
double alpha = plane_alpha (c[], n);
double area = pow(Delta, dimension - 1)*plane_area_center (n, alpha, &p);
return area;
}else{
return 0.;
}
}
double interface_average(scalar g,scalar c){
double inte = 0.;
foreach(reduction(+:inte)){
if (c[]>0. && c[]<1.){//Interface
coord location = interface_centroid(point,c);
double A = cell_interface_area(point,c);
inte+=interpolate(g,location.x,location.y,location.z)*A;
}
}
double A=interface_area(c);
if (A!=0){
return inte/interface_area(c);
}else{
return 0;
}
}
|