diff --git a/alberta/src/3d/level_3d.c b/alberta/src/3d/level_3d.c index e85b1ba412f0f8edee516283dfda9ff38fb390e4..afe37becaebe8c80ee67e81055d74a40703b3707 100644 --- a/alberta/src/3d/level_3d.c +++ b/alberta/src/3d/level_3d.c @@ -46,7 +46,8 @@ REAL level_element_det_3d(const REAL_D coord[N_VERTICES_2D]) { - REAL_D e1, e2, normal; + REAL_D e1, e2; + REAL det; int i; for (i = 0; i < DIM_OF_WORLD; i++) @@ -55,11 +56,19 @@ REAL level_element_det_3d(const REAL_D coord[N_VERTICES_2D]) e2[i] = coord[2][i] - coord[0][i]; } - normal[0] = e1[1]*e2[2] - e1[2]*e2[1]; - normal[1] = e1[2]*e2[0] - e1[0]*e2[2]; - normal[2] = e1[0]*e2[1] - e1[1]*e2[0]; - - return(sqrt(SCP_DOW(normal,normal))); +#if DIM_OF_WORLD == 2 + det = WEDGE_DOW(e1, e2); + det = ABS(det); +#elif DIM_OF_WORLD == 3 + REAL_D normal; + WEDGE_DOW(e1, e2, normal); + det = NORM_DOW(normal); +#else + det = NRM2_DOW(e1) * NRM2_DOW(e2) - SQR(SCP_DOW(e1, e2)); + det = sqrt(det); +#endif + + return det; } void level_coord_to_world_3d(const REAL_D coord[N_VERTICES_2D],