src/myc.h

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
139
140
141
142
143
144
145
146
147
148
149
150
/*-----------------------------------------------------* 
 *MYC - Mixed Youngs and Central Scheme                *
 *-----------------------------------------------------*/
/* 

Known problems: the index [0,0,0], i.e. the central cell
in the block, never occurs: neither in the central scheme
nor in Youngs' method. Therefore an isolated droplet will have
a normal with all components to zero. I took care of the
division-by-zero issue, but not of this one.

Ruben

*/

coord mycs (Point point, scalar c)
{ 
  double m1,m2,m[4][3],t0,t1,t2;
  int cn;

  /* write the plane as: sgn(mx) X =  my Y +  mz Z + alpha 
                             m00 X = m01 Y + m02 Z + alpha */
  m1 = c[-1,0,-1] + c[-1,0,1] + c[-1,-1,0] + c[-1,1,0] + 
       c[-1,0,0];
  m2 = c[1,0,-1] + c[1,0,1] + c[1,-1,0] + c[1,1,0] + 
       c[1,0,0];
  m[0][0] = m1 > m2 ? 1. : -1.;

  m1 = c[-1,-1,0]+ c[1,-1,0]+ c[0,-1,0];
  m2 = c[-1,1,0]+ c[1,1,0]+ c[0,1,0];
  m[0][1] = 0.5*(m1-m2);

  m1 = c[-1,0,-1]+ c[1,0,-1]+ c[0,0,-1];
  m2 = c[-1,0,1]+ c[1,0,1]+ c[0,0,1];
  m[0][2] = 0.5*(m1-m2);

  /* write the plane as: sgn(my) Y =  mx X +  mz Z + alpha 
                             m11 Y = m10 X + m12 Z + alpha */
  m1 = c[-1,-1,0] + c[-1,1,0] + c[-1,0,0];
  m2 = c[1,-1,0] + c[1,1,0] + c[1,0,0];
  m[1][0] = 0.5*(m1-m2);

  m1 = c[0,-1,-1] + c[0,-1,1] + c[1,-1,0] + c[-1,-1,0] +
       c[0,-1,0];
  m2 = c[0,1,-1] + c[0,1,1] + c[1,1,0] + c[-1,1,0] +
       c[0,1,0];
  m[1][1] = m1 > m2 ? 1. : -1.;

  m1 = c[0,-1,-1]+ c[0,0,-1]+ c[0,1,-1];
  m2 = c[0,-1,1]+ c[0,0,1]+ c[0,1,1];
  m[1][2] = 0.5*(m1-m2);

  /* write the plane as: sgn(mz) Z =  mx X +  my Y + alpha 
                             m22 Z = m20 X + m21 Y + alpha */

  m1 = c[-1,0,-1]+ c[-1,0,1]+ c[-1,0,0];
  m2 = c[1,0,-1]+ c[1,0,1]+ c[1,0,0];
  m[2][0] = 0.5*(m1-m2);

  m1 = c[0,-1,-1]+ c[0,-1,1]+ c[0,-1,0];
  m2 = c[0,1,-1]+ c[0,1,1]+ c[0,1,0];
  m[2][1] = 0.5*(m1-m2);

  m1 = c[-1,0,-1] + c[1,0,-1] + c[0,-1,-1] + c[0,1,-1] +
       c[0,0,-1];
  m2 = c[-1,0,1] + c[1,0,1] + c[0,-1,1] + c[0,1,1] +
       c[0,0,1];
  m[2][2] = m1 > m2 ? 1. : -1.;

  /* normalize each set (mx,my,mz): |mx|+|my|+|mz| = 1 */
  t0 = fabs(m[0][0]) + fabs(m[0][1]) + fabs(m[0][2]);
  m[0][0] /= t0;
  m[0][1] /= t0;
  m[0][2] /= t0;

  t0 = fabs(m[1][0]) + fabs(m[1][1]) + fabs(m[1][2]);
  m[1][0] /= t0;
  m[1][1] /= t0;
  m[1][2] /= t0;

  t0 = fabs(m[2][0]) + fabs(m[2][1]) + fabs(m[2][2]);
  m[2][0] /= t0;
  m[2][1] /= t0;
  m[2][2] /= t0;

  /* choose among the three central scheme */ 
  t0 = fabs(m[0][0]);
  t1 = fabs(m[1][1]);
  t2 = fabs(m[2][2]);

  cn = 0;
  if (t1 > t0) {
    t0 = t1;
    cn = 1;
  }
  if (t2 > t0)
    cn = 2;

  /* Youngs-CIAM scheme */  
  m1 = c[-1,-1,-1] + c[-1,1,-1] + c[-1,-1,1] + c[-1,1,1] +
       2.*(c[-1,-1,0] + c[-1,1,0] + c[-1,0,-1] + c[-1,0,1]) +
       4.*c[-1,0,0];
  m2 = c[1,-1,-1] + c[1,1,-1] + c[1,-1,1] + c[1,1,1] +
       2.*(c[1,-1,0] + c[1,1,0] + c[1,0,-1] + c[1,0,1]) +
       4.*c[1,0,0];
  m[3][0] = m1 - m2;

  m1 = c[-1,-1,-1] + c[-1,-1,1] + c[1,-1,-1] + c[1,-1,1] +
       2.*( c[-1,-1,0] + c[1,-1,0] + c[0,-1,-1] + c[0,-1,1]) +
       4.*c[0,-1,0];
  m2 = c[-1,1,-1] + c[-1,1,1] + c[1,1,-1] + c[1,1,1] +
       2.*(c[-1,1,0] + c[1,1,0] + c[0,1,-1] + c[0,1,1]) +
       4.*c[0,1,0];
  m[3][1] = m1 - m2;

  m1 = c[-1,-1,-1] + c[-1,1,-1] + c[1,-1,-1] + c[1,1,-1] +
       2.*(c[-1,0,-1] + c[1,0,-1] + c[0,-1,-1] + c[0,1,-1]) +
       4.*c[0,0,-1];
  m2 = c[-1,-1,1] + c[-1,1,1] + c[1,-1,1] + c[1,1,1] +
       2.*(c[-1,0,1] + c[1,0,1] + c[0,-1,1] + c[0,1,1]) +
       4.*c[0,0,1];
  m[3][2] = m1 - m2;

  /* normalize the set (mx,my,mz): |mx|+|my|+|mz| = 1 */
  t0 = fabs(m[3][0]) + fabs(m[3][1]) + fabs(m[3][2]);  
  if (t0 < 1e-30) {
    coord mxyz = {1., 0., 0.};
    return mxyz;
  }
    
  m[3][0] /= t0;
  m[3][1] /= t0;
  m[3][2] /= t0;

  /* choose between the previous choice and Youngs-CIAM */
  t0 = fabs (m[3][0]);
  t1 = fabs (m[3][1]);
  t2 = fabs (m[3][2]);
  if (t1 > t0)
    t0 = t1;
  if (t2 > t0)
    t0 = t2;

  if (fabs(m[cn][cn]) > t0)
    cn = 3;

  /* components of the normal vector */
  coord mxyz = {m[cn][0], m[cn][1], m[cn][2]};
  return mxyz; 
}