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