Merge branch 'main' into sgs

Author: svchb

Project.toml

@@ -47,7 +47,7 @@ DelimitedFiles = "1"
 DiffEqCallbacks = "4"
 FastPow = "0.1"
 FileIO = "1"
-ForwardDiff = "0.10"
+ForwardDiff = "1"
 GPUArraysCore = "0.2"
 JSON = "0.21"
 KernelAbstractions = "0.9"

docs/src/preprocessing/preprocessing.md

@@ -250,6 +250,51 @@ Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_hormann.jl")
 Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_jacobson.jl")]
 ```
+# [Read geometries from file](@id read_geometries_from_file)
+Geometries can be imported using the [`load_geometry`](@ref) function.
+- For 3D geometries, we support the binary (`.stl`) format.
+- For 2D geometries, the recommended format is DXF (`.dxf`), with optional support for a simple ASCII (`.asc`) format.
+
+## ASCII Format (.asc)
+An .asc file contains a list of 2D coordinates, space-delimited, one point per line,
+where the first column are the x values and the second the y values.
+For example:
+```plaintext
+# ASCII
+0.0 0.0
+1.0 0.0
+1.0 1.0
+0.0 1.0
+```
+It is the user’s responsibility to ensure the points are ordered correctly.
+This format is easy to generate and inspect manually.
+
+## DXF Format (.dxf) – recommended
+The DXF (Drawing Exchange Format) is a widely-used CAD format for 2D and 3D data.
+We recommend this format for defining 2D polygonal.
+Only DXF entities of type `LINE` or `POLYLINE` are supported.
+To create DXF files from scratch, we recommend using the open-source tool [FreeCAD](https://www.freecad.org/).
+For a less technical approach, we recommend using [Inkscape](https://inkscape.org/) to create SVG files and then export them as DXF.
+
+### Creating DXF Files with FreeCAD
+
+- Open FreeCAD and create a new document.
+- Switch to the Sketcher workbench and create a new sketch.
+- Choose the XY-plane orientation and draw your geometry.
+- Select the object to be exported.
+- Go to "File > Export..." and choose "AutoDesk DXF (*.dxf)" as the format.
+- Ensure the following Import-Export options are enabled:
+    - "Use legacy Python exporter".
+    - "Project exported objects along current view direction".
+
+### Creating DXF Files with Inkscape
+SVG vector graphics can also be used as a basis for geometry.
+Use Inkscape to open or create the SVG.
+You can simply draw a Bezier curve by pressing "B" on your keyboard.
+We reommend the mode "Create spiro paths" for a smooth curve.
+Select the relevant path and export it as DXF:
+- Go to "File > Save As...".
+- Choose "Desktop Cutting Plotter (AutoCAD DXF R12)(*.dxf)" format.
 
 ```@autodocs
 Modules = [TrixiParticles]

src/preprocessing/geometries/io.jl

@@ -2,7 +2,9 @@
     load_geometry(filename; element_type=Float64)
 
 Load file and return corresponding type for [`ComplexShape`](@ref).
-Supported file formats are `.stl` and `.asc`.
+Supported file formats are `.stl`, `.asc` and `dxf`.
+For comprehensive information about the supported file formats, refer to the documentation at
+[Read geometries from file](@ref read_geometries_from_file).
 
 # Arguments
 - `filename`: Name of the file to be loaded.
@@ -17,10 +19,12 @@ function load_geometry(filename; element_type=Float64)
 
     if file_extension == ".asc"
         geometry = load_ascii(filename; ELTYPE, skipstart=1)
+    elseif file_extension == ".dxf"
+        geometry = load_dxf(filename; ELTYPE)
     elseif file_extension == ".stl"
         geometry = load(FileIO.query(filename); ELTYPE)
     else
-        throw(ArgumentError("Only `.stl` and `.asc` files are supported (yet)."))
+        throw(ArgumentError("Only `.stl`, `.asc` and `.dxf` files are supported (yet)."))
     end
 
     return geometry
@@ -35,6 +39,97 @@ function load_ascii(filename; ELTYPE=Float64, skipstart=1)
     return Polygon(copy(points'))
 end
 
+function load_dxf(filename; ELTYPE=Float64)
+    points = Tuple{ELTYPE, ELTYPE}[]
+
+    load_dxf!(points, filename)
+
+    return Polygon(stack(points))
+end
+
+function load_dxf!(points::Vector{Tuple{T, T}}, filename) where {T}
+    open(filename) do io
+        lines = readlines(io)
+
+        # Check if the DXF file contains line entities ("LINE") or polyline entities ("POLYLINE")
+        idx_first_line = findfirst(x -> strip(x) == "LINE", lines)
+        idx_first_polyline = findfirst(x -> strip(x) == "POLYLINE", lines)
+
+        if !isnothing(idx_first_line) && isnothing(idx_first_polyline)
+            # The file contains only simple line entities ("LINE")
+            i = idx_first_line
+
+            while i <= length(lines)
+                if strip(lines[i]) == "LINE"
+                    if idx_first_line == i
+                        # For a polygon, we only need to store the start point of the first line entity.
+                        # For subsequent lines, the end point of the previous edge is the start point of the current edge,
+                        # so storing all start points would result in duplicate vertices.
+                        # Therefore, we only push the start point for the very first line.
+
+                        # Search for the coordinates of the start point:
+                        # Group codes "10", "20", "30" represent x, y, z of the start point
+                        while i <= length(lines) && strip(lines[i]) != "10"
+                            i += 1
+                        end
+                        x1 = parse(T, strip(lines[i + 1]))
+                        @assert strip(lines[i + 2]) == "20"
+                        y1 = parse(T, strip(lines[i + 3]))
+
+                        push!(points, (x1, y1))
+                    end
+
+                    # Search for the end point coordinates:
+                    # Group codes "11", "21", "31" represent x, y, z of the end point
+                    while i <= length(lines) && strip(lines[i]) != "11"
+                        i += 1
+                    end
+                    x2 = parse(T, strip(lines[i + 1]))
+                    @assert strip(lines[i + 2]) == "21"
+                    y2 = parse(T, strip(lines[i + 3]))
+
+                    # Add end point of the line to the point list
+                    push!(points, (x2, y2))
+                end
+                i += 1
+            end
+
+        elseif isnothing(idx_first_line) && !isnothing(idx_first_polyline)
+            # The file contains only polyline entities ("POLYLINE")
+            i = idx_first_polyline
+
+            while i <= length(lines)
+                line = strip(lines[i])
+                if line == "VERTEX"
+                    # Search for the coordinates of the current vertex:
+                    # Group codes "10", "20", "30" represent x, y, z of the vertex
+                    while i <= length(lines) && strip(lines[i]) != "10"
+                        i += 1
+                    end
+                    x = parse(T, strip(lines[i + 1]))
+                    @assert strip(lines[i + 2]) == "20"
+                    y = parse(T, strip(lines[i + 3]))
+
+                    # Add the vertex to the point list
+                    push!(points, (x, y))
+
+                elseif line == "SEQEND"
+                    # End of the polyline reached
+                    break
+                end
+                i += 1
+            end
+        else
+            throw(ArgumentError("This entity type is not supported. Only 'LINE' OR 'POLYLINE' are allowed."))
+        end
+    end
+
+    # Remove duplicate points from the list
+    unique(points)
+
+    return points
+end
+
 # FileIO.jl docs:
 # https://juliaio.github.io/FileIO.jl/stable/implementing/#All-at-once-I/O:-implementing-load-and-save
 function load(fn::FileIO.File{FileIO.format"STL_BINARY"}; element_types...)

test/general/smoothing_kernels.jl

@@ -86,7 +86,7 @@
 
                         # This should work with very tight tolerances
                         @test isapprox(analytic_deriv, automatic_deriv,
-                                       rtol=5e-15, atol=2eps())
+                                       rtol=5e-15, atol=4 * eps())
                     end
                 end
             end