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