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
| #define NOT_ZERO 1e-30
/*-----------------------------------------------------*
*MYC - Mixed Youngs and Central Scheme (2D) *
*-----------------------------------------------------*/
coord mycs (Point point, scalar c)
{
int ix;
double c_t,c_b,c_r,c_l;
double mx0,my0,mx1,my1,mm1,mm2;
/* top, bottom, right and left sums of c values */
c_t = c[-1,1] + c[0,1] + c[1,1];
c_b = c[-1,-1] + c[0,-1] + c[1,-1];
c_r = c[1,-1] + c[1,0] + c[1,1];
c_l = c[-1,-1] + c[-1,0] + c[-1,1];
/* consider two lines: sgn(my) Y = mx0 X + alpha,
and: sgn(mx) X = my0 Y + alpha */
mx0 = 0.5*(c_l - c_r);
my0 = 0.5*(c_b - c_t);
/* minimum coefficient between mx0 and my0 wins */
if (fabs(mx0) <= fabs(my0)) {
my0 = my0 > 0. ? 1. : -1.;
ix = 1;
}
else {
mx0 = mx0 > 0. ? 1. : -1.;
ix = 0;
}
/* Youngs' normal to the interface */
mm1 = c[-1,-1] + 2.0*c[-1,0] + c[-1,1];
mm2 = c[1,-1] + 2.0*c[1,0] + c[1,1];
mx1 = mm1 - mm2 + NOT_ZERO;
mm1 = c[-1,-1] + 2.0*c[0,-1] + c[1,-1];
mm2 = c[-1,1] + 2.0*c[0,1] + c[1,1];
my1 = mm1 - mm2 + NOT_ZERO;
/* choose between the best central and Youngs' scheme */
if (ix) {
mm1 = fabs(my1);
mm1 = fabs(mx1)/mm1;
if (mm1 > fabs(mx0)) {
mx0 = mx1;
my0 = my1;
}
}
else {
mm1 = fabs(mx1);
mm1 = fabs(my1)/mm1;
if (mm1 > fabs(my0)) {
mx0 = mx1;
my0 = my1;
}
}
/* normalize the set (mx0,my0): |mx0|+|my0|=1 and
write the two components of the normal vector */
mm1 = fabs(mx0) + fabs(my0);
coord n = {mx0/mm1, my0/mm1};
return n;
}
|
