Merge pull request #154 from adbayonao/Andres_Bayona_2

New function for generating tectonic plates with custom geometries and constant depth

Author: aelligp

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!, add_fault!,
+export add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_plate!, add_slab!, add_stripes!, add_volcano!, add_fault!,
     make_volc_topo,
     ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp, LinearWeightedTemperature,
     McKenzie_subducting_slab,
@@ -299,7 +299,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid;       # required input
                 Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
             end
             if segments !== nothing
-                Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T, segments)
+                Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T, segments)
             else
                 Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T)
             end
@@ -768,6 +768,91 @@ function add_polygon!(
     return nothing
 end
 
+"""
+        add_plate!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim=(), zlim::Tuple = (0.0,0.8), phase = ConstantPhase(1), T=nothing, segments=nothing, cell=false )
+Adds a tectonic plate with phase and temperature structure to a 3D model setup.
+This function enables the definition of tectonic plates in the xy plane and projects them along the z-axis, providing a flexible approach to model complex plate geometries.
+Parameters
+==========
+- `Phase`  - Phase array (consistent with Grid)
+- `Temp`   - Temperature array (consistent with Grid)
+- `Grid`   - Grid structure (usually obtained with read_LaMEM_inputfile)
+- `xlim`   - `x`-coordinate of the polygon points, same ordering as ylim, number of points unlimited
+- `ylim`   - `y`-coordinate of the polygon points, same ordering as xlim, number of points unlimited
+- `zlim`   - `z`-coordinate range for projecting the polygon (start and stop, two values)
+- `phase`  - Specifies the phase of the plate. See `ConstantPhase()`
+- `T`      - Specifies the temperature of the plate. See `ConstantTemp()`, `LinearTemp()`, `HalfspaceCoolingTemp()`, `SpreadingRateTemp()`
+- `segments` - Optional. Allows for thermal segmentation within the polygon. Useful for ridge systems or complex thermal structures.
+- `cell`   - If true, `Phase` and `Temp` are defined on cell centers
+Example
+========
+Tectonic plate in the xy plane with phase and temperature structure:
+```julia-repl
+julia> Grid = CartData(xyz_grid(x, y, z))
+Grid:
+  nel         : (512, 512, 128)
+  marker/cell : (1, 1, 1)
+  markers     : (512, 512, 128)
+  x           ϵ [-1000.0 : 0.0]
+  y           ϵ [-1000.0 : 1000.0]
+  z           ϵ [-660.0 : 0.0]
+julia> Phases = zeros(Int32,   size(Grid.X))
+julia> Temp   = zeros(Float64, size(Grid.X))
+julia> segments = [
+           ((-500.0, -1000.0), (-500.0, 0.0)),  # Segment 1
+           ((-250.0, 0.0), (-250.0, 200.0)),    # Segment 2
+           ((-750.0, 200.0), (-750.0, 1000.0))  # Segment 3
+       ]
+julia> lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
+julia> add_plate!(Phases, Temp, Grid; 
+           xlim=(-1000.0, -750.0, -250.0, 0.0, -250.0, -750.0), 
+           ylim=(0.0, 500.0, 500.0, 0.0, -500.0, -500.0), 
+           zlim=(-150.0, 0.0), 
+           phase=lith, 
+           T=SpreadingRateTemp(SpreadingVel=3), 
+           segments=segments)
+julia> Grid = addfield(Grid, (; Phases, Temp))  # Add fields
+julia> write_paraview(Grid, "Plate")  # Save model to Paraview
+1-element Vector{String}:
+ "Plate.vts"
+"""
+
+function add_plate!(Phase, Temp, Grid::AbstractGeneralGrid;   
+    xlim=(), ylim=(), zlim::Tuple = (0.0,0.8),           
+    phase = ConstantPhase(1),                                   
+    T=nothing, segments=nothing, cell=false )                            
+
+    xlim_ = collect(xlim)
+    ylim_ = collect(ylim)
+    zlim_ = collect(zlim)
+
+    X, Y, Z = coordinate_grids(Grid, cell=cell)
+    ind = zeros(Bool, size(X))
+    ind_slice = zeros(Bool, size(X[:, :, 1]))
+
+    for k = 1:size(Z, 3)
+        if zlim_[1] <= Z[1, 1, k] <= zlim_[2]
+            inpolygon!(ind_slice, xlim_, ylim_, X[:, :, k], Y[:, :, k])
+            @views ind[:, :, k] = ind_slice
+        else
+            @views ind[:, :, k] = zeros(size(X[:, :, 1]))
+        end
+    end
+
+    if !isempty(ind)
+        if T != nothing
+            if segments !== nothing
+                Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T, segments)
+            else
+                Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
+            end
+        end
+        Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
+    end
+
+    return nothing
+end
+
 """
     xrot, yrot, zrot = Rot3D(X::Number,Y::Number,Z::Number, cosStrikeAngle, sindStrikeAngle, cosDipAngle, sinDipAngle)
 
@@ -1147,7 +1232,7 @@ Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid divisio
     Adiabat = 0        # Adiabatic gradient in K/km
     MORside = "left"   # side of box where the MOR is located
     SpreadingVel = 3   # spreading velocity [cm/yr]
-    AgeRidge = 0       # Age of the ridge [Myrs
+    AgeRidge = 0       # Age of the ridge [Myrs]
     maxAge  = 60       # maximum thermal age of plate [Myrs]
 end
 

test/runtests.jl

@@ -78,6 +78,10 @@ using Test
     	include("test_ridge_segments.jl")
     end
 
+    @testset "Plate Tests" begin
+        include("test_plate.jl")
+    end
+
     @testset "Waterflow" begin
         include("test_WaterFlow.jl")
     end

test/test_plate.jl

@@ -0,0 +1,47 @@
+using GeophysicalModelGenerator
+using Test
+
+@testset "Plate Tests" begin
+    # Grid parameters
+    nx, ny, nz = 512, 512, 128
+    x = range(-1000, 0, length=nx)
+    y = range(-1000, 1000, length=ny)
+    z = range(-660, 0, length=nz)
+    Grid = CartData(xyz_grid(x, y, z))
+
+    # Phases and temperature
+    Phases = fill(2, nx, ny, nz)
+    Temp = fill(1350.0, nx, ny, nz)
+
+    # Ridge Segments
+    segments = [
+        ((-500.0, -1000.0), (-500.0, 0.0)),  # Segment 1
+        ((-250.0, 0.0), (-250.0, 200.0)),    # Segment 2
+        ((-750.0, 200.0), (-750.0, 1000.0))  # Segment 3
+    ]
+
+    # Lithospheric phases
+    lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
+
+    # Replace add_box with add_polygon_xy to allow for arbitrary-shaped ridges
+    add_plate!(Phases, Temp, Grid; 
+        xlim=(-1000.0, -750.0, -250.0, 0.0, -250.0, -750.0), 
+        ylim=(0.0, 500.0, 500.0, 0.0, -500.0, -500.0), 
+        zlim=(-150.0, 0.0), 
+        phase=lith, 
+        T=SpreadingRateTemp(SpreadingVel=3), 
+        segments=segments)
+
+    # Add and save results
+    Grid = addfield(Grid, (; Phases, Temp))
+    write_paraview(Grid, "Plate")
+
+    @test minimum(Temp) >= 0.0  # Minimum temperature
+    @test maximum(Temp) <= 1350.0  # Maximum temperature
+    @test all(≥(0),    Temp) # Ensure no negative temperatures
+    @test all(≤(1350), Temp) # Ensure no temperatures above max
+
+    # Check if phases are correctly assigned in expected regions
+    @test first(Phases) == 2  # Example: Verify a point's phase
+    @test last(Phases)  == 2  # Example: Verify another point's phase.
+end