Merge pull request #30 from ZenanH/main

🎉 Release new version v0.1.23

Author: ZenanH

CondaPkg.toml

@@ -7,7 +7,7 @@ numpy = "==2.2.4"
 rtree = "==1.4.0"
 python-gmsh = "==4.13.1"
 trimesh = "==4.6.6"
-shapely = "==2.0.7"
+shapely = "==2.1.0"
 
 [pip.deps]
 rasterio = "==1.4.1"

Project.toml

@@ -1,13 +1,12 @@
 name = "MaterialPointGenerator"
 uuid = "18deba1c-5063-41ee-a388-01f264c5a914"
 authors = ["ZenanH <zenan.huo@outlook.com>"]
-version = "0.1.22"
+version = "0.1.23"
 
 [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"
+LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
@@ -19,8 +18,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 [compat]
 CondaPkg = "0.2"
 DelimitedFiles = "1"
-Gmsh = "0.3"
-GMT = "1.28"
+LibGEOS = "0.9"
 LiveServer = "1.5"
 NearestNeighbors = "0.4"
 PrecompileTools = "1.2"

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.22-pink)]()
+[![Version](https://img.shields.io/badge/version-v0.1.23-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. 

docs/src/assets/showcase/slbl_gui.png

@@ -14,12 +14,12 @@ The DEM file is a simple three-column array. Each DEM must be rasterized to ensu
 
 ```@docs
 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]
+    k         ::T1 = 10, 
+    p         ::T1 = 2, 
+    trimbounds::AbstractMatrix{T2} = [0.0 0.0], 
+    dembounds ::AbstractVector{T2} = [0.0, 0.0]
 ) where {T1, T2}
 ```
 
@@ -31,9 +31,9 @@ Suppose we have a DEM and we want to close it with a base plane, for example, at
 
 ```@docs
 dem2particle(
-    dem   ::Matrix{T2}, 
-    h     ::T2, 
-    bottom::T2
+    dem   ::AbstractMatrix{T2}, 
+    h     ::Real, 
+    bottom::Real
 ) where T2
 ```
 
@@ -43,9 +43,9 @@ If the base used to close DEM-1 is not a flat surface, we can designate another
 
 ```@docs
 dem2particle(
-    dem   ::Matrix{T2}, 
-    h     ::T2, 
-    bottom::Matrix{T2}
+    dem   ::AbstractMatrix{T2}, 
+    h     ::Real, 
+    bottom::AbstractMatrix{T2}
 ) where T2
 ```
 
@@ -63,10 +63,10 @@ To prepare for this, you need to have layered surface files in DEM format, which
 
 ```@docs
 dem2particle(
-    dem   ::Matrix{T2}, 
-    h     ::T2, 
-    bottom::T2,
-    layer ::Vector{Matrix{T2}}
+    dem   ::AbstractMatrix{T2}, 
+    h     ::Real, 
+    bottom::Real,
+    layer ::AbstractVector{AbstractMatrix{T2}}
 ) where T2
 ```
 
@@ -76,9 +76,9 @@ This workflow also supports the case where a bottom DEM is provided:
 
 ```@docs
 dem2particle(
-    dem   ::Matrix{T2}, 
-    h     ::T2, 
-    bottom::Matrix{T2},
+    dem   ::AbstractMatrix{T2}, 
+    h     ::Real, 
+    bottom::AbstractMatrix{T2},
     layer ::Vector{Matrix{T2}}
 ) where T2
 ```

docs/src/workflow/DEM.md

@@ -11,7 +11,7 @@ For the second scenario, to support complex polygons in any situation, we need t
 We are preparing to discretize a two-dimensional pentagon, noting that the specified distance is between the material points in both the x and y directions.
 
 ```@docs
-polygon2particle(polygon, lpx, lpy)
+polygon2particle(polygon::AbstractMatrix{T}, lpx, lpy) where T<:Real
 ```
 
 ```julia
@@ -39,7 +39,7 @@ julia> pts = polygon2particle(polygon, 0.01, 0.01)
 Please use pre-processing tools like Gmsh or MeshLab to ensure that the current mesh is closed in advance.
 
 ```@docs
-polygon2particle(stl_file::String, output_file::String, h; verbose::Bool=true)
+polygon2particle(stl_file::String, output_file::String, h; verbose::Bool=false)
 ```
 
 ## Advanced
@@ -49,7 +49,11 @@ This involves some advanced operations for partitioning the generated material p
 For the first case, we can utilize a practical function to check if a point is inside the polygon. 
 
 ```@docs
-particle_in_polygon(px::T, py::T, polygon) where T
+particle_in_polygon(
+    polygon::AbstractMatrix{T}, 
+    px     ::Real, 
+    py     ::Real
+) where T<:Real
 ```
 
 For the second case, in addition to the `.stl` file, we also need to provide a `.msh` file.
@@ -61,6 +65,6 @@ polygon2particle(
     output_file::String,
     nid_file   ::String,
     h; 
-    verbose    ::Bool=true
+    verbose    ::Bool=false
 )
 ```

docs/src/workflow/polygon.md

@@ -11,21 +11,20 @@
 
 module MaterialPointGenerator
 
-using CondaPkg, DelimitedFiles, Gmsh, NearestNeighbors, Printf, PythonCall
+using CondaPkg, DelimitedFiles, LibGEOS, 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 rtree        = PythonCall.pynew()
+const shapely      = 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()
@@ -36,29 +35,26 @@ const pygmsh       = PythonCall.pynew()
 const pyKDTree     = PythonCall.pynew()
 const pytime       = PythonCall.pynew()
 const embreex      = PythonCall.pynew()
-const ConvexHull   = PythonCall.pynew()
-const v_contains   = PythonCall.pynew()
 
 function __init__()
     @info "checking environment..."
     # import Python modules
     PythonCall.pycopy!(trimesh , pyimport("trimesh" ))
     PythonCall.pycopy!(np      , pyimport("numpy"   ))
+    PythonCall.pycopy!(rtree   , pyimport("rtree"   ))
+    PythonCall.pycopy!(shapely , pyimport("shapely" ))
     PythonCall.pycopy!(rasterio, pyimport("rasterio"))
     PythonCall.pycopy!(pygmsh  , pyimport("gmsh"    ))
     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)
-    PythonCall.pycopy!(ConvexHull  , pyimport("scipy.spatial"     ).ConvexHull  )
-    PythonCall.pycopy!(v_contains  , pyimport("shapely.vectorized").contains    )
     if !Sys.isapple()
         try 
             if !haskey(CondaPkg.current_pip_packages(), "embreex")
@@ -94,6 +90,8 @@ include(joinpath(@__DIR__, "utils.jl"        ))
     @compile_workload begin
         # functions in utils.jl
         getnormals(points)
+        # function in dem.jl
+        getpolygon(rand(10, 2))
     end
 end
 

src/MaterialPointGenerator.jl

@@ -38,19 +38,20 @@ end
     areaPBC = cross2D(v2x - px, v2y - py, v3x - px , v3y - py ) / areaABC
     areaAPC = cross2D(px - v1x, py - v1y, v3x - v1x, v3y - v1y) / areaABC
 
-    alpha = areaPBC
-    beta  = areaAPC
-    gamma = 1 - alpha - beta
+    α = areaPBC
+    β = areaAPC
+    γ = 1 - α - β
 
     # 3) 用 (alpha,beta,gamma) 插值三顶点的 z 值
-    z = alpha*v1z + beta*v2z + gamma*v3z
+    z = α*v1z + β*v2z + γ*v3z
 
-    return floor((z - mzmin) / h) * h + mzmin
+    return z#, floor((z - mzmin) / h) * h + mzmin
 end
 
-function check_projection!(meshdata::GmshMesh{T1, T2}, h::T2, p2t::Vector, mxmin::T2, 
+function check_projection!(meshdata::DataMesh{T1, T2}, h::T2, p2t::Vector, mxmin::T2, 
     mymin::T2, niy::T1) where {T1, T2}
 
+    p2t = copy(p2t)
     @inbounds for i in axes(meshdata.data, 1)
         v1x, v1y = meshdata.data[i][1], meshdata.data[i][2]
         v2x, v2y = meshdata.data[i][4], meshdata.data[i][5] 
@@ -77,7 +78,7 @@ function check_projection!(meshdata::GmshMesh{T1, T2}, h::T2, p2t::Vector, mxmin
     return nothing
 end
 
-function fill_voxel!(pts::Vector{Bool}, p2t::Vector, meshdata::GmshMesh{T1, T2}, niy::T1, 
+function fill_voxel!(pts::Vector{Bool}, p2t::Vector, meshdata::DataMesh{T1, T2}, niy::T1, 
     mxmin::T2, mymin::T2, mzmin::T2, mxmax::T2, mymax::T2, mzmax::T2, h::T2) where {T1, T2}
 
     @inbounds for i in axes(p2t, 1)
@@ -93,8 +94,8 @@ function fill_voxel!(pts::Vector{Bool}, p2t::Vector, meshdata::GmshMesh{T1, T2},
                 push!(tri_order, z)
             end
             sort!(unique!(tri_order))
-            tri_order[  1] = ceil( T1, (tri_order[  1] - mzmin) / h) * h + mzmin
-            tri_order[end] = floor(T1, (tri_order[end] - mzmin) / h) * h + mzmin
+            #tri_order[  1] = ceil( T1, (tri_order[  1] - mzmin) / h) * h + mzmin
+            #tri_order[end] = floor(T1, (tri_order[end] - mzmin) / h) * h + mzmin
             fill_layers!(pts, px, py, tri_order, mxmin, mymin, mzmin, mxmax, mymax, mzmax, 
                 h)
         end

src/_polyhedron/_utils.jl

@@ -2,51 +2,26 @@
     readmesh(mesh_file::String; precision="FP32")
 Description:
 ---
-Read the mesh file (`.stl`) generated by Gmsh. The function will return the vertices and 
-faces for triangle. It's necessary to check the mesh qulity before using this function. 
-You can do it in MeshLab or Gmsh GUI. `precision` can be "FP64" or "FP32" in `String`.
+Read the mesh file (`.stl`). The function will return the vertices and faces for triangle. 
+It's necessary to check the mesh qulity before using this function. You can do it in MeshLab 
+or Gmsh GUI. `precision` can be "FP64" or "FP32" in `String`.
 """
-function readmesh(mesh_file; precision="FP32")::GmshMesh
+function readmesh(mesh_file; precision="FP32")::DataMesh
     T1 = precision == "FP32" ? Int32 : Int64
     T2 = precision == "FP32" ? Float32 : Float64
 
-    # Initialize gmsh
-    gmsh.initialize()
-    gmsh.option.setNumber("General.Terminal", 0)
-    gmsh.option.setNumber("Mesh.Optimize", 1)
-    gmsh.open(mesh_file)
+    mesh = trimesh.load(mesh_file)
+    pyconvert(Bool, mesh.is_watertight) || error(
+        "The mesh is not watertight. Please check the mesh quality.")
 
-    # Optimize mesh to remove duplicate nodes
-    gmsh.model.mesh.optimize("Netgen")
+    coords = pyconvert(Array{T2}, mesh.vertices)
+    triang = pyconvert(Array{T1}, mesh.faces) .+ 1
 
-    # Get mesh data
-    node_tags, node_coords, _ = gmsh.model.mesh.getNodes()
-    num_nodes = length(node_tags)
-    vertices = reshape(node_coords, 3, num_nodes)
-
-    # check if the mesh is watertight
-    element_types, element_tags, element_node_tags = gmsh.model.mesh.getElements()
-    entities = gmsh.model.getEntities(2)  # get all 2d entities
-    boundary = gmsh.model.getBoundary(entities, true, true, true)
-    isempty(boundary) || error("3D object mesh is not closed")
-
-    gmsh.finalize()
-    @info "Mesh data loaded and checking now"
-
-    # Generate vertices and faces in Julia array
-    triangle_type = 2
-    triangle_index = findfirst(==(triangle_type), element_types)
-    isnothing(triangle_index) && error("No triangles found in the mesh")
-
-    triangle_node_tags = element_node_tags[triangle_index]
-    num_triangles = length(triangle_node_tags) ÷ 3
-
-    # Map node tags to consecutive indices
-    node_map = Dict(tag => i for (i, tag) in enumerate(node_tags))
-    faces = reshape(T1.([node_map[tag] for tag in triangle_node_tags]), 3, num_triangles)
-
-    maximum(faces) ≠ size(vertices, 2) && error("Invalid mesh data.")
-    @info "Mesh data checked"
+    v_tmp = Array{T2}(coords)
+    zoffset = minimum(v_tmp[:, 3]) > 10 ? 0 : (10-minimum(v_tmp[:, 3]))
+    v_tmp[:, 3] .+= zoffset
+    vertices = Array{T2}(v_tmp')
+    faces = Array{T1}(triang')
 
     data = [Vector{T2}() for _ in 1:size(faces, 2)]
     @inbounds for i in axes(faces, 2), j in 1:3
@@ -59,18 +34,17 @@ function readmesh(mesh_file; precision="FP32")::GmshMesh
     xmax, ymax, zmax = map(x -> maximum(x), eachrow(vertices))
     bounds = [xmin, ymin, zmin, xmax, ymax, zmax]
 
-    return GmshMesh{T1, T2}(vertices, faces, data, bounds)
+    return DataMesh{T1, T2}(vertices, faces, data, bounds, zoffset)
 end
 
-function _voxelize(meshdata::GmshMesh{T1, T2}, h::Real) where {T1, T2}
+function _voxelize(meshdata::DataMesh{T1, T2}, h::Real) where {T1, T2}
     h > 0 || error("h must be positive")
     h = T2(h)
 
     # get bbox of the entire mesh
-    tmp = 4 .+ rand(T2, 6)
-    mxmin, mxmax = T2(meshdata.bounds[1] - tmp[1]*h), T2(meshdata.bounds[4] + tmp[4]*h)
-    mymin, mymax = T2(meshdata.bounds[2] - tmp[2]*h), T2(meshdata.bounds[5] + tmp[5]*h)
-    mzmin, mzmax = T2(meshdata.bounds[3] - tmp[3]*h), T2(meshdata.bounds[6] + tmp[6]*h)
+    mxmin, mxmax = T2(meshdata.bounds[1])-6, T2(meshdata.bounds[4]+6)
+    mymin, mymax = T2(meshdata.bounds[2])-6, T2(meshdata.bounds[5]+6)
+    mzmin, mzmax = T2(meshdata.bounds[3])-6, T2(meshdata.bounds[6]+6)
 
     # get the number of elements in x, y, z directions
     nex  = ceil(T1, (mxmax - mxmin) / h)
@@ -104,7 +78,7 @@ function _voxelize(meshdata::GmshMesh{T1, T2}, h::Real) where {T1, T2}
 end
 
 function _voxelize(
-    meshdata  ::GmshMesh{T1, T2}, 
+    meshdata  ::DataMesh{T1, T2}, 
     msh_file  ::String, 
     output_xyz::String, 
     output_nid::String,