/*-----------------------------------------------------*
*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;
}
