diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 4ac02106..ba29e034 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -33,7 +33,7 @@ function flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}}) else ind2D = ind_vec end - + return ind2D end @@ -69,7 +69,7 @@ Parameters Example ======== - + Example: Box with striped phase and constant temperature & a dip angle of 10 degrees: ```julia julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat") @@ -99,7 +99,7 @@ function add_stripes!(Phase, Grid::AbstractGeneralGrid; # require DipAngle = 0, # dip angle phase = ConstantPhase(3), # phase to be striped stripePhase = ConstantPhase(4), # stripe phase - cell = false ) # if true, Phase and Temp are defined on cell centers + cell = false ) # if true, Phase and Temp are defined on cell centers # warnings if stripeWidth >= stripeSpacing/2.0 @@ -140,7 +140,7 @@ function add_stripes!(Phase, Grid::AbstractGeneralGrid; # require end Phase[ph_ind[ind]] .= stripePhase.phase; - + return nothing end @@ -151,7 +151,7 @@ end Origin=nothing, StrikeAngle=0, DipAngle=0, phase = ConstantPhase(1), T=nothing, - cell=false ) + cell=false ) Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -210,7 +210,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input Origin=nothing, StrikeAngle=0, DipAngle=0, # origin & dip/strike phase = ConstantPhase(1), # Sets the phase number(s) in the box T=nothing, # Sets the thermal structure (various functions are available) - cell=false ) # if true, Phase and Temp are defined on cell centers + cell=false ) # if true, Phase and Temp are defined on cell centers # Retrieve 3D data arrays for the grid X,Y,Z = coordinate_grids(Grid, cell=cell) @@ -251,15 +251,15 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input ind_flat = flatten_index_dimensions(Phase, ind) # Compute thermal structure accordingly. See routines below for different options - if T != nothing + if T != nothing if isa(T,LithosphericTemp) Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end - # Set the phase. Different routines are available for that - see below. - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + # Set the phase. Different routines are available for that - see below. + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end @@ -268,7 +268,7 @@ end """ add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (1,100), [ylim::Tuple = (0,20)], zlim::Tuple = (0,-100), phase = ConstantPhase(1), - T=nothing, cell=false ) + T=nothing, cell=false ) Adds a layer with phase & temperature structure to a 3D model setup. The most common use would be to add a lithospheric layer to a model setup. @@ -325,8 +325,8 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input xlim=nothing, ylim=nothing, zlim=nothing, # limits of the layer phase = ConstantPhase(1), # Sets the phase number(s) in the box T=nothing, # Sets the thermal structure (various functions are available) - cell = false ) # if true, Phase and Temp are defined on cell centers - + cell = false ) # if true, Phase and Temp are defined on cell centers + # Retrieve 3D data arrays for the grid X,Y,Z = coordinate_grids(Grid, cell=cell) @@ -369,9 +369,9 @@ end """ - add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen::NTuple{3, _T} = (0,0,-1), radius::Number, + add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen::Tuple = (0,0,-1), radius::Number, phase = ConstantPhase(1). - T=nothing, cell=false ) where _T + T=nothing, cell=false ) Adds a sphere with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -412,9 +412,9 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para ``` """ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - cen::NTuple{3, _T} = (0,0,-1), radius::Number, # center and radius of the sphere + cen::Tuple = (0,0,-1), radius::Number, # center and radius of the sphere phase = ConstantPhase(1), # Sets the phase number(s) in the sphere - T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) + T=nothing, cell=false ) # Sets the thermal structure (various functions are available) # Retrieve 3D data arrays for the grid X,Y,Z = coordinate_grids(Grid, cell=cell) @@ -436,10 +436,10 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input end """ - add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen::NTuple{3, _T} = (-1,-1,-1), axes::NTuple{3, _T} = (0.2,0.1,0.5), + add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen::Tuple = (-1,-1,-1), axes::Tuple = (0.2,0.1,0.5), Origin=nothing, StrikeAngle=0, DipAngle=0, phase = ConstantPhase(1). - T=nothing, cell=false ) where _T + T=nothing, cell=false ) Adds an Ellipsoid with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -481,10 +481,10 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para ``` """ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - cen::NTuple{3, _T} = (-1,-1,-1), axes::NTuple{3, _T} = (0.2,0.1,0.5), # center and semi-axes of the ellpsoid + cen::Tuple = (-1,-1,-1), axes::Tuple = (0.2,0.1,0.5), # center and semi-axes of the ellpsoid Origin=nothing, StrikeAngle=0, DipAngle=0, # origin & dip/strike phase = ConstantPhase(1), # Sets the phase number(s) in the box - T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) + T=nothing, cell=false ) # Sets the thermal structure (various functions are available) if Origin==nothing Origin = cen # center @@ -522,9 +522,9 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i end """ - add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::NTuple{3} = (-1,-1,-1.5), cap::NTuple{3} = (-1,-1,-0.5), radius::Number, + add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::Tuple = (-1,-1,-1.5), cap::Tuple = (-1,-1,-0.5), radius::Number, phase = ConstantPhase(1), - T=nothing, cell=false ) + T=nothing, cell=false ) Adds a cylinder with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -566,7 +566,7 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para ``` """ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - base::NTuple{3} = (-1,-1,-1.5), cap::NTuple{3} = (-1,-1,-0.5), radius::Number, # center and radius of the sphere + base::Tuple = (-1,-1,-1.5), cap::Tuple = (-1,-1,-0.5), radius::Number, # center and radius of the sphere phase = ConstantPhase(1), # Sets the phase number(s) in the sphere T=nothing, cell=false ) # Sets the thermal structure (various functions are available) @@ -622,7 +622,7 @@ end """ - add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::Tuple = (0.0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false ) + add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::Tuple = (0.0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false ) Adds a polygon with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -682,7 +682,7 @@ ind_slice = zeros(Bool,size(X[:,1,:])) # find points within the polygon, only in 2D for i = 1:size(Y)[2] - if Y[1,i,1] >= ylim_[1] && Y[1,i,1]<=ylim_[2] + if Y[1,i,1] >= ylim_[1] && Y[1,i,1]<=ylim_[2] inpolygon!(ind_slice, xlim_,zlim_, X[:,i,:], Z[:,i,:]) ind[:,i,:] = ind_slice else @@ -717,7 +717,7 @@ function Rot3D(X::_T,Y::_T,Z::_T, cosStrikeAngle::_T, sindStrikeAngle::_T, cosDi CoordVec = @SVector [X, Y, Z] CoordRot = rotz*CoordVec; CoordRot = roty*CoordRot; - + return CoordRot[1], CoordRot[2], CoordRot[3] end @@ -754,8 +754,8 @@ Optional Parameters - background - this allows loading in a topography and only adding the volcano on top (also allows stacking of several cones to get a volcano with different slopes) """ function add_volcano!( - Phases, - Temp, + Phases, + Temp, Grid::CartData; volcanic_phase = 1, center = (0,0,0), @@ -781,7 +781,7 @@ function add_volcano!( for k in axes(ind, 3) for j in axes(ind, 2), i in axes(ind, 1) - depth[i, j, k] = max(H[i, j] - Grid.z.val[i, j, k], 0) + depth[i, j, k] = max(H[i, j] - Grid.z.val[i, j, k], 0) if Grid.z.val[i, j, k] < H[i, j] && Grid.z.val[i, j, k] ≥ base Phases[i, j, k] = volcanic_phase @@ -796,7 +796,7 @@ function add_volcano!( # @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) - + return nothing end @@ -1117,14 +1117,14 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp) end """ - LithosphericTemp(Tsurface=0.0, Tpot=1350.0, dTadi=0.5, - ubound="const", lbound="const, utbf = 50.0e-3, ltbf = 10.0e-3, - age = 120.0, dtfac = 0.9, nz = 201, - rheology = example_CLrheology() + LithosphericTemp(Tsurface=0.0, Tpot=1350.0, dTadi=0.5, + ubound="const", lbound="const, utbf = 50.0e-3, ltbf = 10.0e-3, + age = 120.0, dtfac = 0.9, nz = 201, + rheology = example_CLrheology() ) -Calculates a 1D temperature profile [C] for variable thermal parameters including radiogenic heat source and - linearly interpolates the temperature profile onto the box. The thermal parameters are defined in +Calculates a 1D temperature profile [C] for variable thermal parameters including radiogenic heat source and + linearly interpolates the temperature profile onto the box. The thermal parameters are defined in rheology and the structure of the lithosphere is define by LithosphericPhases(). @@ -1133,7 +1133,7 @@ Parameters - Tsurface : surface temperature [C] - Tpot : potential mantle temperature [C] - dTadi : adiabatic gradient [K/km] -- ubound : Upper thermal boundary condition ["const","flux"] +- ubound : Upper thermal boundary condition ["const","flux"] - lbound : Lower thermal boundary condition ["const","flux"] - utbf : Upper thermal heat flux [W/m]; if ubound == "flux" - ltbf : Lower thermal heat flux [W/m]; if lbound == "flux" @@ -1153,7 +1153,7 @@ Parameters ltbf = 10.0e-3 # q [W/m^2]; if lbound = "flux" age = 120.0 # Lithospheric age [Ma] dtfac = 0.9 # Diffusion stability criterion - nz = 201 + nz = 201 rheology = example_CLrheology() end @@ -1174,20 +1174,20 @@ struct Thermal_parameters{A} end function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp) - @unpack Tsurface, Tpot, dTadi, ubound, lbound, utbf, ltbf, age, + @unpack Tsurface, Tpot, dTadi, ubound, lbound, utbf, ltbf, age, dtfac, nz, rheology = s # Create 1D depth profile within the box z = LinRange(round(maximum(Z)),round(minimum(Z)),nz) # [km] - z = @. z*1e3 # [m] + z = @. z*1e3 # [m] dz = z[2] - z[1] # Gride resolution # Initialize 1D arrays for explicit solver - T = zeros(nz) + T = zeros(nz) phase = Int64.(zeros(nz)) # Assign phase id from Phase to 1D phase array - phaseid = (minimum(Phase):1:maximum(Phase)) + phaseid = (minimum(Phase):1:maximum(Phase)) ztop = round(maximum(Z[findall(Phase .== phaseid[1])])) zlayer = zeros(length(phaseid)) for i = 1:length(phaseid) @@ -1197,20 +1197,20 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp) end for i = 1:length(phaseid) # Assign phase ids - ind = findall((z .>= zlayer[i]) .& (z .<= ztop)) + ind = findall((z .>= zlayer[i]) .& (z .<= ztop)) phase[ind] .= phaseid[i] ztop = zlayer[i] end # Setup initial T-profile - Tpot = Tpot + 273.15 # Potential temp [K] + Tpot = Tpot + 273.15 # Potential temp [K] Tsurface = Tsurface + 273.15 # Surface temperature [ K ] - T = @. Tpot + abs.(z./1.0e3)*dTadi # Initial T-profile [ K ] - T[1] = Tsurface - + T = @. Tpot + abs.(z./1.0e3)*dTadi # Initial T-profile [ K ] + T[1] = Tsurface + args = (;) - thermal_parameters = Thermal_parameters(nz) - + thermal_parameters = Thermal_parameters(nz) + ## Update thermal parameters ======================================== # compute_density!(thermal_parameters.ρ,rheology,phase,args) compute_heatcapacity!(thermal_parameters.Cp,rheology,phase,args) @@ -1219,7 +1219,7 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp) compute_radioactive_heat!(thermal_parameters.H,rheology,phase,args) # Thermal diffusivity [ m^2/s ] - κ = maximum(thermal_parameters.k) / + κ = maximum(thermal_parameters.k) / minimum(thermal_parameters.ρ) / minimum(thermal_parameters.Cp) ## =================================================================== # ## Time stability criterion ========================================= # @@ -1229,7 +1229,7 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp) dt = dtfac*dtexp # [s] nit = Int64(ceil(age/dt)) # Number of iterations time = zeros(nit) # Time array - + for i = 1:nit if i > 1 time[i] = time[i-1] + dt @@ -1245,18 +1245,18 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp) interp_linear_T = linear_interpolation(-z./1.0e3, T.-273.15) # create interpolation object Temp = interp_linear_T(-Z) - + return Temp end -function SolveDiff1Dexplicit_vary!( +function SolveDiff1Dexplicit_vary!( T, thermal_parameters, ubound,lbound, utbf,ltbf, di, dt -) +) nz = length(T) T0 = T @@ -1266,9 +1266,9 @@ function SolveDiff1Dexplicit_vary!( kB = (thermal_parameters.k[2] + thermal_parameters.k[1])/2.0 kA = (thermal_parameters.k[1] + thermal_parameters.k[1])/2.0 a = (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[1]) - b = 1 - (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[1]) + b = 1 - (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[1]) c = (dt*2.0*utbf)/(di * thermal_parameters.ρCp[1]) - T[1] = a*T0[2] + b*T0[1] + c + + T[1] = a*T0[2] + b*T0[1] + c + thermal_parameters.H[1]*dt/thermal_parameters.ρCp[1] end if lbound == "const" @@ -1287,12 +1287,12 @@ function SolveDiff1Dexplicit_vary!( ai = @. (kBi*dt)/(di^2.0*thermal_parameters.ρCp[2:end-1]) bi = @. 1.0 - (dt*(kAi + kBi))/(di^2.0*thermal_parameters.ρCp[2:end-1]) ci = @. (kAi*dt)/(di^2.0*thermal_parameters.ρCp[2:end-1]) - T[2:end-1] = @. ai*T0[3:end] + bi*T0[2:end-1] + ci*T0[1:end-2] + + T[2:end-1] = @. ai*T0[3:end] + bi*T0[2:end-1] + ci*T0[1:end-2] + thermal_parameters.H[2:end-1]*dt/thermal_parameters.ρCp[2:end-1] - return T + return T end -function example_CLrheology(; +function example_CLrheology(; ρM=3.0e3, # Density [ kg/m^3 ] CpM=1.0e3, # Specific heat capacity [ J/kg/K ] kM=2.3, # Thermal conductivity [ W/m/K ] @@ -1305,7 +1305,7 @@ function example_CLrheology(; CpLC=1.0e3, # Specific heat capacity [ J/kg/K ] kLC=2.0, # Thermal conductivity [ W/m/K ] HLC=43.0e-12, # Radiogenic heat source per mass [H] = W/kg; [H] = [Q/rho] -) +) rheology = ( # Name = "UpperCrust", @@ -1448,37 +1448,37 @@ Parameters it::Int64 = 36 # number of harmonic summation (look Mckenzie formula) end -""" +""" compute_thermal_structure(Temp, X, Y, Z, Phase, s::McKenzie_subducting_slab) Compute the temperature field of a `McKenzie_subducting_slab`. Uses the analytical solution of McKenzie (1969) ["Speculations on the consequences and causes of plate motions"]. The functions assumes -that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate. +that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate. Parameters ============================= Temp Temperature array -- `X` X Array -- `Y` Y Array -- `Z` Z Array -- `Phase` Phase array +- `X` X Array +- `Y` Y Array +- `Z` Z Array +- `Phase` Phase array - `s` McKenzie_subducting_slab """ function compute_thermal_structure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_slab) @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s - # Thickness of the layer: + # Thickness of the layer: Thickness = (maximum(Z)-minimum(Z)); Zshift = Z .- Z[end] # McKenzie model is defined with Z = 0 at the bottom of the slab - # Convert subduction velocity from cm/yr -> m/s; + # Convert subduction velocity from cm/yr -> m/s; convert_velocity = 1/(100.0*365.25*60.0*60.0*24.0); v_s = v_cm_yr*convert_velocity; - + # calculate the thermal Reynolds number Re = (v_s*Thickness*1000)/2/κ; # factor 1000 to transfer Thickness from km to m - + # McKenzie model sc = 1/Thickness σ = ones(size(Temp)); @@ -1488,12 +1488,12 @@ function compute_thermal_structure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_s b = (Re .- (Re.^2 .+ i^2.0 .* pi^2.0).^(0.5)) .*X .*sc; c = sin.(i .*pi .*(1 .- abs.(Zshift .*sc))) ; e = exp.(b); - σ .+= 2*a.*e.*c + σ .+= 2*a.*e.*c end Temp .= Tsurface .+ (Tmantle-Tsurface).*σ; Temp .= Temp + (Adiabat*abs.(Z)) - + return Temp end @@ -1512,41 +1512,41 @@ Parameters - F2: Second temperature field """ -@with_kw_noshow mutable struct LinearWeightedTemperature <: AbstractThermalStructure - w_min::Float64 = 0.0; - w_max::Float64 = 1.0; +@with_kw_noshow mutable struct LinearWeightedTemperature <: AbstractThermalStructure + w_min::Float64 = 0.0; + w_max::Float64 = 1.0; crit_dist::Float64 = 100.0; - dir::Symbol =:X; + dir::Symbol =:X; F1::AbstractThermalStructure = ConstantTemp(); F2::AbstractThermalStructure = ConstantTemp(); end """ compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature) - + Weight average along distance -Do a weight average between two field along a specified direction -Given a distance {could be any array, from X,Y} -> it increase from the origin the weight of -F1, while F2 decreases. -This function has been conceived for averaging the solution of Mckenzie and half space cooling model, but in -can be used to smooth the temperature field from continent ocean: --> Select the boundary to apply; +Do a weight average between two field along a specified direction +Given a distance {could be any array, from X,Y} -> it increase from the origin the weight of +F1, while F2 decreases. +This function has been conceived for averaging the solution of Mckenzie and half space cooling model, but in +can be used to smooth the temperature field from continent ocean: +-> Select the boundary to apply; -> transform the coordinate such that dist represent the perpendicular direction along which you want to apply -this smoothening and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0; --> Select the points that belongs to this area -> compute the thermal fields {F1} {F2} -> then modify F. +this smoothening and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0; +-> Select the points that belongs to this area -> compute the thermal fields {F1} {F2} -> then modify F. """ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature) - @unpack w_min, w_max, crit_dist,dir = s; - @unpack F1, F2 = s; - + @unpack w_min, w_max, crit_dist,dir = s; + @unpack F1, F2 = s; + if dir === :X - dist = X; - elseif dir ===:Y - dist = Y; + dist = X; + elseif dir ===:Y + dist = Y; else - dist = Z; + dist = Z; end - + # compute the 1D thermal structures Temp1 = zeros(size(Temp)); Temp2 = zeros(size(Temp)); @@ -1559,12 +1559,12 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemper ind_1 = findall(weight .>w_max); ind_2 = findall(weight .0 println(io," Weakzone phase : $(g.WeakzonePhase)") end - + return nothing end """ Top, Bot = compute_slab_surface(trench::Trench) -Computes the (`x`,`z`) coordinates of the slab top, bottom surface using the mid surface of the slab as reference. +Computes the (`x`,`z`) coordinates of the slab top, bottom surface using the mid surface of the slab as reference. Parameters -=== +=== - `trench` - `Trench` structure that contains the relevant parameters Method === -It computes it by discretizing the slab surface in `n_seg` segments, and computing the average bending angle (which is a function of the current length of the slab). +It computes it by discretizing the slab surface in `n_seg` segments, and computing the average bending angle (which is a function of the current length of the slab). Next, it compute the coordinates assuming that the trench is at 0.0, and assuming a positive `θ_max` angle. """ function compute_slab_surface(trench::Trench) @@ -1687,7 +1687,7 @@ function compute_slab_surface(trench::Trench) θ = compute_bending_angle(θ_max, Lb, l , type_bending) θ_n = compute_bending_angle(θ_max, Lb, l+dl, type_bending) θ_mean = (θ + θ_n)/2; - + # Mid surface coordinates (x,z) sinθ, cosθ = sincos(θ_mean) @@ -1713,8 +1713,8 @@ function compute_slab_surface(trench::Trench) Top[:,1] *= direction Bottom[:,1] *= direction WeakZone[:,1] *= direction - - return Top, Bottom, WeakZone + + return Top, Bottom, WeakZone end """ @@ -1731,7 +1731,7 @@ Parameters """ function compute_bending_angle(θ_max::Float64,Lb::Float64,l::Float64,type::Symbol) - + if l>Lb θ = θ_max elseif type === :Ribe @@ -1756,12 +1756,12 @@ function find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench::Trench) # Perform rotation of 3D coordinates along the angle from Start -> End: Xrot = X .- Start[1]; Yrot = Y .- Start[2]; - + StrikeAngle = -atand((End[2]-Start[2])/(End[1]-Start[1])) Rot3D!(Xrot,Yrot,Z, StrikeAngle, 0.0) xb = Rot3D(End[1]-Start[1],End[2]-Start[2], 0.0, cosd(StrikeAngle), sind(StrikeAngle), 1.0, 0.0) - + # dl dl = trench.Length/n_seg; l = 0 # length at the trench position @@ -1804,7 +1804,7 @@ function find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench::Trench) d[ip] = -distance_to_linesegment(point_, pa, pd) ls[ip] = distance_to_linesegment(point_, pb, pa) + l end - + #Update l l = ln; end @@ -1826,7 +1826,7 @@ function distance_to_linesegment(p::NTuple{2,_T}, v::NTuple{2,_T}, w::NTuple{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. + # 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 @@ -1848,7 +1848,7 @@ end """ add_slab!(Phase, Temp, Grid::AbstractGeneralGrid, trench::Trench; phase = ConstantPhase(1), T = nothing, cell=false) -Adds a curved slab with phase & temperature structure to a 3D model setup. +Adds a curved slab with phase & temperature structure to a 3D model setup. Parameters ==== @@ -1882,26 +1882,26 @@ julia> add_slab!(Phase, Temp, Cart, trench, phase = phase, T = TsHC) function add_slab!(Phase, Temp, Grid::AbstractGeneralGrid, trench::Trench; # required input phase::AbstractPhaseNumber = ConstantPhase(1), # Sets the phase number(s) in the slab T::Union{AbstractThermalStructure,Nothing} = nothing, cell=false ) # Sets the thermal structure (various functions are available), - + # Retrieve 3D data arrays for the grid X,Y,Z = coordinate_grids(Grid, cell=cell) # Compute top and bottom of the slab - Top,Bottom, WeakZone = compute_slab_surface(trench); - + Top,Bottom, WeakZone = compute_slab_surface(trench); + # Find the distance to the slab (along & perpendicular) - d = fill(NaN,size(Grid)); # -> d = distance perpendicular to the slab + d = fill(NaN,size(Grid)); # -> d = distance perpendicular to the slab ls = fill(NaN,size(Grid)); # -> l = length from the trench along the slab - find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench); + find_slab_distance!(ls, d, X,Y,Z, Top, Bottom, trench); - # Function to fill up the temperature and the phase. + # Function to fill up the temperature and the phase. ind = findall((-trench.Thickness .<= d .<= 0.0)); - + if isa(T, LinearWeightedTemperature) l_decouplingind = findall(Top[:,2].<=-trench.d_decoupling); if !isempty(l_decouplingind) l_decoupling = Top[l_decouplingind[1],1]; - T.crit_dist = abs(l_decoupling); + T.crit_dist = abs(l_decoupling); end end @@ -1915,13 +1915,13 @@ function add_slab!(Phase, Temp, Grid::AbstractGeneralGrid, trench::Trench; # Add a weak zone on top of the slab (indicated by a phase number but not by temperature) if trench.WeakzoneThickness>0.0 - d_weakzone = fill(NaN,size(Grid)); # -> d = distance perpendicular to the slab + d_weakzone = fill(NaN,size(Grid)); # -> d = distance perpendicular to the slab ls_weakzone = fill(NaN,size(Grid)); # -> l = length from the trench along the slab - find_slab_distance!(ls_weakzone, d_weakzone, X,Y,Z, WeakZone, Top, trench); + find_slab_distance!(ls_weakzone, d_weakzone, X,Y,Z, WeakZone, Top, trench); ind = findall( (-trench.WeakzoneThickness .<= d_weakzone .<= 0.0) .& (Z .>-trench.d_decoupling) ); Phase[ind] .= trench.WeakzonePhase end return nothing -end \ No newline at end of file +end