Polygon function compatible with new inPolygon

Author: TatjanaWeiler

docs/make.jl

@@ -82,24 +82,26 @@ makedocs(;
         "Home" => "index.md",
         "Tutorials" => Any[
             "Overview" =>  "man/tutorials.md",
-            "1 - 3D seismic tomography from ASCII" =>  "man/tutorial_load3DSeismicData.md",
-            "2 - 3D seismic tomography from netCDF" =>  "man/tutorial_loadregular3DSeismicData_netCDF.md",
-            "3 - Visualize Moho topography" =>  "man/tutorial_MohoTopo.md",
-            "4 - Create GMT-based topography" =>  "man/tutorial_GMT_Topography.md",
-            "5 - Coastlines" =>  "man/tutorial_Coastlines.md",
-            "6 - Import screenshots" =>  "man/tutorial_Screenshot_To_Paraview.md",
-            "7 - Interpolate irregular 3D seismic tomography" =>  "man/tutorial_loadirregular3DSeismicData.md",
-            "8 - ETOPO1 Topography and geological maps" =>  "man/tutorial_GMT_Topography_GeologicalMap.md",
-            "9 - ISC earthquake data" =>  "man/tutorial_ISC_data.md",
-            "10 - Plot GPS vectors" =>  "man/tutorial_GPS.md",
-            "11 - Read UTM data" =>  "man/tutorial_UTM.md",
-            "12 - VoteMaps" =>  "man/Tutorial_Votemaps.md",
-            "13 - Campi Flegrei" =>  "man/tutorial_local_Flegrei.md",
-            "14 - La Palma volcano Model" =>  "man/Tutorial_LaPalma.md",
-            "15 - Create movies" =>  "man/tutorial_time_Seismicity.md",
-            "16 - Fault Density Map" =>  "man/tutorial_Fault_Map.md",
-            "17 - Jura tutorial" =>  "man/Tutorial_Jura.md",
-            "18 - Geometry setup " => "man/tutorial_Polygon_structures.md"
+            "1 - Getting started" =>  "man/Tutorial_Basic.md",
+            "2 - 3D seismic tomography from ASCII" =>  "man/tutorial_load3DSeismicData.md",
+            "3 - 3D seismic tomography from netCDF" =>  "man/tutorial_loadregular3DSeismicData_netCDF.md",
+            "4 - Visualize Moho topography" =>  "man/Tutorial_MohoTopo_Spada.md",
+            "5 - Create GMT-based topography" =>  "man/tutorial_GMT_Topography.md",
+            "6 - Coastlines" =>  "man/tutorial_Coastlines.md",
+            "7 - Import screenshots" =>  "man/tutorial_Screenshot_To_Paraview.md",
+            "8 - Interpolate irregular 3D seismic tomography" =>  "man/tutorial_loadirregular3DSeismicData.md",
+            "9 - ETOPO1 Topography and geological maps" =>  "man/tutorial_GMT_Topography_GeologicalMap.md",
+            "10 - ISC earthquake data" =>  "man/tutorial_ISC_data.md",
+            "11 - Plot GPS vectors" =>  "man/tutorial_GPS.md",
+            "12 - Read UTM data" =>  "man/tutorial_UTM.md",
+            "13 - VoteMaps" =>  "man/Tutorial_Votemaps.md",
+            "14 - Campi Flegrei" =>  "man/tutorial_local_Flegrei.md",
+            "15 - La Palma volcano Model" =>  "man/Tutorial_LaPalma.md",
+            "16 - Create movies" =>  "man/tutorial_time_Seismicity.md",
+            "17 - Fault Density Map" =>  "man/Tutorial_FaultDensity.md",
+            "18 - Alpine data integration" =>  "man/Tutorial_AlpineData.md",
+            "19 - Jura tutorial" =>  "man/Tutorial_Jura.md", 
+            "20 - Build geometry setup with polygon.md" =>  "man/tutorial_Polygon_structures.md"
         ],
         "User Guide" => Any[
             "Installation" =>  "man/installation.md",

docs/src/man/tutorial_Polygon_structures.md

@@ -8,7 +8,7 @@ This tutorial visualizes simplified geological as it is done here for a passive
 ## Steps
 
 #### 1. Set up your simplified background model
-Before adding specific geological features, a general simplified model setup is necessary. The construction is made by using the AddBox function. For the model the discontinuities are in 15, 45, 145, and 945 km depth.
+Before adding specific geological features, a general simplified model setup is necessary. The construction is made by using the addBox function. For the model the discontinuities are in 15, 45, 145, and 945 km depth.
 
 ```julia
 using GeophysicalModelGenerator
@@ -33,12 +33,14 @@ Temp    = ones(Float64,size(X))*1350
 # add different phases: crust->2, Mantle Lithosphere->3 Mantle->1
 AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(-800.0,0.0), phase = LithosphericPhases(Layers=[15 30 100 800], Phases=[2 3 1 5], Tlab=1300 ), T=LinearTemp(Ttop=20, Tbot=1600))
 
+
+# add air phase 0
 AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(0.0,50.0), phase = ConstantPhase(0), T=ConstantTemp(20.0))
 
  ```
 
 
-#### 2. Add polygone sturture
+#### 2. Add polygon struture
 For including the geological structures of a passive margin into the model via the polygon function, depths down to 150 km are focused on. With the polygon function, it is possible to create different shapes. In the example, the sediment basin shows a more trapezial (2D in x-/z-direction) shape, while the thinning of the plate is a more triangular (2D in x-/z-direction). More complex structures are possible to build in the background model due to the non-limitation number of points in the x- and z-direction. In y-direction only the length can be varied and is set by two values. The shape is not changeable. The x- and z-values of the points need to be in the same order for selecting the correct point (P1(1/3), P2(2/2) --> xlim(1,2), ylim(3,2)).
 
 
@@ -48,18 +50,18 @@ For including the geological structures of a passive margin into the model via t
 # ylim: limits the object within the two ylim values
 # unlimited number of points possible to create the polygon
 
-# add sediment basin # depending on resolution show and angle if it the edge is visible in paraview
-addPolygon!(Phase, Temp, Cart; xlim=(0.0,0.0, 160.0, 200.0),ylim=(100.0,300.0), zlim=(0.0,-10.0,-20.0,0.0), phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+# add sediment basin 
+addPolygon!(Phase, Temp, Cart; xlim=[0.0,0.0, 160.0, 200.0],ylim=[100.0,300.0], zlim=[0.0,-10.0,-20.0,0.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
 
 # add thinning of the continental crust attached to the slab and its thickness 
-addPolygon!(Phase, Temp, Cart; xlim=(0.0, 200.0, 0.0),ylim=(500.0,800.0), zlim=(-100.0,-150.0,-150.0), phase = ConstantPhase(5), T=LinearTemp(Ttop=1000, Tbot=1100))
+addPolygon!(Phase, Temp, Cart; xlim=[0.0, 200.0, 0.0],ylim=[500.0,800.0], zlim=[-100.0,-150.0,-150.0], phase = ConstantPhase(5), T=LinearTemp(Ttop=1000, Tbot=1100))
 
 # add accretionary prism 
-addPolygon!(Phase, Temp, Cart; xlim=(800.1, 600.0, 800.1),ylim=(100.0,800.0), zlim=(0.0,0.0,-60.0), phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+addPolygon!(Phase, Temp, Cart; xlim=[800.0, 600.0, 800.0],ylim=[100.0,800.0], zlim=[0.0,0.0,-60.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
 
  ```
 
-#### 3. Export the final model setup to a paraview file
+#### 3. Export final model setup to a paraview file
 For visualisation and comparison to actual measured data, the mode setup is saved to a paraview file.
 
 ```julia

src/Setup_geometry.jl

@@ -9,7 +9,7 @@ using GeoParams
 # These are routines that help to create input geometries, such as slabs with a given angle
 #
 
-export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addPolygon!
+export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addPolygon!,
         makeVolcTopo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
@@ -472,33 +472,11 @@ end
 
 
 
-
-function inPoly(PolyX, PolyY, x, y)
-    inside = false
-    iSteps = collect(eachindex(PolyX))
-    jSteps = [length(PolyX); collect(1:length(PolyX)-1)]
-
-    for (i,j) in zip(iSteps, jSteps)
-        xi = PolyX[i]
-        yi = PolyY[i]
-        xj = PolyX[j]
-        yj = PolyY[j]
-
-        if ((yi > y) != (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi + eps()) + xi)
-
-            inside = !inside
-        end
-    end
-
-    return inside
-end
-
-
 """
-    addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid;                 # required input
-            xlim=Tuple{}, ylim=Tuple{2}, zlim=Tuple{},     # limits of the box
-            phase = ConstantPhase(1),                       # Sets the phase number(s) in the box
-            T=nothing )    
+    function addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid;    # required input
+        xlim=Vector(), ylim=Vector(2), zlim=Vector(),               # limits of the box
+        phase = ConstantPhase(1),                                   # Sets the phase number(s) in the box
+        T=nothing )   
 Adds a polygon with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
 
 
@@ -510,7 +488,7 @@ Parameters
 - xlim  - x-cooridnate of the polygon points, same ordering as zlim, number of points unlimited
 - ylim  - y-cooridnate, limition in length possible (two values (start and stop))
 - zlim  - z-cooridnate of the polygon points, same ordering as xlim, number of points unlimited
-- phase - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()`
+- phase - specifies the phase of the box. See `ConstantPhase()`
 - T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`
 
 
@@ -540,29 +518,28 @@ julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 """
 
 
-function addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid;                 # required input
-    xlim=Tuple{}, ylim=Tuple{2}, zlim=Tuple{},     # limits of the box
-    phase = ConstantPhase(1),                       # Sets the phase number(s) in the box
-    T=nothing )                                     # Sets the thermal structure (various functions are available)
+function addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid;    # required input
+    xlim=Vector(), ylim=Vector(2), zlim=Vector(),               # limits of the box
+    phase = ConstantPhase(1),                                   # Sets the phase number(s) in the box
+    T=nothing )                                                 # Sets the thermal structure (various functions are available)
 
 # Retrieve 3D data arrays for the grid
 X,Y,Z = coordinate_grids(Grid)
 
+ind = zeros(Bool,size(X))
+ind_slice = zeros(Bool,size(X[:,1,:]))
 
-indx = zeros(length(X))
+# find points within the polygon, only in 2D
+for i = 1:size(Y)[2]
 
-# find points of the total script within the polygone, only in 2D due to the symetric structures and index of y
-for i = 1:length(X)#(indy)
-
-    # working but not the fastest
-    if Y[i] >= ylim[1] && Y[i]<=ylim[2] 
-        indx[i] = inPoly(xlim,zlim, X[i], Z[i])
+    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
+        ind[:,i,:] = zeros(size(X[:,1,:]))
     end
 end
 
-# get all indices which are in the polygone separated
-ind = findall(x->x>0,indx)
-
 
 # Compute thermal structure accordingly. See routines below for different options
 if T != nothing

test/test_setup_geometry.jl

@@ -39,9 +39,6 @@ AddBox!(Phases,Temp,Grid, xlim=(2,4), zlim=(4,8), phase=ConstantPhase(3), DipAng
 
 @test maximum(Phases) == 3
 
-# add accretionary prism 
-addPolygon!(Phase, Temp, Cart; xlim=(1, 2, 1),ylim=(0.0,1.0), zlim=(2.0,3.0,3.0), phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
-@test maximum(Phases) == 8
 
 # Create a CartData structure from it
 Data = CartData(Grid, (T=Temp, Phases=Phases))
@@ -217,3 +214,22 @@ AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0),
 
 
 Data_Final =   AddField(Cart,"Temp",Temp)
+
+# test polygon structure
+
+x        = LinRange(0.0,1200.0,64);
+y        = LinRange(0.0,1200.0,64);
+z        = LinRange(-660,50,64);
+Cart     = CartData(XYZGrid(x, y, z));
+
+# initialize phase and temperature matrix
+Phase   = ones(Int32,(length(x),length(y),length(z)));
+Temp    = ones(Float64,(length(x),length(y),length(z)))*1350;
+
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T=ConstantTemp(120.0));
+
+# add accretionary prism 
+addPolygon!(Phase, Temp, Cart; xlim=[500.0, 200.0, 500.0],ylim=[100.0,400.0], zlim=[0.0,0.0,-60.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+
+@test maximum(Phase) == 8
+@test minimum(Temp) == 21.40845070422536

tutorials/Tutorial_polygon_geometry.jl

@@ -21,23 +21,22 @@ Temp    = ones(Float64,size(X))*1350
 AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(-800.0,0.0), phase = LithosphericPhases(Layers=[15 30 100 800], Phases=[2 3 1 5], Tlab=1300 ), T=LinearTemp(Ttop=20, Tbot=1600) )#T=HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=120, Adiabat=0.4)
 
 
-
-# xlim: x-coordinates of the points, same ordering as zlim
+# xlim: x-coordiates of the points, same ordering as zlim
 # zlim: z-coordinates of the points, same ordering as xlim
 # ylim: limits the object within the two ylim values
 # unlimited number of points possible to create the polygon
-# add sediment basin # depending on resolution show and angle if it the edge is visible in paraview
-addPolygon!(Phase, Temp, Cart; xlim=(0.0,0.0, 160.0, 200.0),ylim=(100.0,300.0), zlim=(0.0,-10.0,-20.0,0.0), phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+# add sediment basin # depending on the resolution and angle if it the edge is visible in paraview
+addPolygon!(Phase, Temp, Cart; xlim=[0.0,0.0, 160.0, 200.0],ylim=[100.0,300.0], zlim=[0.0,-10.0,-20.0,0.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
 
 # add thinning of the continental crust attached to the slab and its thickness 
-addPolygon!(Phase, Temp, Cart; xlim=(0.0, 200.0, 0.0),ylim=(500.0,800.0), zlim=(-100.0,-150.0,-150.0), phase = ConstantPhase(5), T=LinearTemp(Ttop=1000, Tbot=1100))
+addPolygon!(Phase, Temp, Cart; xlim=[0.0, 200.0, 0.0],ylim=[500.0,800.0], zlim=[-100.0,-150.0,-150.0], phase = ConstantPhase(5), T=LinearTemp(Ttop=1000, Tbot=1100))
 
 # add accretionary prism 
-addPolygon!(Phase, Temp, Cart; xlim=(800.1, 600.0, 800.1),ylim=(100.0,800.0), zlim=(0.0,0.0,-60.0), phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+addPolygon!(Phase, Temp, Cart; xlim=[800.0, 600.0, 800.0],ylim=[100.0,800.0], zlim=[0.0,0.0,-60.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
 
 
 # add air phase 0
-addBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(0.0,50.0), phase = ConstantPhase(0), T=ConstantTemp(20.0))
+AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(0.0,50.0), phase = ConstantPhase(0), T=ConstantTemp(20.0))
 
 # # Save data to paraview:
 Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp))