From 8d5f3fca837524b70bddb62607e7bcf9ead23018 Mon Sep 17 00:00:00 2001
From: Claus-Justus Heine <Claus-Justus.Heine@IANS.Uni-Stuttgart.DE>
Date: Fri, 22 Jul 2016 09:58:16 +0000
Subject: [PATCH] Remove by-hand vector arithmetic.

---
 alberta/src/3d/level_3d.c | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/alberta/src/3d/level_3d.c b/alberta/src/3d/level_3d.c
index e85b1ba..afe37be 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],
-- 
GitLab