Merge pull request #95 from NicolasRiel/main

Function to add stripes on shapes

Author: boriskaus

docs/src/man/geodynamic_setups.md

@@ -14,6 +14,7 @@ GeophysicalModelGenerator.AddLayer!
 GeophysicalModelGenerator.AddSphere!
 GeophysicalModelGenerator.AddEllipsoid!
 GeophysicalModelGenerator.AddCylinder!
+GeophysicalModelGenerator.addStripes!
 GeophysicalModelGenerator.makeVolcTopo
 GeophysicalModelGenerator.ConstantTemp
 GeophysicalModelGenerator.LinearTemp

src/Setup_geometry.jl

@@ -9,14 +9,121 @@ using GeoParams
 # These are routines that help to create input geometries, such as slabs with a given angle
 #
 
-export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!,
+export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addStripes!,
         makeVolcTopo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
         Compute_ThermalStructure, Compute_Phase,
         McKenzie_subducting_slab, LinearWeightedTemperature
 
 
+
+"""
+    addStripes!(Phase, Grid::AbstractGeneralGrid;
+        stripAxes       = (1,1,0),
+        stripeWidth     =  0.2,
+        stripeSpacing   =  1,
+        Origin          =  nothing,
+        StrikeAngle     =  0,
+        DipAngle        =  10,
+        phase           =  ConstantPhase(3),
+        stripePhase     =  ConstantPhase(4))
+
+    Adds stripes to a pre-defined phase (e.g. added using AddBox!)
+
+
+    Parameters
+    ====
+    - Phase - Phase array (consistent with Grid)
+    - Grid -  grid structure (usually obtained with ReadLaMEM_InputFile, but can also be other grid types)
+    - stripAxes - sets the axis for which we want the stripes. Default is (1,1,0) i.e. X, Y and not Z
+    - stripeWidth - width of the stripe
+    - stripeSpacing - space between two stripes
+    - Origin - the origin, used to rotate the box around. Default is the left-front-top corner
+    - StrikeAngle - strike angle
+    - DipAngle - dip angle
+    - phase - specifies the phase we want to apply stripes to
+    - stripePhase - specifies the stripe phase
+    
+
+    Example
+    ========
+    
+    Example: Box with striped phase and constant temperature & a dip angle of 10 degrees:
+    ```julia
+    julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
+    LaMEM Grid:
+      nel         : (32, 32, 32)
+      marker/cell : (3, 3, 3)
+      markers     : (96, 96, 96)
+      x           ϵ [-3.0 : 3.0]
+      y           ϵ [-2.0 : 2.0]
+      z           ϵ [-2.0 : 0.0]
+    julia> Phases = zeros(Int32,   size(Grid.X));
+    julia> Temp   = zeros(Float64, size(Grid.X));
+    julia> AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
+    julia> addStripes!(Phases, Grid, stripAxes=(1,1,1), stripeWidth=0.2, stripeSpacing=1, Origin=nothing, StrikeAngle=0, DipAngle=10, phase=ConstantPhase(3), stripePhase=ConstantPhase(4))
+    julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
+    julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
+    1-element Vector{String}:
+     "LaMEM_ModelSetup.vts"
+    ```
+"""
+function addStripes!(Phase, Grid::AbstractGeneralGrid;                # required input
+    stripAxes       = (1,1,0),                          # activate stripes along dimensions x, y and z when set to 1
+    stripeWidth     =  0.2,                             # full width of a stripe
+    stripeSpacing   =  1,                               # spacing between two stripes centers
+    Origin          =  nothing,                         # origin
+    StrikeAngle     =  0,                               # strike
+    DipAngle        =  0,                               # dip angle
+    phase           =  ConstantPhase(3),                # phase to be striped
+    stripePhase     =  ConstantPhase(4))                # stripe phase
+
+    # warnings
+    if stripeWidth >= stripeSpacing/2.0
+        print("WARNING: stripeWidth should be strictly < stripeSpacing/2.0, otherwise phase is overwritten by the stripePhase\n")
+    elseif sum(stripAxes .== 0) == 3
+        print("WARNING: at least one axis should be set to 1 e.g. stripAxes = (1,0,0), otherwise no stripes will be added\n")
+    end
+
+    # Retrieve 3D data arrays for the grid
+    X,Y,Z = coordinate_grids(Grid)
+
+    # sets origin
+    if isnothing(Origin)
+        Origin = (maximum(X), maximum(Y), maximum(Z))  # upper-left corner
+    end
+
+    # Perform rotation of 3D coordinates:
+    Xrot = X .- Origin[1];
+    Yrot = Y .- Origin[2];
+    Zrot = Z .- Origin[3];
+
+    Rot3D!(Xrot,Yrot,Zrot, StrikeAngle, DipAngle)
+
+    ph_ind  = findall(Phase .== phase.phase);
+
+    ind = Int64[]
+    if stripAxes[1] == 1
+        indX     = findall( abs.(Xrot[ph_ind] .% stripeSpacing) .<= stripeWidth/2.0);
+        ind = vcat(ind,indX);
+    end
+    if stripAxes[2] == 1
+        indY     = findall( abs.(Yrot[ph_ind] .% stripeSpacing) .<= stripeWidth/2.0);
+        ind = vcat(ind,indY);
+    end
+    if stripAxes[3] == 1
+        indZ     = findall( abs.(Zrot[ph_ind] .% stripeSpacing) .<= stripeWidth/2.0);
+        ind = vcat(ind,indZ);
+    end
+
+    Phase[ph_ind[ind]] .= stripePhase.phase;
+    
+    return nothing
+end
+
+
+
 """
     AddBox!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
             Origin=nothing, StrikeAngle=0, DipAngle=0,

test/test_setup_geometry.jl

@@ -34,9 +34,11 @@ Temp                =   ones(Float64, Grid.N...)*1350;
 Phases              =   zeros(Int32,  Grid.N...);
 
 AddBox!(Phases,Temp,Grid, xlim=(2,4), zlim=(4,8), phase=ConstantPhase(3), DipAngle=10, T=LinearTemp(Tbot=1350, Ttop=200))
-
 @test maximum(Phases) == 3
 
+addStripes!(Phases, Grid,stripAxes=(1,1,1),stripeWidth=0.2,stripeSpacing=1,Origin=nothing, StrikeAngle=0, DipAngle=10,phase = ConstantPhase(3),stripePhase = ConstantPhase(4))
+@test maximum(Phases) == 4
+
 # Create a CartData structure from it
 Data = CartData(Grid, (T=Temp, Phases=Phases))