diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index aa7fa64b..c5f1f79a 100644
--- a/src/Setup_geometry.jl
+++ b/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!,
+export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, AddPolygon!,
         makeVolcTopo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
@@ -450,6 +450,65 @@ function AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid;   # required input
     return nothing
 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
+
+
+function AddPolygon!(Phase, Temp, Grid::AbstractGeneralGrid;                 # required input
+    xlim=Tuple{}, ylim=Tuple{2}, zlim=Tuple{},     # limits of the box
+    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)
+
+# Retrieve 3D data arrays for the grid
+X,Y,Z = coordinate_grids(Grid)
+
+
+indx = zeros(length(X))
+
+# 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])
+    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
+Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], T)
+end
+
+# Set the phase. Different routines are available for that - see below.
+Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
+
+return nothing
+end
+
+
+
 # Internal function that rotates the coordinates
 function Rot3D!(X,Y,Z, StrikeAngle, DipAngle)
 
diff --git a/tutorials/Tutorial_polygon_geometry b/tutorials/Tutorial_polygon_geometry
new file mode 100644
index 00000000..f7b10de9
--- /dev/null
+++ b/tutorials/Tutorial_polygon_geometry
@@ -0,0 +1,46 @@
+using GeophysicalModelGenerator
+
+
+# number of cells in every direction
+nx = 100
+ny = 100
+nz = 200
+
+# define domain size
+x        = LinRange(0.0,800.0,nx)
+y        = LinRange(0.0,800.0,ny)
+z        = LinRange(-660,50,nz)
+X,Y,Z    = XYZGrid(x, y, z);
+Cart     = CartData(X,Y,Z, (Data=Z,))
+
+# initialize phase and temperature matrix
+Phase   = ones(Int32,size(X))
+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) )#T=HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=120, Adiabat=0.4)
+
+
+
+# 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 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))
+
+# 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))
+
+
+# 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))
+
+# # Save data to paraview:
+Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp)) 
+Write_Paraview(Data_Final, "Sedimentary_basin")
+
+