ScatteredInterpolation gives the wrong answer, so replaced that with distance point/line

Author: boriskaus

Project.toml

@@ -22,7 +22,6 @@ MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-ScatteredInterpolation = "3f865c0f-6dca-5f4d-999b-29fe1e7e3c92"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -54,7 +53,6 @@ JLD2 = "0.4"
 MeshIO = "0.1 - 0.4"
 NearestNeighbors = "0.2 - 0.4"
 Parameters = "0.9 - 0.12"
-ScatteredInterpolation = "0.3"
 SpecialFunctions = "1.0, 2"
 StaticArrays = "1"
 WriteVTK = "1"

src/Setup_geometry.jl

@@ -3,7 +3,6 @@ using Printf
 using Parameters        # helps setting default parameters in structures
 using SpecialFunctions: erfc
 using GeoParams
-using ScatteredInterpolation
 using StaticArrays
 
 # Setup_geometry
@@ -1475,7 +1474,7 @@ function find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench::Trench)
     xb = Rot3D(End[1]-Start[1],End[2]-Start[2], 0.0, cosd(StrikeAngle), sind(StrikeAngle), 1.0, 0.0)
 
     # dl
-    dl = Length/n_seg;
+    dl = trench.Length/n_seg;
     l = 0  # length at the trench position
     D = @SVector [Top[1,2], Bottom[1,2], Bottom[1,2],Top[1,2] ]
 
@@ -1510,20 +1509,11 @@ function find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench::Trench)
         # indexes of the segment
         ind_seg = ind_s[ind]
 
-        # Prepare the variable to interpolate {I put here because to allow also a variation of thickness of the slab}
-        L = @SVector [l,l,ln,ln];
-
-        # Interpolations
-        points = @SMatrix [pa[1] pb[1] pc[1] pd[1]; pa[2] pb[2] pc[2] pd[2]]
-
-        itp1 = ScatteredInterpolation.interpolate(Shepard(), points, D);
-        itp2 = ScatteredInterpolation.interpolate(Shepard(), points, L);
-        
         # Loop over the chosen particles and interpolate the current value of L and D.
         for ip in ind_seg
-            point_ = [Yrot[ip], Z[ip]];
-            d[ip]  = ScatteredInterpolation.evaluate(itp1,point_)[1];
-            ls[ip] = ScatteredInterpolation.evaluate(itp2,point_)[1];
+            point_ = (Yrot[ip], Z[ip]);
+            d[ip]  = -distance_to_linesegment(point_, pa, pd)
+            ls[ip]  = distance_to_linesegment(point_, pb, pa) + l
         end
 
         #Update l
@@ -1532,6 +1522,40 @@ function find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench::Trench)
 end
 
 
+"""
+    distance_to_linesegment(p::NTuple{2,_T}, v::NTuple{2,_T}, w::NTuple{2,_T})
+
+Computes the distance normal distance from a point `p` to a line segment defined by the points `v` and `w`.
+"""
+function distance_to_linesegment(p::NTuple{2,_T}, v::NTuple{2,_T}, w::NTuple{2,_T})  where _T<:Number
+    dx = w[1] - v[1]
+    dy = w[2] - v[2]
+    l2 = dx*dx + dy*dy  # i.e. |w-v|^2 -  avoid a sqrt
+    if l2 == 0.0
+        dx = p[1] - v[1]
+        dy = p[2] - v[2]
+        return sqrt(dx*dx + dy*dy)   # v == w case
+    end
+    # Consider the line extending the segment, parameterized as v + t (w - v).
+    # We find projection of point p onto the line. 
+    # It falls where t = [(p-v) . (w-v)] / |w-v|^2
+    t = ((p[1] - v[1])*dx + (p[2] - v[2])*dy) / l2
+    if t < 0.0
+        dx = p[1] - v[1]
+        dy = p[2] - v[2]
+        return sqrt(dx*dx + dy*dy)       # Beyond the 'v' end of the segment
+    elseif t > 1.0
+        dx = p[1] - w[1]
+        dy = p[2] - w[2]
+        return sqrt(dx*dx + dy*dy)  # Beyond the 'w' end of the segment
+    end
+    projection_x = v[1] + t * dx
+    projection_y = v[2] + t * dy
+    dx = p[1] - projection_x
+    dy = p[2] - projection_y
+    return sqrt(dx*dx + dy*dy)
+end
+
 """
     addSlab!(Phase, Temp, Grid::AbstractGeneralGrid,  trench::Trench; phase = ConstantPhase(1), T = nothing)
 

test/test_setup_geometry.jl

@@ -237,7 +237,7 @@ TsHC = HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=30, Adiabat=0.4)
 temp = TsHC;
 
 addSlab!(Phase,Temp,Cart, t1, phase=phase, T = TsHC)
-@test Temp[84,84,110]  ≈ 1142.25814954244
+@test Temp[84,84,110]  ≈ 1045.1322688510577
 @test extrema(Phase) == (1, 4)
 
 # with weak zone
@@ -264,12 +264,10 @@ phase = LithosphericPhases(Layers=[5 7 88], Phases = [2 3 4], Tlab=nothing)
 t1 = Trench(Start = (400.0,400.0), End = (800.0,800.0),θ_max = 90.0, direction = 1.0, n_seg = 50, Length = 600.0, Thickness = 80.0, Lb = 500.0,d_decoupling = 100.0, type_bending =:Ribe,   WeakzoneThickness=10, WeakzonePhase=9)
 
 addSlab!(Phase,Temp,Cart, t1, phase=phase, T = T_slab)
-@test Temp[84,84,110]  ≈ 718.8406936737412
+@test Temp[84,84,110]  ≈ 624.6682008876219
 
 Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp)) 
 
-
-
 # 2D slab:
 nx,nz = 512,128
 x = range(-1000,1000, nx);