tmp2

src/Evaluate4D.fs

@@ -20,6 +20,7 @@ let private toArray3 (a, b, c) = [| a; b; c |]
 
 let private getCell (grid: ExtendedGrid) (idx: ElemIdx) = (grid :> IGrid).getCell idx
 let private getVertex (grid: ExtendedGrid) (idx: NodeIdx) = (grid :> IGrid).getVertex idx
+let private getCellVertices (grid: ExtendedGrid) (idx: NodeIdx) = (grid :> IGrid).getCellVertices idx
 let private getCentroid (grid: ExtendedGrid) (idx: ElemIdx) = (grid.getCentroids ())[idx]
 
 let private containsIndex (idx: NodeIdx) (cell: Cell) =
@@ -277,7 +278,7 @@ let private calcSigmaLayers (field: Field4D) (sigma: float) : int * int =
     | _ -> (slayer - 1, slayer)
 
 // Linear interpolation in sigma and time in given cell
-let private interpolateST (field: Field4D) (s: float) (t: float) (eIdx: ElemIdx) =
+let private interpolateST (grid: ExtendedGrid) (field: Field4D) (s: float) (t: float) (eIdx: ElemIdx) =
     let slayer0, slayer1 = calcSigmaLayers field s
     let s0 = field.Siglays[slayer0] |> float
     let s1 = field.Siglays[slayer1] |> float
@@ -314,8 +315,23 @@ let private interpolateST (field: Field4D) (s: float) (t: float) (eIdx: ElemIdx)
     let w1 = interpolateTime (float w01) (float w11)
     let w = interpolateSigma w0 w1
 
-    let h = field.Bath[eIdx] |> float
-    (u, v, w / h)
+    let centroids = grid.getCentroids ()
+    let c0 = centroids[eIdx]
+    let h0 = (field.BathCenter[eIdx] |> float)
+
+    // let neighbors = getAdjacentNeighbors grid eIdx
+    // let cn = neighbors |> Array.map (Array.get centroids)
+    // let hn = neighbors |> Array.map (fun n -> (field.BathCenter[n] |> float))
+
+    let vertices = (getCell grid eIdx) |> (fun (a,b,c) -> [| a; b; c |])
+    let cn = getCellVertices grid eIdx |> (fun (a,b,c) -> [| a; b; c |])
+    let hn = vertices |> Array.map (fun n -> (field.BathVertex[n] |> float))
+
+    let thrs = 1.0e-6
+    let hx, hy = calcGradientFromNeighbors c0 h0 cn hn thrs
+    //printfn $"dh = {hx}, {hy}, {h0}"
+
+    u, v, (w - s*hx*u - s*hy*v)/h0
 
 // Evaluate 4D field based on constant horizontal extrapolation from centroid field value
 // Linear interpolation in vertical and time
@@ -324,7 +340,7 @@ let evaluateConstantCentroid4D (grid: ExtendedGrid) (field: Field4D) (t: float)
     | None -> None
     | Some e ->
         // Set up linear sigma and time interpolation
-        let readUVW eIdx = interpolateST field s t eIdx
+        let readUVW eIdx = interpolateST grid field s t eIdx
 
         let u = readUVW e |> fst3
         let v = readUVW e |> snd3
@@ -340,7 +356,7 @@ let evaluateLinearCentroid4D (grid: ExtendedGrid) (field: Field4D) (t: float) ((
     | None -> None
     | Some e ->
         // Set up linear sigma and time interpolation
-        let readUVW eIdx = interpolateST field s t eIdx
+        let readUVW eIdx = interpolateST grid field s t eIdx
 
         let centroids = grid.getCentroids ()
         let c0 = centroids[e]
@@ -381,7 +397,7 @@ let evaluateCloughTocherCentroid4D (grid: ExtendedGrid) (field: Field4D) (t: flo
     | None -> None
     | Some e ->
         // Set up linear sigma and time interpolation
-        let readUVW eIdx = interpolateST field s t eIdx
+        let readUVW eIdx = interpolateST grid field s t eIdx
 
         // Position of vertices
         let idx = (getCell grid e) |> toArray3
@@ -443,7 +459,7 @@ let evaluateCloughTocherVertex4D (grid: ExtendedGrid) (field: Field4D) (t: float
     | None -> None
     | Some e ->
         // Set up linear sigma and time interpolation
-        let readUVW eIdx = interpolateST field s t eIdx
+        let readUVW eIdx = interpolateST grid field s t eIdx
         let readU = fst3 << readUVW
         let readV = snd3 << readUVW
         let readW = trd3 << readUVW
@@ -493,4 +509,4 @@ let evaluateCloughTocherVertex4D (grid: ExtendedGrid) (field: Field4D) (t: float
         let u1 = interpolateTriangleCT p_n u_n' du_n' (x, y)
         let v1 = interpolateTriangleCT p_n v_n' dv_n' (x, y)
         let w1 = interpolateTriangleCT p_n w_n' dw_n' (x, y)
-        Some(u1, v1, w1)
+        Some(u1, v1, w1)

src/Grid.fs

@@ -31,7 +31,8 @@ type Field4D =
         Field0: (single * single * single) [] []
         Field1: (single * single * single) [] []
         Siglays: single []
-        Bath: single []
+        BathCenter: single []
+        BathVertex: single []
     }
 
 type BBox =
@@ -374,7 +375,7 @@ module Util =
             let Sy = det Tx T2 E |> (*) 0.5
             let a = det Tx Ty E
             let b = det Tx Ty T2
-            if abs a < 1.0e-12 then failwith "co-linear vertices" 
+            if abs a < 1.0e-12 then failwith "co-linear vertices"
             let S2 = (Sx, Sy) |> square
             let r = (b / a + S2 / (a * a)) |> sqrt
             Sx / a, Sy / a, r
@@ -701,4 +702,4 @@ type ExtendedGrid(grid: IGrid) =
             elementTree <- binarySerializer.UnPickle<IdxTree> pickle |> Some
             true
         else
-            false
+            false