From 2c477f13ad1c437a526190773159b1f2d9306a09 Mon Sep 17 00:00:00 2001
From: "guo.2154" <guo.2154@cse-anc240738l.mynetworksettings.com>
Date: Sun, 29 Oct 2023 12:13:39 -0400
Subject: [PATCH] WIP: normal velocity interp

---
 include/ftk/mesh/mpas_mesh.hh             | 20 +++++++++++++++-----
 include/ftk/numeric/vector_dot_product.hh |  7 +++++++
 include/ftk/numeric/vector_norm.hh        |  8 ++++++--
 3 files changed, 28 insertions(+), 7 deletions(-)

diff --git a/include/ftk/mesh/mpas_mesh.hh b/include/ftk/mesh/mpas_mesh.hh
index 7a71d359..2f7ce415 100644
--- a/include/ftk/mesh/mpas_mesh.hh
+++ b/include/ftk/mesh/mpas_mesh.hh
@@ -266,10 +266,11 @@ void mpas_mesh<I, F>::initialize_cell_tangent_plane()
   if (edgeNormalVectors.empty())
     initialize_edge_normals();
 
-  cellTangentPlane.reshape(2, n_cells());
+  cellTangentPlane.reshape(3, 2, n_cells());
   F p[2], // tangent plane
     rhat[3], // unit vector for a cell,
     xhat[3],
+    yhat[3],
     n[3]; // edge normal of the first edge
 
   for (auto ci = 0; ci < n_cells(); ci ++) {
@@ -284,9 +285,15 @@ void mpas_mesh<I, F>::initialize_cell_tangent_plane()
 
     for (int i = 0; i < 3; i ++)
       xhat[i] = n[i] - normal_dot_r * rhat[i];
-   
-    cellTangentPlane(0, ci) = xhat[0];
-    cellTangentPlane(1, ci) = xhat[1];
+    vector_normalization2_3(xhat);
+
+    cross_product(rhat, xhat, yhat);
+    vector_normalization2_3(yhat); // not necesary but just make sure
+ 
+    for (int i = 0; i < 3; i ++) {
+      cellTangentPlane(i, 0, ci) = xhat[i];
+      cellTangentPlane(i, 1, ci) = yhat[i];
+    }
   }
 }
 
@@ -317,7 +324,10 @@ void mpas_mesh<I, F>::initialize_coeffs_reconstruct()
       alpha += vector_dist_2norm_3(x, Xe[i]);
     }
 
-    const F t[2] = {cellTangentPlane(0, ci), cellTangentPlane(1, ci)};
+    const F t[2][3] = {
+      {cellTangentPlane(0, 0, ci), cellTangentPlane(1, 0, ci), cellTangentPlane(2, 0, ci)},
+      {cellTangentPlane(0, 1, ci), cellTangentPlane(1, 1, ci), cellTangentPlane(2, 1, ci)}
+    };
 
     F coeffs[max_edges];
     rbf3d_plane_vec_const_dir(
diff --git a/include/ftk/numeric/vector_dot_product.hh b/include/ftk/numeric/vector_dot_product.hh
index 29343ae6..fb8a224e 100644
--- a/include/ftk/numeric/vector_dot_product.hh
+++ b/include/ftk/numeric/vector_dot_product.hh
@@ -3,6 +3,13 @@
 
 namespace ftk {
 
+template <typename T>
+__device__ __host__
+inline T vector_dot_product2(const T A[2], const T B[2])
+{
+  return A[0]*B[0] + A[1]*B[1];
+}
+
 template <typename T>
 __device__ __host__
 inline T vector_dot_product3(const T A[3], const T B[3])
diff --git a/include/ftk/numeric/vector_norm.hh b/include/ftk/numeric/vector_norm.hh
index 03a20f56..616edd13 100644
--- a/include/ftk/numeric/vector_norm.hh
+++ b/include/ftk/numeric/vector_norm.hh
@@ -68,8 +68,12 @@ template <int n, typename T>
 __device__ __host__
 inline T vector_dist_2norm2(const T v[], const T w[])
 {
-  T u[3] = {v[0] - w[0], v[1] - w[1], v[2] - w[2]};
-  return vector_2norm2<n, T>(u);
+  T d2(0);
+  for (int i=0; i<n; i++) {
+    T d = v[i] - w[i];
+    d2 += d * d;
+  }
+  return d2;
 }
 
 }