oceanbox/Oceanbox.FvcomKit
6
0
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 80 Wiki Activity

feat: Add ElemsAroundElem to NeighborIndex

Browse Source 
Also format with Fantomas
This commit is contained in:
Moritz Jörg
2025-11-05 11:39:57 +01:00
parent bec03ed5ec
commit 71e861e417
4 changed files with 185 additions and 168 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
.editorconfig
View File

@@ -11,6 +11,10 @@ insert_final_newline = false
indent_size = 2
max_line_length= 80
[*.nix]
indent_size = 2
max_line_length= 80
[*.fs]
max_line_length= 120
@@ -21,6 +25,7 @@ fsharp_space_before_uppercase_invocation = true
fsharp_blank_lines_around_nested_multiline_expressions = false
fsharp_newline_between_type_definition_and_members = false
fsharp_multiline_bracket_style = stroustrup
fsharp_multi_line_lambda_closing_newline = true
fsharp_array_or_list_multiline_formatter = character_width
fsharp_max_array_or_list_width = 70

16
shell.nix
View File

@@ -2,15 +2,19 @@ with import <nixpkgs> { };
let
dotnet-sdk = dotnet-sdk_9;
in
mkShell {
buildInputs = [
netcdf
mkShell rec {
packages = [
bun
dotnet-sdk
fsautocomplete
fantomas
];
buildInputs = [
netcdf
stdenv.cc.cc.lib
];
DOTNET_ROOT = "${dotnet-sdk}/share/dotnet";
shellHook = ''
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${stdenv.cc.cc.lib}/lib;
'';
LD_LIBRARY_PATH = pkgs.lib.makeLibraryPath buildInputs;
}

233
src/Grid.fs
View File

@@ -9,7 +9,7 @@ open ProjNet.FSharp
open MessagePack
open MBrace.FsPickler
//open FsKDTree
open KdTree // C# version
open KdTree // NOTE: C# version
open Types
@@ -22,7 +22,7 @@ type Node = float * float
type Pos = float * float
type Leaf<'a> = { Pos: Pos; Data: 'a }
type Field = (float * float) []
type Field = (float * float) array
type Cell = NodeIdx * NodeIdx * NodeIdx
@@ -30,16 +30,15 @@ type IGrid =
abstract getVertex: int -> Vertex
abstract getCell: int -> Cell
abstract getCellVertices: int -> Vertex * Vertex * Vertex
abstract getVertices: unit -> Vertex []
abstract getCells: unit -> Cell []
abstract getVertices: unit -> Vertex array
abstract getCells: unit -> Cell array
abstract getBoundingBox: unit -> BBox
type Grid =
{
type Grid = {
Elem: Elem array
Nodes: Node array
BBox: BBox
}
} with
interface IGrid with
member this.getVertex n = this.Nodes[n]
member this.getCell n = this.Elem[n]
@@ -49,29 +48,25 @@ type Grid =
member this.getVertices() = this.Nodes
member this.getCells() = this.Elem
member this.getBoundingBox() = this.BBox
static member empty =
{
Elem = Array.empty
Nodes = Array.empty
BBox = BBox.empty
}
static member empty = { Elem = Array.empty; Nodes = Array.empty; BBox = BBox.empty }
type ElemsAroundNode = Map<NodeIdx, ElemIdx []>
type NodesAroundNode = Map<NodeIdx, NodeIdx []>
type ElemsAroundNode = Map<NodeIdx, ElemIdx array>
type NodesAroundNode = Map<NodeIdx, NodeIdx array>
type ElemsAroundElem = Map<ElemIdx, ElemIdx array>
type NeighborIndex =
{
type NeighborIndex = {
ElemsAroundNode: ElemsAroundNode
NodesAroundNode: NodesAroundNode
}
static member empty = { ElemsAroundNode = Map.empty; NodesAroundNode = Map.empty }
ElemsAroundElem: ElemsAroundElem
} with
static member empty = { ElemsAroundNode = Map.empty; NodesAroundNode = Map.empty; ElemsAroundElem = Map.empty }
type private Ean = Map<NodeIdx, ElemIdx list>
// NOTE(SimenLK): The amount of items to be stored in the trees leafs
// let treeLeafSize = LeafNodeSize 64
let private createTree (points: Leaf<int> []) =
let private createTree (points: Leaf<int> array) =
let tree = KdTree<float, int> (2, KdTree.Math.DoubleMath ())
points
|> Array.iter (fun a -> tree.Add ([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
@@ -89,12 +84,9 @@ let private makeElemsSurroundingNodeMap (elem: Elem array) : ElemsAroundNode =
elem
|> Array.fold
(fun (n, acc) (a, b, c) ->
let acc' =
acc
|> addElIdx a n
|> addElIdx b n
|> addElIdx c n
n + 1, acc')
let acc' = acc |> addElIdx a n |> addElIdx b n |> addElIdx c n
n + 1, acc'
)
(0, Map.empty)
|> snd
|> Map.mapValues toArray
@@ -108,12 +100,9 @@ let private makeElemsSurroundingNodeMap' (elem: Elem array) : ElemsAroundNode =
elem
|> Array.fold
(fun (n, acc) (a, b, c) ->
let acc' =
acc
|> addElemIdx a n
|> addElemIdx b n
|> addElemIdx c n
n + 1, acc')
let acc' = acc |> addElemIdx a n |> addElemIdx b n |> addElemIdx c n
n + 1, acc'
)
(0, Map.empty)
|> snd
|> Map.mapValues toArray
@@ -124,13 +113,26 @@ let private makeNodesSurroudingNodeMap (n2e: ElemsAroundNode) (elem: Elem array)
n
|> Array.collect (fun x ->
let n1, n2, n3 = elem[x]
[| n1; n2; n3 |])
|> Array.distinct)
[| n1; n2; n3 |]
)
|> Array.distinct
)
let getSurrounding (idx: Map<int, int []>) (a, b, c) =
[| idx[a]; idx[b]; idx[c] |]
let private makeElemsSurroundingElemMap (ean: ElemsAroundNode) (elem: Elem array) : ElemsAroundElem =
elem
|> Array.mapi (fun elemIdx (n1, n2, n3) ->
// For each element, find all elements that share any of its nodes
let surroundingElems =
[| ean[n1]; ean[n2]; ean[n3] |]
|> Array.concat
|> Array.distinct
|> Array.filter (fun x -> x <> elemIdx) // Remove self
elemIdx, surroundingElems
)
|> Map.ofArray
let getSurrounding (idx: Map<int, int[]>) (a, b, c) =
[| idx[a]; idx[b]; idx[c] |] |> Array.concat |> Array.distinct
let makeNeighborIndex (grid: IGrid) =
let elem = grid.getCells ()
@@ -138,14 +140,17 @@ let makeNeighborIndex (grid: IGrid) =
{
ElemsAroundNode = ean
NodesAroundNode = makeNodesSurroudingNodeMap ean elem
ElemsAroundElem = makeElemsSurroundingElemMap ean elem
}
let getElemsSurroundingNode (idx: NeighborIndex) n = idx.ElemsAroundNode[n]
let getNodesSurroundingNode (idx: NeighborIndex) n = idx.NodesAroundNode[n]
let getNodesSurroundingElem (idx: NeighborIndex) (grid: Grid) e =
getSurrounding idx.NodesAroundNode grid.Elem[e]
let getNodesSurroundingElem (idx: NeighborIndex) e = idx.ElemsAroundElem[e]
// let getNodesSurroundingElem (idx: NeighborIndex) (grid: Grid) e =
// getSurrounding idx.NodesAroundNode grid.Elem[e]
let getElemsSurroundingElem (idx: NeighborIndex) (grid: Grid) e =
getSurrounding idx.ElemsAroundNode grid.Elem[e]
@@ -182,10 +187,16 @@ let bboxToLngLat (coordsys: CoordinateSystem) b =
let projectBBox proj b =
let x0, y0 = proj (b.minX, b.minY)
let x1, y1 = proj (b.maxX, b.maxY)
{ minX = x0; maxX = x1; minY = y0; maxY = y1; center = proj b.center }
{
minX = x0
maxX = x1
minY = y0
maxY = y1
center = proj b.center
}
let projectGrid proj (grid: Grid) : Grid =
{ grid with
let projectGrid proj (grid: Grid) : Grid = {
grid with
Nodes = grid.Nodes |> Array.Parallel.map proj
BBox = projectBBox proj grid.BBox
}
@@ -193,18 +204,12 @@ let projectGrid proj (grid: Grid) : Grid =
let rescaleGrid (factor: float) (grid: Grid) : Grid =
let nodes' = grid.Nodes |> Array.Parallel.map (fun (x, y) -> factor * x, factor * y)
let bbox' = calcBBox nodes'
{ grid with
Nodes = nodes'
BBox = bbox'
}
{ grid with Nodes = nodes'; BBox = bbox' }
let translateGrid (x0, y0) grid =
let nodes' = grid.Nodes |> Array.Parallel.map (fun (x, y) -> x + x0, y + y0)
let bbox' = calcBBox nodes'
{ grid with
Nodes = nodes'
BBox = bbox'
}
{ grid with Nodes = nodes'; BBox = bbox' }
let toWebMercator (coordsys: CoordinateSystem) (grid: Grid) =
let toWebMercator = makeTransform coordsys CoordSys.EPSG3857
@@ -231,20 +236,19 @@ let private readGrdHeader (h: string []) =
let nodes = parse h[0]
let nele = parse h[1]
Some (nodes, nele)
with
| _ -> None
with _ ->
None
let private readObcHeader (h: string) =
try
(chomp h)[4] |> int |> Some
with
| _ -> None
with _ ->
None
let private reader (parser: string[] -> 'a) (f: string array) =
try
f |> Array.map (chomp >> parser) |> Some
with
| e ->
with e ->
Log.Error e.Message
None
@@ -282,13 +286,13 @@ let readGrdFile (filename: string) =
let elem = readElem els
let nodes = readNodes nds
let toGrid e n = { Elem = e; Nodes = n; BBox = calcBBox n }
toGrid <!> elem <*> nodes)
toGrid <!> elem <*> nodes
)
let readObcFile (filename: string) =
let f = System.IO.File.ReadAllLines filename
let hdr, rest = Array.splitAt 1 f
readObcHeader hdr[0]
|> Option.bind (fun _ -> readObc rest)
readObcHeader hdr[0] |> Option.bind (fun _ -> readObc rest)
module Boundary =
let normalizeElement (a, b, c) =
@@ -303,10 +307,7 @@ module Boundary =
| Some v -> Map.add edge (n :: v) a
| None -> Map.add edge [ n ] a
let appendEdges a (n, x: Edge[]) =
a
|> appendEdge (n, x[0])
|> appendEdge (n, x[1])
|> appendEdge (n, x[2])
a |> appendEdge (n, x[0]) |> appendEdge (n, x[1]) |> appendEdge (n, x[2])
grid.Elem
|> Array.mapi normElIdx
|> Array.fold appendEdges Map.empty
@@ -324,8 +325,7 @@ module Boundary =
|> Map.mapValues (fun x -> List.head x, List.last x)
let makeBoundaryByElementMap (edgeMap: Map<Edge, ElemIdx>) =
edgeMap
|> Map.fold (fun a k v -> Map.add v k a) Map.empty
edgeMap |> Map.fold (fun a k v -> Map.add v k a) Map.empty
let getBoundaryNodesArray (edgeMap: Map<Edge, ElemIdx>) =
edgeMap
@@ -345,8 +345,7 @@ module Util =
(x0 + x1 + x2) / 3f, (y0 + y1 + y2) / 3f
static member calcArea((x0, y0), (x1, y1), (x2, y2)) =
x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
|> (*) 0.5
x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1) |> (*) 0.5
// static member calcArea((x0, y0), (x1, y1), (x2, y2)) =
// x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
@@ -364,7 +363,8 @@ 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
@@ -373,9 +373,8 @@ module Util =
[| 0 .. nIdx.NodesAroundNode.Count - 1 |] // total number of nodes
|> Array.Parallel.map (fun i ->
let ns = nIdx.ElemsAroundNode[i]
ns
|> Array.fold (fun a n -> s[n] + a) 0.
|> (*) (1. / float ns.Length))
ns |> Array.fold (fun a n -> s[n] + a) 0. |> (*) (1. / float ns.Length)
)
static member speedToNodal (nIdx: NeighborIndex) (s: (float * float)[]) =
[| 0 .. nIdx.NodesAroundNode.Count - 1 |] // total number of nodes
@@ -385,9 +384,11 @@ module Util =
|> Array.fold
(fun a n ->
let u, v = s[n]
sqrt (u * u + v * v) + a)
sqrt (u * u + v * v) + a
)
0.
|> (*) (1. / float ns.Length))
|> (*) (1. / float ns.Length)
)
static member velocityToNodal (nIdx: NeighborIndex) (s: (float * float)[]) =
[| 0 .. nIdx.NodesAroundNode.Count - 1 |] // total number of nodes
@@ -398,9 +399,11 @@ module Util =
|> Array.fold
(fun (au, av) n ->
let u, v = s[n]
u + au, v + av)
u + au, v + av
)
(0., 0.)
|> fun (u, v) -> u / n, v / n)
|> fun (u, v) -> u / n, v / n
)
type Node =
static member calcNodeControlArea (idx: NeighborIndex) (grid: IGrid) =
@@ -419,8 +422,10 @@ module Util =
let centroid = Element.calcCentroid (p0, p1, p2)
let a1 = Element.calcArea (p0', p1', p2')
let a2 = Element.calcArea (p1', centroid, p2')
a1 + a2)
|> Array.sum)
a1 + a2
)
|> Array.sum
)
static member calcNodeArea (idx: NeighborIndex) (grid: IGrid) =
let nodes = grid.getVertices ()
@@ -428,7 +433,8 @@ module Util =
|> Array.Parallel.map (fun n ->
getElemsSurroundingNode idx n
|> Array.map (grid.getCellVertices >> Element.calcArea)
|> Array.sum)
|> Array.sum
)
let calcCentroids (grid: IGrid) =
let n = grid.getVertices ()
@@ -437,12 +443,12 @@ module Util =
let p0 = n[a]
let p1 = n[b]
let p2 = n[c]
Element.calcCentroid (p0, p1, p2))
Element.calcCentroid (p0, p1, p2)
)
let inline isInsideTriangle (x, y, z) p =
let sign (p1x, p1y) (p2x, p2y) (p3x, p3y) =
(p1x - p3x) * (p2y - p3y)
- (p2x - p3x) * (p1y - p3y)
(p1x - p3x) * (p2y - p3y) - (p2x - p3x) * (p1y - p3y)
let d1 = sign p x y
let d2 = sign p y z
@@ -458,7 +464,8 @@ module Util =
grid.getVertices ()
|> Array.mapi (fun i v ->
let x, y = v
{ Pos = x, y; Data = i })
{ Pos = x, y; Data = i }
)
// |> create2DTree treeLeafSize
|> createTree
@@ -466,12 +473,10 @@ module Util =
let buildNearestElementTree (grid: IGrid) =
grid.getCells ()
|> Array.mapi (fun i _ ->
let pos =
i
|> grid.getCellVertices
|> Element.calcCentroid
let pos = i |> grid.getCellVertices |> Element.calcCentroid
// |> fun (x, y) -> { X = x; Y = y }
{ Pos = pos; Data = i })
{ Pos = pos; Data = i }
)
// |> create2DTree treeLeafSize
|> createTree
@@ -494,7 +499,8 @@ module Util =
if isInsideTriangle vx (p0, p1) then
Some leaf.Value
else
None)
None
)
// type private IdxTree = Tree<Leaf<single, ElemIdx> array, Node<single>>
// type private NodeIdxTree = Tree<Leaf<single, NodeIdx> array, Node<single>>
@@ -503,16 +509,14 @@ type private NodeIdxTree = KdTree<float, int>
[<MessagePackObject>]
type BinGrid = {
[<Key(0)>] hash : byte[]
[<Key(1)>] vertices : (float * float)[]
[<Key(2)>] cells : (int * int * int)[]
[<Key(0)>]
hash: byte[]
[<Key(1)>]
vertices: (float * float)[]
[<Key(2)>]
cells: (int * int * int)[]
} with
member this.toGrid (): Grid =
{
Nodes = this.vertices
Elem = this.cells
BBox = calcBBox this.vertices
}
member this.toGrid() : Grid = { Nodes = this.vertices; Elem = this.cells; BBox = calcBBox this.vertices }
type ExtendedGrid(grid: IGrid) =
let mutable nodeTree: NodeIdxTree option = None
@@ -542,8 +546,7 @@ type ExtendedGrid(grid: IGrid) =
member this.Grid = grid
member this.ToGrid() =
{
member this.ToGrid() = {
Nodes = this.Grid.getVertices ()
Elem = this.Grid.getCells ()
BBox = this.Grid.getBoundingBox ()
@@ -580,11 +583,15 @@ type ExtendedGrid(grid: IGrid) =
this.nearestNode p
|> Option.bind (fun n ->
this.getElemsSurroundingNode n
|> Array.fold (fun (a: int option) e ->
|> Array.fold
(fun (a: int option) e ->
if a.IsNone then
let vx = this.Grid.getCellVertices e
if Util.isInsideTriangle vx p then Some n else None
else a) None
else
a
)
None
)
member private this.tryFindElementTwice (grid: IGrid) (tree: IdxTree) ((p0, p1): float * float as p) =
@@ -596,16 +603,18 @@ type ExtendedGrid(grid: IGrid) =
this.getElemsSurroundingElem leaf.Value
|> Array.tryFind (fun eIdx ->
let vx = this.Grid.getCellVertices eIdx
Util.isInsideTriangle vx p))
Util.isInsideTriangle vx p
)
)
)
member this.tryGetElement p =
elementTree
|> Option.bind (fun tree ->
this.tryFindElementTwice grid tree p)
|> Option.bind (fun tree -> this.tryFindElementTwice grid tree p)
|> Option.orElseWith (fun () ->
this.initElementTree ()
this.tryFindElementTwice grid elementTree.Value p)
this.tryFindElementTwice grid elementTree.Value p
)
member this.tryGetElementSloppy(p0, p1 as p) =
match elementTree with
@@ -647,12 +656,10 @@ type ExtendedGrid(grid: IGrid) =
Util.Element.calcCircumscribedCircle triangle
member this.getElemsSurroundingNode n =
getNeighborIdx ()
|> fun idx -> idx.ElemsAroundNode[n]
getNeighborIdx () |> fun idx -> idx.ElemsAroundNode[n]
member this.getNodesSurroundingNode n =
getNeighborIdx ()
|> fun idx -> idx.NodesAroundNode[n]
getNeighborIdx () |> fun idx -> idx.NodesAroundNode[n]
member this.getNodesSurroundingElem e =
let idx = getNeighborIdx ()
@@ -674,9 +681,7 @@ type ExtendedGrid(grid: IGrid) =
let binarySerializer = FsPickler.CreateBinarySerializer ()
if IO.File.Exists fname then
let pickle = IO.File.ReadAllBytes fname
neighborIndex <-
binarySerializer.UnPickle<NeighborIndex> pickle
|> Some
neighborIndex <- binarySerializer.UnPickle<NeighborIndex> pickle |> Some
true
else
false
@@ -687,7 +692,8 @@ type ExtendedGrid(grid: IGrid) =
nodeTree
|> Option.defaultWith (fun () ->
this.initNodeTree ()
nodeTree.Value)
nodeTree.Value
)
let pickle = binarySerializer.Pickle tree
IO.File.WriteAllBytes (fname, pickle)
@@ -707,7 +713,8 @@ type ExtendedGrid(grid: IGrid) =
elementTree
|> Option.defaultWith (fun () ->
this.initElementTree ()
elementTree.Value)
elementTree.Value
)
let pickle = binarySerializer.Pickle tree
IO.File.WriteAllBytes (fname, pickle)

1
src/Oceanbox.FvcomKit.fsproj
View File

@@ -5,6 +5,7 @@
<TargetFramework>net9.0</TargetFramework>
<IsPackable>true</IsPackable>
<PackageId>Oceanbox.FvcomKit</PackageId>
<RuntimeIdentifier>linux-x64</RuntimeIdentifier>
<Authors/>
<Company/>
<Version>5.12.2</Version>