This commit is contained in:
Stig Rune Jensen
2022-10-24 08:21:35 +02:00
parent 49e399b20d
commit 3e74452aa0
2 changed files with 19 additions and 11 deletions

View File

@@ -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) =
@@ -312,19 +313,25 @@ let private interpolateST (grid: ExtendedGrid) (field: Field4D) (s: float) (t: f
let w0 = interpolateTime (float w00) (float w10)
let w1 = interpolateTime (float w01) (float w11)
let w = 0.0//interpolateSigma w0 w1
let w = interpolateSigma w0 w1
let centroids = grid.getCentroids ()
let neighbors = getAdjacentNeighbors grid eIdx
let c0 = centroids[eIdx]
let cn = neighbors |> Array.map (Array.get centroids)
let h0 = 1.0 / (field.Bath[eIdx] |> float)
let hn = neighbors |> Array.map (fun e -> 1.0 / (field.Bath[e] |> float))
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*u}, {hy*v}, {h0*w}"
//printfn $"dh = {hx}, {hy}, {h0}"
u, v, 0.0//-u*hx - v*hy - w*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
@@ -502,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)

View File

@@ -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