Merge branch 'readin-complex-shapes' into complex-shapes-docs

NEWS.md

@@ -15,12 +15,16 @@ Particle sampling for complex geometries from `.stl` and `.asc` files.
 
 ### Deprecated
 
+## Version 0.2.1
+
+### Highlights
+Particle sampling for complex geometries from `.stl` and `.asc` files.
+
 ## Version 0.2.0
 
 ### Removed
 - Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
 
-
 ## Development Cycle 0.1
 
 ### Highlights

README.md

@@ -22,6 +22,7 @@ Its features include:
   - Methods:  Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM)
 - Fluid-Structure Interaction
 - Preprocess complex shapes
+  - support for STL geometries
 - Output formats:
   - VTK
 

docs/make.jl

@@ -104,6 +104,9 @@ makedocs(sitename="TrixiParticles.jl",
              "Tutorial" => "tutorial.md",
              "Examples" => "examples.md",
              "Visualization" => "visualization.md",
+             "preprocessing" => [
+                 "Sample Geometries" => joinpath("preprocessing", "preprocessing.md"),
+             ],
              "Components" => [
                  "Overview" => "overview.md",
                  "General" => [
@@ -116,9 +119,6 @@ makedocs(sitename="TrixiParticles.jl",
                      "Neighborhood Search" => joinpath("general", "neighborhood_search.md"),
                      "Util" => joinpath("general", "util.md"),
                  ],
-                 "preprocessing" => [
-                     "Sample Geometries" => joinpath("preprocessing", "preprocessing.md"),
-                 ],
                  "Systems" => [
                      "Discrete Element Method (Solid)" => joinpath("systems",
                                                                    "dem.md"),

examples/preprocessing/complex_shape_2d.jl

@@ -1,16 +1,17 @@
 using TrixiParticles
+using Plots
 
 particle_spacing = 0.05
 
-filename = "hexagon"
+filename = "inverted_open_curve"
 file = joinpath("examples", "preprocessing", "data", filename * ".asc")
 
 geometry = load_geometry(file)
 
 trixi2vtk(geometry)
 
 point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
-                                                    #winding_number_factor=0.4,
+                                                    winding_number_factor=0.4,
                                                     hierarchical_winding=true)
 
 # Returns `InitialCondition`
@@ -19,4 +20,9 @@ shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
                              point_in_geometry_algorithm)
 
 trixi2vtk(shape_sampled.initial_condition)
-trixi2vtk(shape_sampled.grid, w=shape_sampled.winding_numbers)
+
+# Plot the winding number field
+plot(InitialCondition(; coordinates=shape_sampled.grid, density=1.0, particle_spacing),
+     zcolor=shape_sampled.winding_numbers)
+
+# trixi2vtk(shape_sampled.grid, w=shape_sampled.winding_numbers)

src/preprocessing/point_in_poly/hierarchical_winding.jl

@@ -2,29 +2,32 @@
 # It is implementing a binary tree and thus stores the left and right child and also the
 # faces and closing faces which are inside the bounding box.
 struct BoundingBoxTree{MC, NDIMS}
-    faces         :: Vector{Int}
+    faces         :: Vector{NTuple{NDIMS, Int}}
+    closing_faces :: Vector{NTuple{NDIMS, Int}}
     min_corner    :: MC
     max_corner    :: MC
     is_leaf       :: Bool
-    closing_faces :: Vector{NTuple{NDIMS, Int}}
     child_left    :: BoundingBoxTree
     child_right   :: BoundingBoxTree
 
-    function BoundingBoxTree(geometry, faces, directed_edges, min_corner, max_corner)
+    function BoundingBoxTree(geometry, face_ids, directed_edges, min_corner, max_corner)
+        faces = Vector{NTuple{ndims(geometry), Int}}()
         closing_faces = Vector{NTuple{ndims(geometry), Int}}()
 
+        add_faces_in_box!(faces, face_ids, geometry)
+
         max_faces_in_box = ndims(geometry) == 3 ? 100 : 20
-        if length(faces) < max_faces_in_box
+        if length(face_ids) < max_faces_in_box
             return new{typeof(min_corner),
-                       ndims(geometry)}(faces, min_corner, max_corner, true, closing_faces)
+                       ndims(geometry)}(faces, closing_faces, min_corner, max_corner, true)
         end
 
-        determine_closure!(closing_faces, min_corner, max_corner, geometry, faces,
+        determine_closure!(closing_faces, min_corner, max_corner, geometry, face_ids,
                            directed_edges)
 
-        if length(closing_faces) >= length(faces)
+        if length(closing_faces) >= length(face_ids)
             return new{typeof(min_corner),
-                       ndims(geometry)}(faces, min_corner, max_corner, true, closing_faces)
+                       ndims(geometry)}(faces, closing_faces, min_corner, max_corner, true)
         end
 
         # Bisect the box splitting its longest side
@@ -37,20 +40,47 @@ struct BoundingBoxTree{MC, NDIMS}
         max_corner_left = max_corner - 0.5box_edges[split_direction] * uvec
         min_corner_right = min_corner + 0.5box_edges[split_direction] * uvec
 
-        faces_left = is_in_box(geometry, faces, min_corner, max_corner_left)
-        faces_right = is_in_box(geometry, faces, min_corner_right, max_corner)
+        faces_left = is_in_box(geometry, face_ids, min_corner, max_corner_left)
+        faces_right = is_in_box(geometry, face_ids, min_corner_right, max_corner)
 
         child_left = BoundingBoxTree(geometry, faces_left, directed_edges,
                                      min_corner, max_corner_left)
         child_right = BoundingBoxTree(geometry, faces_right, directed_edges,
                                       min_corner_right, max_corner)
 
         return new{typeof(min_corner),
-                   ndims(geometry)}(faces, min_corner, max_corner, false, closing_faces,
+                   ndims(geometry)}(faces, closing_faces, min_corner, max_corner, false,
                                     child_left, child_right)
     end
 end
 
+function add_faces_in_box!(edges, edge_ids, geometry::Polygon)
+    (; edge_vertices_ids) = geometry
+
+    for edge in edge_ids
+        v1 = edge_vertices_ids[edge][1]
+        v2 = edge_vertices_ids[edge][2]
+
+        push!(edges, (v1, v2))
+    end
+
+    return edges
+end
+
+function add_faces_in_box!(faces, face_ids, geometry::TriangleMesh)
+    (; face_vertices_ids) = geometry
+
+    for face in face_ids
+        v1 = face_vertices_ids[face][1]
+        v2 = face_vertices_ids[face][2]
+        v3 = face_vertices_ids[face][3]
+
+        push!(faces, (v1, v2, v3))
+    end
+
+    return faces
+end
+
 struct HierarchicalWinding{BB}
     bounding_box::BB
 
@@ -214,7 +244,7 @@ function determine_closure!(closing_edges, min_corner, max_corner, polygon::Poly
         end
     end
 
-    exterior_vertices = sort(findall(!iszero, vertex_count))
+    exterior_vertices = findall(!iszero, vertex_count)
     resize!(closing_edges, 0)
 
     if !isempty(exterior_vertices)

src/preprocessing/point_in_poly/winding_number_jacobson.jl

@@ -1,21 +1,24 @@
 struct NaiveWinding end
 
 @inline function (winding::NaiveWinding)(polygon::Polygon{2}, query_point)
-    (; edge_vertices) = polygon
+    (; edge_vertices_ids) = polygon
 
-    return naive_winding(polygon, edge_vertices, query_point)
+    return naive_winding(polygon, edge_vertices_ids, query_point)
 end
 
 @inline function (winding::NaiveWinding)(mesh::TriangleMesh{3}, query_point)
-    (; face_vertices) = mesh
+    (; face_vertices_ids) = mesh
 
-    return naive_winding(mesh, face_vertices, query_point)
+    return naive_winding(mesh, face_vertices_ids, query_point)
 end
 
 @inline function naive_winding(polygon::Polygon{2}, edges, query_point)
     winding_number = sum(edges, init=zero(eltype(polygon))) do edge
-        a = edge_vertex(polygon, edge, 1) - query_point
-        b = edge_vertex(polygon, edge, 2) - query_point
+        v1 = polygon.vertices[edge[1]]
+        v2 = polygon.vertices[edge[2]]
+
+        a = v1 - query_point
+        b = v2 - query_point
 
         return atan(det([a b]), (dot(a, b)))
     end
@@ -25,12 +28,15 @@ end
 
 @inline function naive_winding(mesh::TriangleMesh{3}, faces, query_point)
     winding_number = sum(faces, init=zero(eltype(mesh))) do face
+        v1 = mesh.vertices[face[1]]
+        v2 = mesh.vertices[face[2]]
+        v3 = mesh.vertices[face[3]]
 
         # Eq. 6 of Jacobsen et al. Based on A. Van Oosterom (1983),
         # "The Solid Angle of a Plane Triangle" (doi: 10.1109/TBME.1983.325207)
-        a = face_vertex(mesh, face, 1) - query_point
-        b = face_vertex(mesh, face, 2) - query_point
-        c = face_vertex(mesh, face, 3) - query_point
+        a = v1 - query_point
+        b = v2 - query_point
+        c = v3 - query_point
         a_ = norm(a)
         b_ = norm(b)
         c_ = norm(c)
@@ -101,47 +107,3 @@ function (point_in_poly::WindingNumberJacobson)(geometry, points;
 
     return inpoly, winding_numbers
 end
-
-# The following functions distinguish between actual triangles (edges)
-# and reconstructed triangles (edges) in the hierarchical winding approach.
-# When we reconstruct the triangles then we pass the index of the vertices. Otherwise,
-# we pass the coordinates of the vertices.
-
-# This method is used, when `naive_winding` is called with `(winding::NaiveWinding)`.
-# `face` holds the coordinates of each vertex.
-@inline face_vertex(mesh, face, index) = face[index]
-
-# `edge` holds the coordinates of each vertex
-@inline edge_vertex(mesh, edge, index) = edge[index]
-
-# This method is used, when `naive_winding` is called with `(winding::HierarchicalWinding)`
-# and the query point is outside the bounding box. That is, we use the closure of the box.
-# `face` holds the index of each vertex.
-@inline function face_vertex(mesh, face::NTuple{3, Int}, index)
-    v_id = face[index]
-
-    return mesh.vertices[v_id]
-end
-
-# `edge` holds the index of each vertex
-@inline function edge_vertex(mesh, edge::NTuple{2, Int}, index)
-    v_id = edge[index]
-
-    return mesh.vertices[v_id]
-end
-
-# This method is used, when `naive_winding` is called with `(winding::HierarchicalWinding)`
-# and the bounding box is a leaf.
-# `face` is the index of the face.
-@inline function face_vertex(mesh, face::Int, index)
-    (; face_vertices) = mesh
-
-    return face_vertices[face][index]
-end
-
-# `edge` is the index of the edge
-@inline function edge_vertex(mesh, edge::Int, index)
-    (; edge_vertices) = mesh
-
-    return edge_vertices[edge][index]
-end

src/setups/complex_shape.jl

@@ -50,10 +50,13 @@ function ComplexShape(geometry::Union{TriangleMesh, Polygon}; particle_spacing,
         throw(ArgumentError("`WindingNumberHorman` only supports 2D geometries"))
     end
 
+    if grid_offset < 0.0
+        throw(ArgumentError("only a positive `grid_offset` is supported"))
+    end
+
     return sample(geometry; particle_spacing, density, pressure, mass, velocity,
                   point_in_geometry_algorithm, store_winding_number, grid_offset,
-                  max_nparticles,
-                  padding=pad_initial_particle_grid)
+                  max_nparticles, padding=pad_initial_particle_grid)
 end
 
 function sample(geometry; particle_spacing, density, pressure=0.0, mass=nothing,