Merge pull request #149 from JuliaGeodynamics/pa-volcano

New function: `add_fault!()` in 3D

src/Setup_geometry.jl

@@ -11,7 +11,7 @@ import Base: show
 # These are routines that help to create input geometries, such as slabs with a given angle
 #
 
-export  add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano!,
+export  add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano!, add_fault!,
         make_volc_topo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp, LinearWeightedTemperature,
         McKenzie_subducting_slab,
@@ -37,6 +37,24 @@ function flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})
     return ind2D
 end
 
+"""
+    ind2D = flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})
+
+This converts the indices to purely 2D indices if the array `phase` is 2D
+"""
+function flatten_index_dimensions(Phase::AbstractArray{T, N}, ind_vec::Array{Bool, 3}) where {T,N}
+    if N==2
+        ind2D = Vector{CartesianIndex{2}}(undef,length(ind_vec))
+        for (num, ind) in enumerate(ind_vec)
+            ind2D[num] = CartesianIndex(ind[1], ind[3])
+        end
+    else
+        ind2D = ind_vec
+    end
+
+    return ind2D
+end
+
 """
     add_stripes!(Phase, Grid::AbstractGeneralGrid;
         stripAxes       = (1,1,0),
@@ -792,7 +810,7 @@ function add_volcano!(
         end
     end
 
-    ind_flat = flatten_index_dimensions(Temp, ind)
+    ind_flat = flatten_index_dimensions(Phases, ind)
 
     # @views Temp[ind .== false] .= 0.0
     Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T)
@@ -1929,3 +1947,103 @@ function add_slab!(Phase, Temp, Grid::AbstractGeneralGrid,  trench::Trench;
 
     return nothing
 end
+
+"""
+    add_fault!(Phase, Temp, Grid::AbstractGeneralGrid;
+        Start=(20,100), End=(10,80),
+        Fault_thickness=10.0,
+        Depth_extent=nothing,
+        DipAngle=0e0,
+        phase=ConstantPhase(1),
+        T=nothing,
+        cell=false)
+
+Adds a fault to the given 3D grid by modifying the `Phase` and `Temp` arrays.
+For a 2D grid, use `add_box` instead.
+
+# Arguments
+- `Phase`: Phase array
+- `Temp`: Temp array
+- `Grid`: The grid on which the fault is to be added.
+- `Start`: Tuple representing the starting coordinates of the fault (X, Y).
+- `End`: Tuple representing the ending coordinates of the fault (X, Y).
+- `Fault_thickness`: Thickness of the fault.
+- `Depth_extent`: Depth extent of the fault. If `nothing`, the fault extends through the entire domain.
+- `DipAngle`: Dip angle of the fault.
+- `phase`: Phase to be assigned to the fault.
+- `T`: Temperature to be assigned to the fault. If `nothing`, the temperature is not modified.
+
+
+# Example
+```julia
+add_fault!(Phase, Temp, Grid;
+        Start=(20,100), End=(10,80),
+        Fault_thickness=10.0,
+        Depth_extent=(-25.0, 0.0),
+        DipAngle=-10.0,
+        phase=ConstantPhase(1)
+        )
+```
+"""
+function add_fault!(
+    Phase,
+    Temp,
+    Grid::AbstractGeneralGrid;
+    Start=(20,100),
+    End=(10,80),
+    Fault_thickness=10.0,
+    Depth_extent=nothing,
+    DipAngle=0e0,
+    phase=ConstantPhase(1),
+    T=nothing,
+    cell=false
+)
+
+   # Extract the coordinates
+   X, Y, Z = coordinate_grids(Grid, cell=cell)
+
+   # Calculate the direction vector from Start to End
+   direction = (End[1] - Start[1], End[2] - Start[2])
+   length = sqrt(direction[1]^2 + direction[2]^2)
+   unit_direction = (direction[1], direction[2]) ./ length
+
+   # Calculate the fault region based on fault thickness and length
+   fault_half_thickness = Fault_thickness / 2
+
+   # Create a mask for the fault region
+   fault_mask = falses(size(X))
+
+    for k in 1:size(Z, 3), j in 1:size(Y, 2), i in 1:size(X, 1)
+        # Rotate the point using the dip angle
+        x_rot, y_rot, z_rot = Rot3D(X[i, j, k], Y[i, j, k], Z[i, j, k], 1.0, 0.0, cosd(DipAngle), sind(DipAngle))
+
+        # Calculate the projection of the rotated point onto the fault line
+        projection_length = (x_rot - Start[1]) * unit_direction[1] + (y_rot - Start[2]) * unit_direction[2]
+        if 0 ≤ projection_length ≤ length
+            # Calculate the perpendicular distance to the fault line
+            perpendicular_distance = abs((x_rot - Start[1]) * unit_direction[2] - (y_rot - Start[2]) * unit_direction[1])
+            if perpendicular_distance ≤ fault_half_thickness
+                fault_mask[i, j, k] = true
+            end
+        end
+    end
+
+   ind = findall(fault_mask)
+
+   # Apply depth extent if provided
+   if !isnothing(Depth_extent)
+       ind = ind[Z[ind] .≥ Depth_extent[1] .&& Z[ind] .≤ Depth_extent[2]]
+   end
+
+   ind_flat = flatten_index_dimensions(Phase, ind)
+
+   # Compute thermal structure accordingly
+   if T != nothing
+       Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
+   end
+
+   # Set the phase
+   Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
+
+   return nothing
+end