Merge pull request #29 from ZenanH/main

🎉 Release new version v0.1.22

Author: ZenanH

Project.toml

@@ -1,21 +1,29 @@
 name = "MaterialPointGenerator"
 uuid = "18deba1c-5063-41ee-a388-01f264c5a914"
 authors = ["ZenanH <zenan.huo@outlook.com>"]
-version = "0.1.21"
+version = "0.1.22"
 
 [deps]
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
+GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 CondaPkg = "0.2"
 DelimitedFiles = "1"
 Gmsh = "0.3"
+GMT = "1.28"
+LiveServer = "1.5"
 NearestNeighbors = "0.4"
+PrecompileTools = "1.2"
 PythonCall = "0.9"
 julia = "1.11"
 

README.md

@@ -3,7 +3,7 @@
 [![CI](https://github.com/LandslideSIM/MaterialPointGenerator.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/LandslideSIM/MaterialPointGenerator.jl/actions/workflows/ci.yml) 
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg?logo=quicklook)](https://LandslideSIM.github.io/MaterialPointGenerator.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-red.svg?logo=quicklook)](https://LandslideSIM.github.io/MaterialPointGenerator.jl/dev)
-[![Version](https://img.shields.io/badge/version-v0.1.21-pink)]()
+[![Version](https://img.shields.io/badge/version-v0.1.22-pink)]()
 
 During the EGU2023 conference, when I presented a high-performance MPM  (Material Point Method) solver, I was asked, 
 "How do you discretize the computational model for the MPM?" I didn't have a clear answer (I didn't even consider it a problem) because the models were relatively simple and could be generated directly using some straightforward functions. However, as computational models gradually became more complex and diverse, I began to realize that this was indeed a very good question. The preprocessing for MPM should not be a computationally intensive task; it should be fast enough. Yet, I couldn't find a "plug-and-play" generalized code for this purpose. Some literatures have contributed to this issue, and I built upon their work to create a comprehensive and refined julia package. 
@@ -30,6 +30,7 @@ julia> ]
 - [x] Particle generation from a Digital Elevation Model (DEM) file  
 - [x] Automatically interpolate DEM files with support for shape trimming
 - [x] Attach attributions to the particles
+- [x] SLBL and boundary selector interface
 
 ## Showcases 🎲
 
@@ -45,6 +46,10 @@ julia> ]
 |:--------:|:---:|
 | <img src="docs/src/example/image9.png" width="300"> | <img src="docs/src/example/image10.png" width="360"> |
 
+| SLBL |
+|:----:|
+| <img src="docs/src/assets/showcase/slbl_gui.png" width="660"> |
+
 ## Acknowledgement 👍
 
 This project is sponserd by [Risk Group | Université de Lausanne](https://wp.unil.ch/risk/) and [China Scholarship Council [中国国家留学基金管理委员会]](https://www.csc.edu.cn/).

docs/src/assets/showcase/slbl_gui.png

@@ -11,15 +11,25 @@
 
 module MaterialPointGenerator
 
-using DelimitedFiles, CondaPkg, Gmsh, NearestNeighbors, Printf, PythonCall
+using CondaPkg, DelimitedFiles, Gmsh, NearestNeighbors, Printf, PythonCall
+using LinearAlgebra: mul!, eigen, Symmetric, normalize
+using Statistics: mean
+using GMT: concavehull
+using LiveServer: serve
+using PrecompileTools: @setup_workload, @compile_workload
+
+import LinearAlgebra.LAPACK: syev!
 
 const trimesh      = PythonCall.pynew()
 const np           = PythonCall.pynew()
+const o3d          = PythonCall.pynew()
 const MultiPolygon = PythonCall.pynew()
 const Polygon      = PythonCall.pynew()
+const LineString   = PythonCall.pynew()
 const Point        = PythonCall.pynew()
 const mapping      = PythonCall.pynew()
 const unary_union  = PythonCall.pynew()
+const rotate       = PythonCall.pynew()
 const rasterio     = PythonCall.pynew()
 const rasterize    = PythonCall.pynew()
 const pygmsh       = PythonCall.pynew()
@@ -39,9 +49,11 @@ function __init__()
     PythonCall.pycopy!(pytime  , pyimport("time"    ))
     # import their submudules
     PythonCall.pycopy!(Polygon     , pyimport("shapely.geometry"  ).Polygon     )
+    PythonCall.pycopy!(LineString  , pyimport("shapely.geometry"  ).LineString  )
     PythonCall.pycopy!(Point       , pyimport("shapely.geometry"  ).Point       )
     PythonCall.pycopy!(mapping     , pyimport("shapely.geometry"  ).mapping     )
     PythonCall.pycopy!(unary_union , pyimport("shapely.ops"       ).unary_union )
+    PythonCall.pycopy!(rotate      , pyimport("shapely.affinity"  ).rotate      )
     PythonCall.pycopy!(rasterize   , pyimport("rasterio.features" ).rasterize   )
     PythonCall.pycopy!(pyKDTree    , pyimport("scipy.spatial"     ).cKDTree     )
     PythonCall.pycopy!(MultiPolygon, pyimport("shapely.geometry"  ).MultiPolygon)
@@ -73,6 +85,16 @@ include(joinpath(@__DIR__, "meshgenerator.jl"))
 include(joinpath(@__DIR__, "polygon.jl"      ))
 include(joinpath(@__DIR__, "polyhedron.jl"   ))
 include(joinpath(@__DIR__, "dem.jl"          ))
+include(joinpath(@__DIR__, "slbl/slbl.jl"    ))
+include(joinpath(@__DIR__, "slbl/slbl_gui.jl"))
 include(joinpath(@__DIR__, "utils.jl"        ))
 
+@setup_workload begin
+    points = rand(20, 3)
+    @compile_workload begin
+        # functions in utils.jl
+        getnormals(points)
+    end
+end
+
 end

src/MaterialPointGenerator.jl

@@ -31,9 +31,9 @@ Array of the particles, `tree` is the KDTree of the DEM.
 @views function IDW!(
     k      ::T1, 
     p      ::T1, 
-    dem    ::Matrix{T2}, 
-    idxs   ::Matrix{T1}, 
-    ptslist::Matrix{T2}, 
+    dem    ::AbstractMatrix{T2}, 
+    idxs   ::AbstractMatrix{T1}, 
+    ptslist::AbstractMatrix{T2}, 
     tree   ::KDTree
 ) where {T1, T2}
     @inbounds Threads.@threads for i in axes(ptslist, 1)
@@ -75,12 +75,12 @@ neighbors (10 by default), `p` is the power parameter (2 by default), `trimbound
 boundary of the particles, `dembounds` is the boundary of the DEM.
 """
 @views function rasterizeDEM(
-    dem       ::Matrix{T2},
+    dem       ::AbstractMatrix{T2},
     h         ::T2; 
     k         ::T1        = 10, 
     p         ::T1        = 2, 
-    trimbounds::Matrix{T2}= [0.0 0.0], 
-    dembounds ::Vector{T2}= [0.0, 0.0]
+    trimbounds::AbstractMatrix{T2}= [0.0 0.0], 
+    dembounds ::AbstractVector{T2}= [0.0, 0.0]
 ) where {T1, T2}
     # check input arguments
     size(dem, 2) ≠ 3 && throw(ArgumentError("DEM should have three columns: x, y, z"))
@@ -122,16 +122,17 @@ boundary of the particles, `dembounds` is the boundary of the DEM.
     # move the models to the grid (space = h)
     @. ptslist[:, 3] = ceil(ptslist[:, 3] / h) * h
 
-    return ptslist
+    return sort_pts_xy(ptslist)
 end
 
 @views function rasterizeDEM(
-    dem    ::Matrix{T2},
+    dem    ::AbstractMatrix{T2},
     h      ::T2,
-    polygon::Py; 
-    k      ::T1 = 10, 
-    p      ::T1 = 2
-) where {T1, T2}
+    polygon::AbstractMatrix{T2}; 
+    k      ::Int = 10, 
+    p      ::Int = 2
+) where T2
+    pypolygon = Polygon(np.array(polygon))
     # check input arguments
     size(dem, 2) ≠ 3 && throw(ArgumentError("DEM should have three columns: x, y, z"))
     h > 0 || throw(ArgumentError("h must be positive"))
@@ -144,12 +145,12 @@ end
     ξ0 = meshbuilder(dem_xmin : h : dem_xmax, dem_ymin : h : dem_ymax)
     x_test = np.array(ξ0[:, 1])
     y_test = np.array(ξ0[:, 2])
-    v_in_id = pyconvert(Vector, v_contains(polygon, x_test, y_test))
+    v_in_id = pyconvert(Vector, v_contains(pypolygon, x_test, y_test))
     ptslist = hcat(ξ0[v_in_id, :], zeros(T2, count(v_in_id)))
 
     # create KDTree for DEM
     tree = KDTree(dem[:, 1:2]')
-    idxs = zeros(T1, size(ptslist, 1), k)
+    idxs = zeros(Int, size(ptslist, 1), k)
     IDW!(k, p, dem, idxs, ptslist, tree)
 
     # move the models to the grid (space = h)
@@ -158,19 +159,17 @@ end
     return ptslist
 end
 
-function getpolygon(pts::Matrix)
-    points = size(pts, 2) == 3 ? np.array(pts[:, 1:2]) : np.array(pts)
-
-    @pyexec """
-    def get_polygon(points, ConvexHull, Polygon):
-        hull = ConvexHull(points)
-        hull_points = points[hull.vertices]
-        polygon = Polygon(hull_points)
-        return polygon
-    """ => get_polygon
-
-    polygon = get_polygon(points, ConvexHull, Polygon)
-    return polygon
+function getpolygon(pts::AbstractMatrix; ratio=0.1)
+    pts_col = size(pts, 2)
+    if pts_col == 2
+        points = pts
+    elseif pts_col == 3
+        points = pts[:, 1:2]
+    else
+        throw(ArgumentError("points must be a Nx2 or Nx3 array"))
+    end
+    point_polygon = concavehull(points, ratio, false)
+    return point_polygon.data
 end
 
 """
@@ -183,7 +182,7 @@ Generate particles from a given DEM file. `dem` is a coordinates Array with thre
 size in the MPM simulation. `bottom` is a value, which means the plane `z = bottom`.
 """
 @views function dem2particle(
-    dem   ::Matrix{T2}, 
+    dem   ::AbstractMatrix{T2}, 
     h     ::T2, 
     bottom::T2
 ) where T2
@@ -236,10 +235,10 @@ than the dem. `h` is the space of grid size in `z` direction used in the MPM sim
 surfaces. Note that layers are sorted from top to bottom.
 """
 @views function dem2particle(
-    dem   ::Matrix{T2}, 
+    dem   ::AbstractMatrix{T2}, 
     h     ::T2, 
     bottom::T2,
-    layer ::Vector{Matrix{T2}}
+    layer ::AbstractVector{AbstractMatrix{T2}}
 ) where T2
     # values check
     T1 = T2 == Float32 ? Int32 : Int64
@@ -341,9 +340,9 @@ Array with three columns (x, y, z), which has to be initialized with the struct
 than the dem. `h` is the space of grid size in `z` direction used in the MPM simulation.
 """
 @views function dem2particle(
-    dem   ::Matrix{T2}, 
+    dem   ::AbstractMatrix{T2}, 
     h     ::T2, 
-    bottom::Matrix{T2}
+    bottom::AbstractMatrix{T2}
 ) where T2
     # values check
     T1 = T2 == Float32 ? Int32 : Int64
@@ -401,9 +400,9 @@ Matrix with three columns (x, y, z), which represents the layer surfaces. Note t
 are sorted from top to bottom.
 """
 @views function dem2particle(
-    dem   ::Matrix{T2}, 
+    dem   ::AbstractMatrix{T2}, 
     h     ::T2, 
-    bottom::Matrix{T2},
+    bottom::AbstractMatrix{T2},
     layer ::Vector{Matrix{T2}}
 ) where T2
     # values check