major: remove evaluate and interpolation from FvcomKit
This commit is contained in:
Stig Rune Jensen
@@ -1,100 +0,0 @@
|
||||
module Oceanbox.FvcomKit.CloughTocher
|
||||
|
||||
open Grid
|
||||
|
||||
// Cyclic access of array, (A[A.Length] = A[0]), (A[-1] = A[A.Length-1]), etc
|
||||
let private access (A: float []) (idx: int) : float =
|
||||
let l = A.Length
|
||||
let i = idx % l
|
||||
if i >= 0 then A[i] else A[l + i]
|
||||
|
||||
// Clough-Tocher interpolation of 1D field inside float triangular cell
|
||||
// p_0: position to evaluate, should be within the triangle defined by P_n
|
||||
// p_n, f_n, df_n: Positions, values and gradients at cell vertices
|
||||
// Ref: C.L. Lawson "C1-Compatible Interpolation Over a Triangle" (1976)
|
||||
let interpolateTriangleCT (p_n: Pos []) (f_n: float []) (df_n: Field) (p_0: Pos) : float =
|
||||
let idx = [| 0..2 |]
|
||||
|
||||
// Phase 1:
|
||||
let x_n = p_n |> Array.map fst
|
||||
let y_n = p_n |> Array.map snd
|
||||
|
||||
let u =
|
||||
idx
|
||||
|> Array.map (fun i -> (access x_n (i - 1)) - (access x_n (i + 1)))
|
||||
let v =
|
||||
idx
|
||||
|> Array.map (fun i -> (access y_n (i - 1)) - (access y_n (i + 1)))
|
||||
let l2 = Array.map2 (fun a b -> a * a + b * b) u v
|
||||
|
||||
let delta = u[0] * v[1] - v[0] * u[1]
|
||||
if (abs delta) < 1.0e-12 then failwith "Co-linear vertices!"
|
||||
let d_inv = 1.0 / delta
|
||||
|
||||
let x, y = p_0
|
||||
let x_tilde = x - x_n[0]
|
||||
let y_tilde = y - y_n[0]
|
||||
|
||||
let r_1 = d_inv * (u[1] * y_tilde - v[1] * x_tilde)
|
||||
let r_2 = d_inv * (u[2] * y_tilde - v[2] * x_tilde)
|
||||
let r_0 = 1.0 - (r_1 + r_2)
|
||||
let r = [| r_0; r_1; r_2 |]
|
||||
|
||||
let exterior = r |> Array.exists (fun r_i -> r_i < 0.0)
|
||||
if exterior then failwith "Point is outside the triangle!"
|
||||
|
||||
let phi =
|
||||
idx
|
||||
|> Array.map (fun i -> (access r (i + 1)) * (access r (i - 1)))
|
||||
|
||||
// Phase 2 (Zienkievicz):
|
||||
let rho =
|
||||
idx
|
||||
|> Array.map (fun i ->
|
||||
let r_im1 = access r (i - 1)
|
||||
let r_ip1 = access r (i + 1)
|
||||
(r[i] * r_ip1 * r_ip1 * r_im1 * r_im1)
|
||||
/ ((1.0 - r_ip1) * (1.0 - r_im1)))
|
||||
|
||||
// Phase 3:
|
||||
let fx_n = df_n |> Array.map fst
|
||||
let fy_n = df_n |> Array.map snd
|
||||
|
||||
let h_tilde =
|
||||
idx
|
||||
|> Array.map (fun i ->
|
||||
let fx_ip1 = access fx_n (i + 1)
|
||||
let fy_ip1 = access fy_n (i + 1)
|
||||
u[i] * fx_ip1 + v[i] * fy_ip1)
|
||||
|
||||
let k_tilde =
|
||||
idx
|
||||
|> Array.map (fun i ->
|
||||
let fx_im1 = access fx_n (i - 1)
|
||||
let fy_im1 = access fy_n (i - 1)
|
||||
u[i] * fx_im1 + v[i] * fy_im1)
|
||||
|
||||
let g_tilde =
|
||||
idx
|
||||
|> Array.map (fun i ->
|
||||
let r_ip1 = access r (i + 1)
|
||||
let r_im1 = access r (i - 1)
|
||||
let l2_ip1 = access l2 (i + 1)
|
||||
let l2_im1 = access l2 (i - 1)
|
||||
let rho_ip1 = access rho (i + 1)
|
||||
let rho_im1 = access rho (i - 1)
|
||||
(r_ip1 - r_im1) * phi[i]
|
||||
+ (3.0 * rho[i]) * (l2_ip1 - l2_im1) / l2[i]
|
||||
- rho_ip1
|
||||
+ rho_im1)
|
||||
|
||||
let w =
|
||||
idx
|
||||
|> Array.map (fun i ->
|
||||
let f_im1 = access f_n (i - 1)
|
||||
let f_ip1 = access f_n (i + 1)
|
||||
f_n[i] * r[i]
|
||||
+ 0.5 * (h_tilde[i] - k_tilde[i]) * phi[i]
|
||||
+ g_tilde[i]
|
||||
* (0.5 * (h_tilde[i] + k_tilde[i]) - f_im1 + f_ip1))
|
||||
w |> Array.sum
|
||||
@@ -1,332 +0,0 @@
|
||||
module Oceanbox.FvcomKit.Evaluate2D
|
||||
|
||||
open Grid
|
||||
open Gradient
|
||||
open CloughTocher
|
||||
open EvaluateVertex
|
||||
|
||||
let private (<*>) h (x: float, y: float) = (h * x, h * y)
|
||||
let private (<.>) (ax, ay) (bx, by) = (ax * bx + ay * by)
|
||||
let private (<+>) (ax, ay) (bx, by) = (ax + bx, ay + by)
|
||||
let private (<->) (ax, ay) (bx, by) = (ax - bx, ay - by)
|
||||
let private (><) (ax, ay) (bx, by) = (ax * by - ay * bx)
|
||||
|
||||
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 getCentroid (grid: ExtendedGrid) (idx: ElemIdx) = (grid.getCentroids ())[idx]
|
||||
|
||||
let private containsIndex (idx: NodeIdx) (cell: Cell) =
|
||||
let x, y, z = cell
|
||||
(idx = x) || (idx = y) || (idx = z)
|
||||
|
||||
let private countIndex (idx: NodeIdx) (cell: Cell) =
|
||||
match (containsIndex idx cell) with
|
||||
| true -> 1
|
||||
| false -> 0
|
||||
|
||||
let private countCommonIndices (c1: Cell) (c2: Cell) =
|
||||
let i0, i1, i2 = c1
|
||||
let s0 = countIndex i0 c2
|
||||
let s1 = countIndex i1 c2
|
||||
let s2 = countIndex i2 c2
|
||||
s0 + s1 + s2
|
||||
|
||||
// Elements are adjacent if they have a shared edge, i.e. have two common node indices
|
||||
let private areAdjacent (grid: ExtendedGrid) (i1: ElemIdx) (i0: ElemIdx) =
|
||||
let nCommon =
|
||||
countCommonIndices
|
||||
<| (getCell grid i0)
|
||||
<| (getCell grid i1)
|
||||
nCommon = 2
|
||||
|
||||
// Picks out the three neighbors of an element that are sharing a common edge
|
||||
// (two neighbors for boundary elements)
|
||||
let private getAdjacentNeighbors (grid: ExtendedGrid) (idx: ElemIdx) =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingElem idx
|
||||
|> Array.filter (areAdjacent grid idx)
|
||||
|
||||
// a valid centroid grid should always give two or three neighbors
|
||||
match neighbors.Length with
|
||||
| n when n < 2 -> failwith "Too few neighbors"
|
||||
| n when n > 3 -> failwith "Too many neighbors"
|
||||
| _ -> neighbors
|
||||
|
||||
let private isEdgeBoundary (grid: ExtendedGrid) (idx: ElemIdx) =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingElem idx
|
||||
|> Array.filter (areAdjacent grid idx)
|
||||
neighbors.Length = 2
|
||||
|
||||
let private isPointBoundary (grid: ExtendedGrid) (idx: ElemIdx) =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingElem idx
|
||||
|> Array.filter (areAdjacent grid idx)
|
||||
|
||||
match neighbors.Length with
|
||||
| 2 -> false
|
||||
| 3 ->
|
||||
let neighboundaries = neighbors |> Array.filter (isEdgeBoundary grid)
|
||||
neighboundaries.Length > 0
|
||||
| _ -> failwith "Invalid neighbor length"
|
||||
|
||||
let private getEdgeBoundaryVertices (grid: ExtendedGrid) (e: ElemIdx) : bool [] =
|
||||
let n =
|
||||
grid.getElemsSurroundingElem e
|
||||
|> Array.filter (areAdjacent grid e)
|
||||
match n.Length with
|
||||
| 2 ->
|
||||
let n0 = getCell grid n[0]
|
||||
let n1 = getCell grid n[1]
|
||||
let cell = (getCell grid e) |> toArray3
|
||||
cell
|
||||
|> Array.map (fun ci ->
|
||||
(not (containsIndex ci n0)
|
||||
|| not (containsIndex ci n1)))
|
||||
| _ -> failwith "Invalid number of edge boundary neighbors"
|
||||
|
||||
let private getEdges (grid: ExtendedGrid) (eIdx: ElemIdx) : Edge [] =
|
||||
let c0, c1, c2 = getCell grid eIdx
|
||||
[| (c0, c1); (c1, c2); (c2, c0) |]
|
||||
|
||||
let private pairEdge (e1: Edge) (e2: Edge) = (fst e1 = snd e2) && (snd e1 = fst e2)
|
||||
|
||||
let private sameEdge ((e11, e12): Edge) ((e21, e22): Edge) =
|
||||
(min e11 e12) = (min e21 e22)
|
||||
&& (max e11 e12) = (max e21 e22)
|
||||
|
||||
let rec private containsPair (edges: Edge list) (e: Edge) =
|
||||
match edges with
|
||||
| [] -> false
|
||||
| x :: xs -> if pairEdge e x then true else containsPair xs e
|
||||
|
||||
let rec private containsSame (edges: Edge list) (e: Edge) =
|
||||
match edges with
|
||||
| [] -> false
|
||||
| x :: xs -> if sameEdge e x then true else containsSame xs e
|
||||
|
||||
let private keepBoundary (edges: Edge list) =
|
||||
let rec loop O Y X =
|
||||
match X with
|
||||
| [] -> O
|
||||
| x :: xs ->
|
||||
if (containsPair Y x) then
|
||||
loop O Y xs
|
||||
else
|
||||
loop (x :: O) Y xs
|
||||
loop [] edges edges
|
||||
|
||||
let private commonEdges (e1: Edge list) (e2: Edge list) =
|
||||
let rec loop O Y X =
|
||||
match X with
|
||||
| [] -> O
|
||||
| x :: xs ->
|
||||
if (containsSame Y x) then
|
||||
loop (x :: O) Y xs
|
||||
else
|
||||
loop O Y xs
|
||||
loop [] e1 e2
|
||||
|
||||
let private toNodeList (edges: Edge list) =
|
||||
let rec append O X =
|
||||
match X with
|
||||
| [] -> O
|
||||
| (x1, x2) :: xs -> append (x1 :: x2 :: O) xs
|
||||
append [] edges
|
||||
|
||||
let private onGlobalBoundary (grid: ExtendedGrid) (idx: NodeIdx) : bool =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> toNodeList
|
||||
|> List.contains idx
|
||||
|
||||
let private contains (idx: NodeIdx) (edge: Edge) : bool = fst edge = idx || snd edge = idx
|
||||
|
||||
let private projectOnEdge (grid: ExtendedGrid) (v: float * float) (edge: Edge) =
|
||||
let n0, n1 = edge
|
||||
let x0, y0 = getVertex grid n0
|
||||
let x1, y1 = getVertex grid n1
|
||||
let u = (x1 - x0, y1 - y0)
|
||||
|
||||
// Project only if pointing outward
|
||||
if (v >< u) < 0.0 then (v <.> u) / (u <.> u) <*> u else v
|
||||
|
||||
let private divertOnEdge (grid: ExtendedGrid) (v: float * float) (edge: Edge) =
|
||||
let n0, n1 = edge
|
||||
let x0, y0 = getVertex grid n0
|
||||
let x1, y1 = getVertex grid n1
|
||||
let u = (x1 - x0, y1 - y0)
|
||||
|
||||
// Project only if pointing outward
|
||||
if (v >< u) < 0.0 then
|
||||
let sign = (u <.> v) / abs (u <.> v)
|
||||
sign * sqrt ((v <.> v) / (u <.> u)) <*> u
|
||||
else
|
||||
v
|
||||
|
||||
let private projectOnBoundary (grid: ExtendedGrid) (eIdx: ElemIdx) (nIdx: NodeIdx, v: float * float) =
|
||||
let boundaries =
|
||||
grid.getElemsSurroundingElem eIdx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> commonEdges (getEdges grid eIdx |> Array.toList)
|
||||
|
||||
match boundaries.Length with
|
||||
| 0 -> v
|
||||
| 1 ->
|
||||
let b = boundaries[0]
|
||||
if contains nIdx b then projectOnEdge grid v b else v
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private divertOnBoundary (grid: ExtendedGrid) (eIdx: ElemIdx) (nIdx: NodeIdx, v: float * float) =
|
||||
let boundaries =
|
||||
grid.getElemsSurroundingElem eIdx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> commonEdges (getEdges grid eIdx |> Array.toList)
|
||||
|
||||
match boundaries.Length with
|
||||
| 0 -> v
|
||||
| 1 ->
|
||||
let b = boundaries[0]
|
||||
if contains nIdx b then divertOnEdge grid v b else v
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private projectOnNeighboundaries (grid: ExtendedGrid) (idx: NodeIdx, v: float * float) =
|
||||
let neighboundaries =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> List.filter (contains idx)
|
||||
|
||||
let v_proj = neighboundaries |> List.map (projectOnEdge grid v)
|
||||
match v_proj.Length with
|
||||
| 0 -> v
|
||||
| 1 -> v_proj[0]
|
||||
| 2 -> 0.5 <*> (v_proj[0] <+> v_proj[1])
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private divertOnNeighboundaries (grid: ExtendedGrid) (idx: NodeIdx, v: float * float) =
|
||||
let neighboundaries =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> List.filter (contains idx)
|
||||
|
||||
let v_proj = neighboundaries |> List.map (divertOnEdge grid v)
|
||||
match v_proj.Length with
|
||||
| 0 -> v
|
||||
| 1 -> v_proj[0]
|
||||
| 2 -> 0.5 <*> (v_proj[0] <+> v_proj[1])
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private adjustBoundary (grid: ExtendedGrid) (e: ElemIdx) (uv_n: (float * float) []) =
|
||||
let idx_n = (getCell grid e) |> toArray3
|
||||
if (isEdgeBoundary grid e) then
|
||||
(Array.zip idx_n uv_n)
|
||||
|> Array.map (divertOnBoundary grid e)
|
||||
elif (isPointBoundary grid e) then
|
||||
(Array.zip idx_n uv_n)
|
||||
|> Array.map (divertOnNeighboundaries grid)
|
||||
else
|
||||
uv_n
|
||||
|
||||
// Evaluate 2D field based on constant extrapolation from centroid field value
|
||||
let evaluateConstantCentroid2D (grid: ExtendedGrid) (readUV: ElemIdx -> float * float) ((px, py): Pos) =
|
||||
match grid.tryGetElement (px, py) with
|
||||
| None -> None
|
||||
| Some e -> Some(readUV e)
|
||||
|
||||
// Evaluate 2D field based on linear extrapolation from centroid field value,
|
||||
// using gradient estimated from values in the three adjacent neigboring elements
|
||||
let evaluateLinearCentroid2D (grid: ExtendedGrid) (readUV: ElemIdx -> float * float) ((px, py): Pos) =
|
||||
match grid.tryGetElement (px, py) with
|
||||
| None -> None
|
||||
| Some e ->
|
||||
let centroids = grid.getCentroids ()
|
||||
|
||||
let c0 = centroids[e]
|
||||
let u0 = readUV e |> fst
|
||||
let v0 = readUV e |> snd
|
||||
|
||||
let neighbors = getAdjacentNeighbors grid e
|
||||
let cn = neighbors |> Array.map (Array.get centroids)
|
||||
let un =
|
||||
neighbors
|
||||
|> Array.map (fun en -> readUV en |> fst)
|
||||
let vn =
|
||||
neighbors
|
||||
|> Array.map (fun en -> readUV en |> snd)
|
||||
|
||||
let thrs = 1.0e-3
|
||||
let du_x, du_y = calcGradientFromNeighbors c0 u0 cn un thrs
|
||||
let dv_x, dv_y = calcGradientFromNeighbors c0 v0 cn vn thrs
|
||||
|
||||
let x, y = px, py
|
||||
let x0, y0 = c0
|
||||
|
||||
let u = u0 + du_x * (x - x0) + du_y * (y - y0)
|
||||
let v = v0 + dv_x * (x - x0) + dv_y * (y - y0)
|
||||
Some(u, v)
|
||||
|
||||
// Evaluate 2D field based on Clough-Tocher interpolation within grid element
|
||||
let evaluateCloughTocherCentroid2D (grid: ExtendedGrid) (readUV: ElemIdx -> float * float) ((px, py): Pos) =
|
||||
match grid.tryGetElement (px, py) with
|
||||
| None -> None
|
||||
| Some e ->
|
||||
let idx = (getCell grid e) |> toArray3
|
||||
|
||||
let readU = fst << readUV
|
||||
let readV = snd << readUV
|
||||
|
||||
// Position of vertices
|
||||
let p_n = idx |> Array.map (getVertex grid)
|
||||
|
||||
// Field values at vertices
|
||||
let u_n = idx |> Array.map (calcVertexValue1D grid readU)
|
||||
let v_n = idx |> Array.map (calcVertexValue1D grid readV)
|
||||
|
||||
// Enforce boundary conditions
|
||||
let uv_n = Array.zip u_n v_n
|
||||
let uv_n' = adjustBoundary grid e uv_n
|
||||
let u_n' = uv_n' |> Array.map fst
|
||||
let v_n' = uv_n' |> Array.map snd
|
||||
|
||||
// Gradient values at vertices
|
||||
let thrs = 1.0e-3
|
||||
let du_n =
|
||||
(Array.zip idx u_n')
|
||||
|> Array.map (calcVertexGradient grid readU thrs)
|
||||
let dv_n =
|
||||
(Array.zip idx v_n')
|
||||
|> Array.map (calcVertexGradient grid readV thrs)
|
||||
|
||||
// Rearrange to x and y derivatives: [grad(U), grad(V)] -> [dx(UV), dy(UV)]
|
||||
let dx_uv_n = Array.zip (Array.map fst du_n) (Array.map fst dv_n)
|
||||
let dy_uv_n = Array.zip (Array.map snd du_n) (Array.map snd dv_n)
|
||||
|
||||
// Enforce boundary conditions
|
||||
let dx_uv_n' = adjustBoundary grid e dx_uv_n
|
||||
let dy_uv_n' = adjustBoundary grid e dy_uv_n
|
||||
|
||||
// Rearrange back to gradients: [dx(UV), dy(UV)] -> [grad(U), grad(V)]
|
||||
let du_n' = Array.zip (Array.map fst dx_uv_n') (Array.map fst dy_uv_n')
|
||||
let dv_n' = Array.zip (Array.map snd dx_uv_n') (Array.map snd dy_uv_n')
|
||||
|
||||
// Perform Clough-Tocher interpolation
|
||||
let u1 = interpolateTriangleCT p_n u_n' du_n' (px, py)
|
||||
let v1 = interpolateTriangleCT p_n v_n' dv_n' (px, py)
|
||||
Some(u1, v1)
|
||||
@@ -1,496 +0,0 @@
|
||||
module Oceanbox.FvcomKit.Evaluate4D
|
||||
|
||||
open Grid
|
||||
open Gradient
|
||||
open CloughTocher
|
||||
open EvaluateVertex
|
||||
|
||||
let private (<*>) h (x: float, y: float, z: float) = (h * x, h * y, h * z)
|
||||
let private (<.>) (ax, ay, az) (bx, by, bz) = (ax * bx + ay * by + az * bz)
|
||||
let private (<+>) (ax, ay, az) (bx, by, bz) = (ax + bx, ay + by, az + bz)
|
||||
let private (<->) (ax, ay, az) (bx, by, bz) = (ax - bx, ay - by, az - bz)
|
||||
|
||||
let private (><) (ax, ay, az) (bx, by, bz) =
|
||||
(ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
|
||||
|
||||
let private fst3 (a, _, _) = a
|
||||
let private snd3 (_, b, _) = b
|
||||
let private trd3 (_, _, c) = c
|
||||
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 getCentroid (grid: ExtendedGrid) (idx: ElemIdx) = (grid.getCentroids ())[idx]
|
||||
|
||||
let private containsIndex (idx: NodeIdx) (cell: Cell) =
|
||||
let x, y, z = cell
|
||||
(idx = x) || (idx = y) || (idx = z)
|
||||
|
||||
let private countIndex (idx: NodeIdx) (cell: Cell) =
|
||||
match (containsIndex idx cell) with
|
||||
| true -> 1
|
||||
| false -> 0
|
||||
|
||||
let private countCommonIndices (c1: Cell) (c2: Cell) =
|
||||
let i0, i1, i2 = c1
|
||||
let s0 = countIndex i0 c2
|
||||
let s1 = countIndex i1 c2
|
||||
let s2 = countIndex i2 c2
|
||||
s0 + s1 + s2
|
||||
|
||||
// Elements are adjacent if they have a shared edge, i.e. have two common node indices
|
||||
let private areAdjacent (grid: ExtendedGrid) (i1: ElemIdx) (i0: ElemIdx) =
|
||||
let nCommon =
|
||||
countCommonIndices
|
||||
<| (getCell grid i0)
|
||||
<| (getCell grid i1)
|
||||
nCommon = 2
|
||||
|
||||
// Picks out the three neighbors of an element that are sharing a common edge
|
||||
// (two neighbors for boundary elements)
|
||||
let private getAdjacentNeighbors (grid: ExtendedGrid) (idx: ElemIdx) =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingElem idx
|
||||
|> Array.filter (areAdjacent grid idx)
|
||||
|
||||
// a valid centroid grid should always give two or three neighbors
|
||||
match neighbors.Length with
|
||||
| n when n < 2 -> failwith "Too few neighbors"
|
||||
| n when n > 3 -> failwith "Too many neighbors"
|
||||
| _ -> neighbors
|
||||
|
||||
let private isEdgeBoundary (grid: ExtendedGrid) (idx: ElemIdx) =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingElem idx
|
||||
|> Array.filter (areAdjacent grid idx)
|
||||
neighbors.Length = 2
|
||||
|
||||
let private isPointBoundary (grid: ExtendedGrid) (idx: ElemIdx) =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingElem idx
|
||||
|> Array.filter (areAdjacent grid idx)
|
||||
|
||||
match neighbors.Length with
|
||||
| 2 -> false
|
||||
| 3 ->
|
||||
let neighboundaries = neighbors |> Array.filter (isEdgeBoundary grid)
|
||||
neighboundaries.Length > 0
|
||||
| _ -> failwith "Invalid neighbor length"
|
||||
|
||||
let private getEdgeBoundaryVertices (grid: ExtendedGrid) (e: ElemIdx) : bool [] =
|
||||
let n =
|
||||
grid.getElemsSurroundingElem e
|
||||
|> Array.filter (areAdjacent grid e)
|
||||
match n.Length with
|
||||
| 2 ->
|
||||
let n0 = getCell grid n[0]
|
||||
let n1 = getCell grid n[1]
|
||||
let cell = (getCell grid e) |> toArray3
|
||||
cell
|
||||
|> Array.map (fun ci ->
|
||||
(not (containsIndex ci n0)
|
||||
|| not (containsIndex ci n1)))
|
||||
| _ -> failwith "Invalid number of edge boundary neighbors"
|
||||
|
||||
let private getEdges (grid: ExtendedGrid) (eIdx: ElemIdx) : Edge [] =
|
||||
let c0, c1, c2 = getCell grid eIdx
|
||||
[| (c0, c1); (c1, c2); (c2, c0) |]
|
||||
|
||||
let private pairEdge (e1: Edge) (e2: Edge) = (fst e1 = snd e2) && (snd e1 = fst e2)
|
||||
|
||||
let private sameEdge ((e11, e12): Edge) ((e21, e22): Edge) =
|
||||
(min e11 e12) = (min e21 e22)
|
||||
&& (max e11 e12) = (max e21 e22)
|
||||
|
||||
let rec private containsPair (edges: Edge list) (e: Edge) =
|
||||
match edges with
|
||||
| [] -> false
|
||||
| x :: xs -> if pairEdge e x then true else containsPair xs e
|
||||
|
||||
let rec private containsSame (edges: Edge list) (e: Edge) =
|
||||
match edges with
|
||||
| [] -> false
|
||||
| x :: xs -> if sameEdge e x then true else containsSame xs e
|
||||
|
||||
let private keepBoundary (edges: Edge list) =
|
||||
let rec loop O Y X =
|
||||
match X with
|
||||
| [] -> O
|
||||
| x :: xs ->
|
||||
if (containsPair Y x) then
|
||||
loop O Y xs
|
||||
else
|
||||
loop (x :: O) Y xs
|
||||
loop [] edges edges
|
||||
|
||||
let private commonEdges (e1: Edge list) (e2: Edge list) =
|
||||
let rec loop O Y X =
|
||||
match X with
|
||||
| [] -> O
|
||||
| x :: xs ->
|
||||
if (containsSame Y x) then
|
||||
loop (x :: O) Y xs
|
||||
else
|
||||
loop O Y xs
|
||||
loop [] e1 e2
|
||||
|
||||
let private toNodeList (edges: Edge list) =
|
||||
let rec append O X =
|
||||
match X with
|
||||
| [] -> O
|
||||
| (x1, x2) :: xs -> append (x1 :: x2 :: O) xs
|
||||
append [] edges
|
||||
|
||||
let private onGlobalBoundary (grid: ExtendedGrid) (idx: NodeIdx) : bool =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> toNodeList
|
||||
|> List.contains idx
|
||||
|
||||
let private contains (idx: NodeIdx) (edge: Edge) : bool = fst edge = idx || snd edge = idx
|
||||
|
||||
let private projectOnEdge (grid: ExtendedGrid) (v: float * float * float) (edge: Edge) : float * float * float =
|
||||
let n0, n1 = edge
|
||||
let x0, y0 = getVertex grid n0
|
||||
let x1, y1 = getVertex grid n1
|
||||
let u = (x1 - x0, y1 - y0, 0.0)
|
||||
|
||||
// Project only if pointing outward (negative z-component of cross product)
|
||||
let _, _, k = v >< u
|
||||
if k < 0.0 then (v <.> u) / (u <.> u) <*> u else v
|
||||
|
||||
let private divertOnEdge (grid: ExtendedGrid) (v: float * float * float) (edge: Edge) : float * float * float =
|
||||
let n0, n1 = edge
|
||||
let x0, y0 = getVertex grid n0
|
||||
let x1, y1 = getVertex grid n1
|
||||
let u = (x1 - x0, y1 - y0, 0.0)
|
||||
|
||||
// Project only if pointing outward (negative z-component of cross product)
|
||||
let _, _, k = v >< u
|
||||
if k < 0.0 then
|
||||
let sign = (u <.> v) / abs (u <.> v)
|
||||
sign * sqrt ((v <.> v) / (u <.> u)) <*> u
|
||||
else
|
||||
v
|
||||
|
||||
let private projectOnBoundary
|
||||
(grid: ExtendedGrid)
|
||||
(eIdx: ElemIdx)
|
||||
(nIdx: NodeIdx, v: float * float * float)
|
||||
: float * float * float =
|
||||
let boundaries =
|
||||
grid.getElemsSurroundingElem eIdx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> commonEdges (getEdges grid eIdx |> Array.toList)
|
||||
|
||||
match boundaries.Length with
|
||||
| 0 -> v
|
||||
| 1 ->
|
||||
let b = boundaries[0]
|
||||
if contains nIdx b then projectOnEdge grid v b else v
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private divertOnBoundary
|
||||
(grid: ExtendedGrid)
|
||||
(eIdx: ElemIdx)
|
||||
(nIdx: NodeIdx, v: float * float * float)
|
||||
: float * float * float =
|
||||
let boundaries =
|
||||
grid.getElemsSurroundingElem eIdx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> commonEdges (getEdges grid eIdx |> Array.toList)
|
||||
|
||||
match boundaries.Length with
|
||||
| 0 -> v
|
||||
| 1 ->
|
||||
let b = boundaries[0]
|
||||
if contains nIdx b then divertOnEdge grid v b else v
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private projectOnNeighboundaries
|
||||
(grid: ExtendedGrid)
|
||||
(idx: NodeIdx, v: float * float * float)
|
||||
: float * float * float =
|
||||
let neighboundaries =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> List.filter (contains idx)
|
||||
|
||||
let v_proj = neighboundaries |> List.map (projectOnEdge grid v)
|
||||
match v_proj.Length with
|
||||
| 0 -> v
|
||||
| 1 -> v_proj[0]
|
||||
| 2 -> 0.5 <*> (v_proj[0] <+> v_proj[1])
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private divertOnNeighboundaries
|
||||
(grid: ExtendedGrid)
|
||||
(idx: NodeIdx, v: float * float * float)
|
||||
: float * float * float =
|
||||
let neighboundaries =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.map (getEdges grid)
|
||||
|> Array.reduce Array.append
|
||||
|> Array.toList
|
||||
|> keepBoundary
|
||||
|> List.filter (contains idx)
|
||||
|
||||
let v_proj = neighboundaries |> List.map (divertOnEdge grid v)
|
||||
match v_proj.Length with
|
||||
| 0 -> v
|
||||
| 1 -> v_proj[0]
|
||||
| 2 -> 0.5 <*> (v_proj[0] <+> v_proj[1])
|
||||
| _ -> failwith "Too many boundary edges"
|
||||
|
||||
let private adjustBoundary (grid: ExtendedGrid) (e: ElemIdx) (uvw_n: (float * float * float) []) =
|
||||
let idx_n = (getCell grid e) |> toArray3
|
||||
if (isEdgeBoundary grid e) then
|
||||
(Array.zip idx_n uvw_n)
|
||||
|> Array.map (divertOnBoundary grid e)
|
||||
elif (isPointBoundary grid e) then
|
||||
(Array.zip idx_n uvw_n)
|
||||
|> Array.map (divertOnNeighboundaries grid)
|
||||
else
|
||||
uvw_n
|
||||
|
||||
let private calcSigmaLayers (field: Field4D) (sigma: float) : int * int =
|
||||
let slayers = field.Siglays
|
||||
let slayer =
|
||||
slayers
|
||||
|> Array.filter (fun s -> (float s) > sigma)
|
||||
|> Array.length
|
||||
match slayer with
|
||||
| s when s = 0 -> (0, 0)
|
||||
| s when s = slayers.Length -> (s - 1, s - 1)
|
||||
| _ -> (slayer - 1, slayer)
|
||||
|
||||
// Linear interpolation in sigma and time in given cell
|
||||
let private interpolateST (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
|
||||
let t0 = field.Time0
|
||||
let t1 = field.Time1
|
||||
let frame0 = field.Frame0
|
||||
let frame1 = field.Frame1
|
||||
|
||||
// Interpolate f(x) between fi = f(xi) and fj = f(xj)
|
||||
let fx x i j xi xj fi fj =
|
||||
if (i = j) then
|
||||
fi
|
||||
else
|
||||
let a = (x - xi) / (xj - xi)
|
||||
(a * fj) + (1.0 - a) * fi
|
||||
|
||||
let interpolateSigma fi fj = fx s slayer0 slayer1 s0 s1 fi fj
|
||||
let interpolateTime fi fj = fx t frame0 frame1 t0 t1 fi fj
|
||||
|
||||
let u00, v00, w00 = field.Field0[slayer0][eIdx]
|
||||
let u01, v01, w01 = field.Field0[slayer1][eIdx]
|
||||
let u10, v10, w10 = field.Field1[slayer0][eIdx]
|
||||
let u11, v11, w11 = field.Field1[slayer1][eIdx]
|
||||
|
||||
let u0 = interpolateTime (float u00) (float u10)
|
||||
let u1 = interpolateTime (float u01) (float u11)
|
||||
let u = interpolateSigma u0 u1
|
||||
|
||||
let v0 = interpolateTime (float v00) (float v10)
|
||||
let v1 = interpolateTime (float v01) (float v11)
|
||||
let v = interpolateSigma v0 v1
|
||||
|
||||
let w0 = interpolateTime (float w00) (float w10)
|
||||
let w1 = interpolateTime (float w01) (float w11)
|
||||
let w = interpolateSigma w0 w1
|
||||
|
||||
let h = field.Bath[eIdx] |> float
|
||||
(u, v, w / h)
|
||||
|
||||
// Evaluate 4D field based on constant horizontal extrapolation from centroid field value
|
||||
// Linear interpolation in vertical and time
|
||||
let evaluateConstantCentroid4D (grid: ExtendedGrid) (field: Field4D) (t: float) ((x, y, s): Pos3D) =
|
||||
match grid.tryGetElement (x, y) with
|
||||
| None -> None
|
||||
| Some e ->
|
||||
// Set up linear sigma and time interpolation
|
||||
let readUVW eIdx = interpolateST field s t eIdx
|
||||
|
||||
let u = readUVW e |> fst3
|
||||
let v = readUVW e |> snd3
|
||||
let w = readUVW e |> trd3
|
||||
|
||||
Some(u, v, w)
|
||||
|
||||
// Evaluate 4D field based on linear horizontal extrapolation from centroid field value,
|
||||
// using gradient estimated from values in the three adjacent neigboring elements.
|
||||
// Linear interpolation in vertical and time
|
||||
let evaluateLinearCentroid4D (grid: ExtendedGrid) (field: Field4D) (t: float) ((x, y, s): Pos3D) =
|
||||
match grid.tryGetElement (x, y) with
|
||||
| None -> None
|
||||
| Some e ->
|
||||
// Set up linear sigma and time interpolation
|
||||
let readUVW eIdx = interpolateST field s t eIdx
|
||||
|
||||
let centroids = grid.getCentroids ()
|
||||
let c0 = centroids[e]
|
||||
let u0 = readUVW e |> fst3
|
||||
let v0 = readUVW e |> snd3
|
||||
let w0 = readUVW e |> trd3
|
||||
|
||||
let neighbors = getAdjacentNeighbors grid e
|
||||
let cn = neighbors |> Array.map (Array.get centroids)
|
||||
let un =
|
||||
neighbors
|
||||
|> Array.map (fun en -> readUVW en |> fst3)
|
||||
let vn =
|
||||
neighbors
|
||||
|> Array.map (fun en -> readUVW en |> snd3)
|
||||
let wn =
|
||||
neighbors
|
||||
|> Array.map (fun en -> readUVW en |> trd3)
|
||||
|
||||
let thrs = 1.0e-6
|
||||
let dxU, dyU = calcGradientFromNeighbors c0 u0 cn un thrs
|
||||
let dxV, dyV = calcGradientFromNeighbors c0 v0 cn vn thrs
|
||||
let dxW, dyW = calcGradientFromNeighbors c0 w0 cn wn thrs
|
||||
|
||||
let dx = x - (fst c0)
|
||||
let dy = y - (snd c0)
|
||||
|
||||
let u = u0 + dxU * dx + dyU * dy
|
||||
let v = v0 + dxV * dx + dyV * dy
|
||||
let w = w0 + dxW * dx + dyW * dy
|
||||
Some(u, v, w)
|
||||
|
||||
// Evaluate 4D field based on horizontal Clough-Tocher interpolation within grid element
|
||||
// Linear interpolation in vertical and time
|
||||
// Expects field to be given in centroid points
|
||||
let evaluateCloughTocherCentroid4D (grid: ExtendedGrid) (field: Field4D) (t: float) ((x, y, s): Pos3D) =
|
||||
match grid.tryGetElement (x, y) with
|
||||
| None -> None
|
||||
| Some e ->
|
||||
// Set up linear sigma and time interpolation
|
||||
let readUVW eIdx = interpolateST field s t eIdx
|
||||
|
||||
// Position of vertices
|
||||
let idx = (getCell grid e) |> toArray3
|
||||
let p_n = idx |> Array.map (getVertex grid)
|
||||
|
||||
// Field values at vertices
|
||||
let uvw_n = idx |> Array.map (calcVertexValue3D grid readUVW)
|
||||
let uvw_n' = adjustBoundary grid e uvw_n
|
||||
let u_n' = uvw_n' |> Array.map fst3
|
||||
let v_n' = uvw_n' |> Array.map snd3
|
||||
let w_n' = uvw_n' |> Array.map trd3
|
||||
|
||||
let readU = fst3 << readUVW
|
||||
let readV = snd3 << readUVW
|
||||
let readW = trd3 << readUVW
|
||||
|
||||
// Gradient values at vertices
|
||||
let thrs = 1.0e-6
|
||||
let du_n =
|
||||
(Array.zip idx u_n')
|
||||
|> Array.map (calcVertexGradient grid readU thrs)
|
||||
let dv_n =
|
||||
(Array.zip idx v_n')
|
||||
|> Array.map (calcVertexGradient grid readV thrs)
|
||||
let dw_n =
|
||||
(Array.zip idx w_n')
|
||||
|> Array.map (calcVertexGradient grid readW thrs)
|
||||
|
||||
// Rearrange to x and y derivatives: [grad(U), grad(V)] -> [dx(UV), dy(UV)]
|
||||
let dx_uvw_n =
|
||||
Array.zip3 (Array.map fst du_n) (Array.map fst dv_n) (Array.map fst dw_n)
|
||||
let dy_uvw_n =
|
||||
Array.zip3 (Array.map snd du_n) (Array.map snd dv_n) (Array.map snd dw_n)
|
||||
|
||||
// Enforce boundary conditions
|
||||
let dx_uvw_n' = adjustBoundary grid e dx_uvw_n
|
||||
let dy_uvw_n' = adjustBoundary grid e dy_uvw_n
|
||||
|
||||
// Rearrange back to gradients: [dx(UV), dy(UV)] -> [grad(U), grad(V)]
|
||||
let du_n' =
|
||||
Array.zip (Array.map fst3 dx_uvw_n') (Array.map fst3 dy_uvw_n')
|
||||
let dv_n' =
|
||||
Array.zip (Array.map snd3 dx_uvw_n') (Array.map snd3 dy_uvw_n')
|
||||
let dw_n' =
|
||||
Array.zip (Array.map trd3 dx_uvw_n') (Array.map trd3 dy_uvw_n')
|
||||
|
||||
// Perform Clough-Tocher interpolation
|
||||
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)
|
||||
|
||||
|
||||
// Evaluate 4D field based on Clough-Tocher interpolation within grid element
|
||||
// Expects field to be given in vertex points
|
||||
// Linear interpolation in vertical and time
|
||||
let evaluateCloughTocherVertex4D (grid: ExtendedGrid) (field: Field4D) (t: float) ((x, y, s): Pos3D) =
|
||||
match grid.tryGetElement (x, y) with
|
||||
| None -> None
|
||||
| Some e ->
|
||||
// Set up linear sigma and time interpolation
|
||||
let readUVW eIdx = interpolateST field s t eIdx
|
||||
let readU = fst3 << readUVW
|
||||
let readV = snd3 << readUVW
|
||||
let readW = trd3 << readUVW
|
||||
|
||||
// Position of vertices
|
||||
let idx = (getCell grid e) |> toArray3
|
||||
let p_n = idx |> Array.map (getVertex grid)
|
||||
|
||||
// Field values at vertices
|
||||
let uvw_n = idx |> Array.map readUVW
|
||||
let uvw_n' = adjustBoundary grid e uvw_n
|
||||
let u_n' = uvw_n' |> Array.map fst3
|
||||
let v_n' = uvw_n' |> Array.map snd3
|
||||
let w_n' = uvw_n' |> Array.map trd3
|
||||
|
||||
// Gradient values at vertices
|
||||
let thrs = 1.0e-6
|
||||
let du_n =
|
||||
(Array.zip idx u_n')
|
||||
|> Array.map (calcVertexGradient' grid readU thrs)
|
||||
let dv_n =
|
||||
(Array.zip idx v_n')
|
||||
|> Array.map (calcVertexGradient' grid readV thrs)
|
||||
let dw_n =
|
||||
(Array.zip idx w_n')
|
||||
|> Array.map (calcVertexGradient' grid readW thrs)
|
||||
|
||||
// Rearrange to x and y derivatives: [grad(U), grad(V)] -> [dx(UV), dy(UV)]
|
||||
let dx_uvw_n =
|
||||
Array.zip3 (Array.map fst du_n) (Array.map fst dv_n) (Array.map fst dw_n)
|
||||
let dy_uvw_n =
|
||||
Array.zip3 (Array.map snd du_n) (Array.map snd dv_n) (Array.map snd dw_n)
|
||||
|
||||
// Enforce boundary conditions
|
||||
let dx_uvw_n' = adjustBoundary grid e dx_uvw_n
|
||||
let dy_uvw_n' = adjustBoundary grid e dy_uvw_n
|
||||
|
||||
// Rearrange back to gradients: [dx(UV), dy(UV)] -> [grad(U), grad(V)]
|
||||
let du_n' =
|
||||
Array.zip (Array.map fst3 dx_uvw_n') (Array.map fst3 dy_uvw_n')
|
||||
let dv_n' =
|
||||
Array.zip (Array.map snd3 dx_uvw_n') (Array.map snd3 dy_uvw_n')
|
||||
let dw_n' =
|
||||
Array.zip (Array.map trd3 dx_uvw_n') (Array.map trd3 dy_uvw_n')
|
||||
|
||||
// Perform Clough-Tocher interpolation
|
||||
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)
|
||||
@@ -1,109 +0,0 @@
|
||||
module Oceanbox.FvcomKit.EvaluateVertex
|
||||
|
||||
open Grid
|
||||
|
||||
let private (<->) (ax, ay) (bx, by) = (ax - bx, ay - by)
|
||||
let private inv_norm (px, py) = 1.0 / sqrt (px * px + py * py)
|
||||
|
||||
let private fst3 (u, _, _) = u
|
||||
let private snd3 (_, v, _) = v
|
||||
let private trd3 (_, _, w) = w
|
||||
|
||||
let private getVertex (grid: ExtendedGrid) (idx: NodeIdx) = (grid :> IGrid).getVertex idx
|
||||
let private getCentroid (grid: ExtendedGrid) (idx: ElemIdx) = (grid.getCentroids ())[idx]
|
||||
|
||||
let calcVertexValue1D (grid: ExtendedGrid) (readU: ElemIdx -> float) (idx: NodeIdx) : float =
|
||||
// Fetch surrounding elements
|
||||
let e_n = grid.getElemsSurroundingNode idx
|
||||
|
||||
// Fetch points
|
||||
let p_0 = getVertex grid idx
|
||||
let p_n = e_n |> Array.map (getCentroid grid)
|
||||
|
||||
// Fetch function values
|
||||
let f_n = e_n |> Array.map readU
|
||||
|
||||
// Compute inverse distance weights
|
||||
let w_n =
|
||||
p_n
|
||||
|> Array.map (fun p -> p <-> p_0)
|
||||
|> Array.map inv_norm
|
||||
let W = w_n |> Array.sum
|
||||
|
||||
// Compute mean function value on vertex P_0
|
||||
(Array.map2 (*) w_n f_n |> Array.sum) / W
|
||||
|
||||
let centroidToVertexValues1D (grid: ExtendedGrid) (readU: ElemIdx -> float) : float [] =
|
||||
let N = ((grid :> IGrid).getVertices ()).Length
|
||||
[| 0 .. (N - 1) |]
|
||||
|> Array.Parallel.map (calcVertexValue1D grid readU)
|
||||
|
||||
|
||||
|
||||
// Given 2D centroid evaluator "readUV", compute field value in single vertex
|
||||
// by a weighted average of values in all surrounding elements (within sigma layer)
|
||||
let calcVertexValue2D (grid: ExtendedGrid) (readUV: ElemIdx -> float * float) (idx: NodeIdx) =
|
||||
// Fetch surrounding elements
|
||||
let e_n = grid.getElemsSurroundingNode idx
|
||||
|
||||
// Fetch points
|
||||
let p_0 = getVertex grid idx
|
||||
let p_n = e_n |> Array.map (getCentroid grid)
|
||||
|
||||
// Fetch function values
|
||||
let f_n = e_n |> Array.map readUV
|
||||
let u_n = f_n |> Array.map fst
|
||||
let v_n = f_n |> Array.map snd
|
||||
|
||||
// Compute inverse distance weights
|
||||
let wgt_n =
|
||||
p_n
|
||||
|> Array.map (fun p -> p <-> p_0)
|
||||
|> Array.map inv_norm
|
||||
let W = wgt_n |> Array.sum
|
||||
|
||||
// Compute mean function value on vertex P_0
|
||||
let u_0 = (Array.map2 (*) wgt_n u_n |> Array.sum) / W
|
||||
let v_0 = (Array.map2 (*) wgt_n v_n |> Array.sum) / W
|
||||
u_0, v_0
|
||||
|
||||
let centroidToVertexValues2D (grid: ExtendedGrid) (readUV: ElemIdx -> float * float) =
|
||||
let N = ((grid :> IGrid).getVertices ()).Length
|
||||
[| 0 .. (N - 1) |]
|
||||
|> Array.Parallel.map (calcVertexValue2D grid readUV)
|
||||
|
||||
|
||||
|
||||
// Given 3D centroid evaluator "readUVW", compute field value in single vertex
|
||||
// by a weighted average of values in all surrounding elements (within sigma layer)
|
||||
let calcVertexValue3D (grid: ExtendedGrid) (readUVW: ElemIdx -> float * float * float) (idx: NodeIdx) =
|
||||
// Fetch surrounding elements
|
||||
let e_n = grid.getElemsSurroundingNode idx
|
||||
|
||||
// Fetch points
|
||||
let p_0 = getVertex grid idx
|
||||
let p_n = e_n |> Array.map (getCentroid grid)
|
||||
|
||||
// Fetch function values
|
||||
let f_n = e_n |> Array.map readUVW
|
||||
let u_n = f_n |> Array.map fst3
|
||||
let v_n = f_n |> Array.map snd3
|
||||
let w_n = f_n |> Array.map trd3
|
||||
|
||||
// Compute inverse distance weights
|
||||
let wgt_n =
|
||||
p_n
|
||||
|> Array.map (fun p -> p <-> p_0)
|
||||
|> Array.map inv_norm
|
||||
let W = wgt_n |> Array.sum
|
||||
|
||||
// Compute mean function value on vertex P_0
|
||||
let u_0 = (Array.map2 (*) wgt_n u_n |> Array.sum) / W
|
||||
let v_0 = (Array.map2 (*) wgt_n v_n |> Array.sum) / W
|
||||
let w_0 = (Array.map2 (*) wgt_n w_n |> Array.sum) / W
|
||||
u_0, v_0, w_0
|
||||
|
||||
let centroidToVertexValues3D (grid: ExtendedGrid) (readUVW: ElemIdx -> float * float * float) =
|
||||
let N = ((grid :> IGrid).getVertices ()).Length
|
||||
[| 0 .. (N - 1) |]
|
||||
|> Array.Parallel.map (calcVertexValue3D grid readUVW)
|
||||
122
src/Gradient.fs
122
src/Gradient.fs
@@ -1,122 +0,0 @@
|
||||
module Oceanbox.FvcomKit.Gradient
|
||||
|
||||
open Grid
|
||||
|
||||
let private getCell (grid: ExtendedGrid) (idx: ElemIdx) = (grid :> IGrid).getCell idx
|
||||
let private getVertex (grid: ExtendedGrid) (idx: NodeIdx) = (grid :> IGrid).getVertex idx
|
||||
let private getCentroid (grid: ExtendedGrid) (idx: ElemIdx) = (grid.getCentroids ())[idx]
|
||||
|
||||
// Compute contributions to Q and s from a single neighboring point, given current estimate of gradient df
|
||||
let private calcContributionQandS (df: float * float) (p0: Pos) (f0: float) (p1: Pos) (f1: float) =
|
||||
let x0, y0 = p0
|
||||
let x1, y1 = p1
|
||||
|
||||
let ex = x1 - x0
|
||||
let ey = y1 - y0
|
||||
|
||||
let L = sqrt (ex * ex + ey * ey)
|
||||
let L3 = L * L * L
|
||||
|
||||
let dfx, dfy = df
|
||||
let df2 = -ex * dfx - ey * dfy
|
||||
|
||||
let q1 = 4.0 * ex * ex / L3
|
||||
let q2 = 4.0 * ex * ey / L3
|
||||
let q3 = 4.0 * ey * ey / L3
|
||||
|
||||
let s1 = (6.0 * (f0 - f1) - 2.0 * df2) * ex / L3
|
||||
let s2 = (6.0 * (f0 - f1) - 2.0 * df2) * ey / L3
|
||||
|
||||
let Q = [| q1; q2; q2; q3 |]
|
||||
let s = [| s1; s2 |]
|
||||
Q, s
|
||||
|
||||
// Set up linear system to solve (Qx = s), accumulated from N neighboring points
|
||||
let private setupQandS (p0: Pos) (f0: float) (p_n: Pos []) (f_n: float []) (df: float * float) =
|
||||
let calcQS_n = calcContributionQandS df p0 f0
|
||||
|
||||
let (<+>) a b = Array.map2 (+) a b
|
||||
let rec accumulateQandS Q s idx =
|
||||
match idx with
|
||||
| [] -> Q, s
|
||||
| x :: xs ->
|
||||
let Q', s' = calcQS_n p_n[x] f_n[x]
|
||||
accumulateQandS (Q <+> Q') (s <+> s') xs
|
||||
|
||||
let Q0 = [| 0.0; 0.0; 0.0; 0.0 |]
|
||||
let s0 = [| 0.0; 0.0 |]
|
||||
let idx = [ 0 .. (p_n.Length - 1) ]
|
||||
accumulateQandS Q0 s0 idx
|
||||
|
||||
// Estimates gradient in point P_0 based on neighboring function values
|
||||
// p0: Point of interest
|
||||
// f0: Function value in POI
|
||||
// p_n: Neighboring points
|
||||
// f_n: Function values in neighboring points
|
||||
// thrs: Convergence threshold in iterative procedure
|
||||
//
|
||||
// Ref: Based on SciPy implementation which uses the following
|
||||
// G. Nielson "A method for interpolating scattered data based upon a minimum norm network" (1983)
|
||||
// R. J. Renka and A. K. Cline "A Triangle-based C1 interpolation method" (1984)
|
||||
let calcGradientFromNeighbors (p0: Pos) (f0: float) (p_n: Pos []) (f_n: float []) (thrs: float) : float * float =
|
||||
let setupQS = setupQandS p0 f0 p_n f_n
|
||||
|
||||
let rec convergeGradient (error: float) (df: float * float) =
|
||||
if error < thrs then
|
||||
df
|
||||
else
|
||||
let Q, s = setupQS df
|
||||
let detQ = Q[0] * Q[3] - Q[1] * Q[2]
|
||||
|
||||
if (abs detQ) < thrs * thrs then
|
||||
df
|
||||
else
|
||||
let r0 = (Q[3] * s[0] - Q[1] * s[1]) / detQ
|
||||
let r1 = (-Q[2] * s[0] + Q[0] * s[1]) / detQ
|
||||
let dfx, dfy = df
|
||||
let change = max (abs (dfx + r0)) (abs (dfy + r1))
|
||||
convergeGradient change (-r0, -r1)
|
||||
|
||||
let err0 = 1.0
|
||||
let grad0 = 0.0, 0.0
|
||||
convergeGradient err0 grad0
|
||||
|
||||
// Given 1D centroid evaluator "evalCentroid", compute gradient in single vertex
|
||||
let calcVertexGradient
|
||||
(grid: ExtendedGrid)
|
||||
(evalCentroid: int -> float)
|
||||
(thrs: float)
|
||||
(idx: NodeIdx, f0: float)
|
||||
: float * float =
|
||||
let neighbors =
|
||||
grid.getElemsSurroundingNode idx
|
||||
|> Array.filter (fun i -> i <> idx)
|
||||
|
||||
// Get data from vertex point of interest
|
||||
let p0 = getVertex grid idx
|
||||
|
||||
// Get data from neighbor elements of POI
|
||||
let p_n = neighbors |> Array.map (getCentroid grid)
|
||||
let f_n = neighbors |> Array.map evalCentroid
|
||||
|
||||
calcGradientFromNeighbors p0 f0 p_n f_n thrs
|
||||
|
||||
// Given 1D vertex evaluator "evalVertex", compute gradient in single vertex
|
||||
let calcVertexGradient'
|
||||
(grid: ExtendedGrid)
|
||||
(evalVertex: int -> float)
|
||||
(thrs: float)
|
||||
(idx: NodeIdx, f0: float)
|
||||
: float * float =
|
||||
let neighbors =
|
||||
grid.getNodesSurroundingNode idx
|
||||
|> Array.filter (fun i -> i <> idx)
|
||||
|
||||
// Get data from vertex point of interest
|
||||
let p0 = getVertex grid idx
|
||||
|
||||
// Get data from neighbor elements of POI
|
||||
let p_n = neighbors |> Array.map (getVertex grid)
|
||||
let f_n = neighbors |> Array.map evalVertex
|
||||
|
||||
calcGradientFromNeighbors p0 f0 p_n f_n thrs
|
||||
18
src/Grid.fs
18
src/Grid.fs
@@ -20,20 +20,6 @@ type Node = float * float
|
||||
type Pos = float * float
|
||||
type Field = (float * float) []
|
||||
|
||||
type Pos3D = float * float * float
|
||||
|
||||
type Field4D =
|
||||
{
|
||||
Frame0: int
|
||||
Frame1: int
|
||||
Time0: float
|
||||
Time1: float
|
||||
Field0: (single * single * single) [] []
|
||||
Field1: (single * single * single) [] []
|
||||
Siglays: single []
|
||||
Bath: single []
|
||||
}
|
||||
|
||||
type BBox =
|
||||
{
|
||||
minX: float
|
||||
@@ -374,7 +360,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 +687,4 @@ type ExtendedGrid(grid: IGrid) =
|
||||
elementTree <- binarySerializer.UnPickle<IdxTree> pickle |> Some
|
||||
true
|
||||
else
|
||||
false
|
||||
false
|
||||
|
||||
@@ -22,11 +22,6 @@
|
||||
<Compile Include="NorKyst.fs"/>
|
||||
<Compile Include="NorShelf.fs"/>
|
||||
<Compile Include="Smoothing.fs"/>
|
||||
<Compile Include="Gradient.fs"/>
|
||||
<Compile Include="CloughTocher.fs"/>
|
||||
<Compile Include="EvaluateVertex.fs"/>
|
||||
<Compile Include="Evaluate2D.fs"/>
|
||||
<Compile Include="Evaluate4D.fs"/>
|
||||
</ItemGroup>
|
||||
<ItemGroup>
|
||||
<PackageReference Update="FSharp.Core" Version="6.0.5"/>
|
||||
|
||||
Reference in New Issue
Block a user