src/myc2d.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
  
#define NOT_ZERO 1.e-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;
  }