diff --git a/.typos.toml b/.typos.toml
index 353f06a0..5e69ca29 100644
--- a/.typos.toml
+++ b/.typos.toml
@@ -2,6 +2,8 @@
 MOR = "MOR"
 dum = "dum"
 Shepard = "Shepard"
+arange = "arange"
+iy = "iy"
 nin = "nin"
 
 [files]
diff --git a/docs/src/man/Tutorial_NumericalModel_2D.md b/docs/src/man/Tutorial_NumericalModel_2D.md
index 2789312b..6e767171 100644
--- a/docs/src/man/Tutorial_NumericalModel_2D.md
+++ b/docs/src/man/Tutorial_NumericalModel_2D.md
@@ -47,13 +47,13 @@ Temp = fill(1350.0, nx, 1, nz);
 We will start with a simple subduction setup, which consists of a horizontal part:
 
 ```julia
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1));
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1));
 ```
 
 And with the inclined part:
 
 ```julia
-add_box!(Phases, Temp, Grid2D; xlim=(0,300), zlim=(-80.0, 0.0), phase = ConstantPhase(1), DipAngle=30);
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,300.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1), DipAngle=30);
 ```
 
 Add them to the `CartData` dataset:
@@ -100,8 +100,8 @@ LithosphericPhases([15 55], [1 2], nothing)
 and set the slab again:
 
 ```julia
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = lith);
-add_box!(Phases, Temp, Grid2D; xlim=(0,300), zlim=(-80.0, 0.0), phase = lith, DipAngle=30);
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = lith);
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,300.0), zlim=(-80.0, 0.0), phase = lith, DipAngle=30);
 ```
 
 Which looks like:
@@ -124,8 +124,8 @@ We can do that by specifying a thermal structure. For example, we can use the ha
 
 ```julia
 therm = HalfspaceCoolingTemp(Age=40)
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = lith, T=therm);
-add_box!(Phases, Temp, Grid2D; xlim=(0,300), zlim=(-80.0, 0.0), phase = lith, T = therm, DipAngle=30);
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = lith, T=therm);
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,300.0), zlim=(-80.0, 0.0), phase = lith, T = therm, DipAngle=30);
 ```
 
 Which looks like:
@@ -159,13 +159,13 @@ a spreading velocity (note that this simply relates to the thermal structure and
 
 ```julia
 lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
 ```
 
 For the subduction we use a thermal structure of a slab heated by hot asthenosphere
 
 ```julia
-add_box!(Phases, Temp, Grid2D; xlim=(0,300), zlim=(-80.0, 0.0), phase = lith, T = McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=30);
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,300.0), zlim=(-80.0, 0.0), phase = lith, T = McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=30);
 ```
 
 We can set the mantle lithosphere that is hotter > 1250 C to mantle:
@@ -186,25 +186,25 @@ Saved file: Grid2D_SubductionRidge.vts
 ![Mechanical2D_Tutorial_4](../assets/img/Mechanical2D_Tutorial_4.png)
 
 #### Overriding slab and weak layer
-Ok, lets add an overriding slab as well. For this, we use the `AddLayer!` function
+Ok, lets add an overriding slab as well. For this, we use the `add_layer!` function
 
 ```julia
 lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)
-add_box!(Phases, Temp, Grid2D; xlim=(0,1000), zlim=(-80.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,1000.0), zlim=(-80.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
 ```
 
 The oceanic plate is as before
 
 ```julia
 lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
 ```
 
 For the inclined part, we set a layer above the slab (the "weak" layer to facilitate subduction initiation )
 
 ```julia
 lith = LithosphericPhases(Layers=[10 15 55], Phases=[6 1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid2D; xlim=(0,300), zlim=(-80.0, 10.0), phase = lith, T = McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=30);
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,300.0), zlim=(-80.0, 10.0), phase = lith, T = McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=30);
 ```
 
 Lithosphere-asthenosphere boundary:
@@ -236,7 +236,7 @@ z = range(-660,0,    nz);
 Grid2D = CartData(xyz_grid(x,0,z))
 Phases = zeros(Int64, nx, 1, nz);
 Temp = fill(1350.0, nx, 1, nz);
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1));
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1));
 ```
 
 Next, we should define a `Trench` structure, which contains info about the trench which goes in 3D from `Start` - `End` coordinates (`x`,`y`)-coordinates respectively. As we are dealing with a 2D model, we set the `y`-coordinates to -100.0 and 100.0 respectively.
@@ -284,8 +284,8 @@ LithosphericPhases([15 20 55], [3 4 5], 1250)
 Lets start with defining the horizontal part of the overriding plate. Note that we define this twice with different thickness to deal with the bending subduction area:
 
 ```julia
-add_box!(Phases, Temp, Grid2D; xlim=(200,1000), zlim=(-150.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
-add_box!(Phases, Temp, Grid2D; xlim=(0,200), zlim=(-50.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+add_box!(Phases, Temp, Grid2D; xlim=(200.0,1000.0), zlim=(-150.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,200.0), zlim=(-50.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
 ```
 
 The horizontal part of the oceanic plate is as before:
@@ -293,7 +293,7 @@ The horizontal part of the oceanic plate is as before:
 ```julia
 v_spread_cm_yr = 3      #spreading velocity
 lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-150.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=v_spread_cm_yr));
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-150.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=v_spread_cm_yr));
 ```
 
 Yet, now we add a trench as well. The starting thermal age at the trench is that of the horizontal part of the oceanic plate:
diff --git a/docs/src/man/Tutorial_NumericalModel_3D.md b/docs/src/man/Tutorial_NumericalModel_3D.md
index 01531e1d..4e26e292 100644
--- a/docs/src/man/Tutorial_NumericalModel_3D.md
+++ b/docs/src/man/Tutorial_NumericalModel_3D.md
@@ -43,7 +43,7 @@ Note that if the lowermost layer has the same phase as the mantle, you can defin
 
 ```julia
 lith = LithosphericPhases(Layers=[15 45 10], Phases=[0 1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid; xlim=(-800,0.0), ylim=(-400, 400.0), zlim=(-80.0, 0.0), phase = lith,
+add_box!(Phases, Temp, Grid; xlim=(-800.0,0.0), ylim=(-400, 400.0), zlim=(-80.0, 0.0), phase = lith,
         Origin=(-0,0,0),
         T=SpreadingRateTemp(SpreadingVel=3, MORside="right"), StrikeAngle=30);
 ```
@@ -51,7 +51,7 @@ add_box!(Phases, Temp, Grid; xlim=(-800,0.0), ylim=(-400, 400.0), zlim=(-80.0, 0
 And an an inclined part:
 
 ```julia
-add_box!(Phases, Temp, Grid; xlim=(0,300), ylim=(-400, 400.0), zlim=(-80.0, 0.0), phase = lith,
+add_box!(Phases, Temp, Grid; xlim=(0.0,300.0), ylim=(-400.0, 400.0), zlim=(-80.0, 0.0), phase = lith,
         Origin=(-0,0,0),
         T=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=30, StrikeAngle=30);
 ```
@@ -104,12 +104,12 @@ Overriding plate with a 30 km crust and mantle lithosphere that where T<1250 cel
 
 ```julia
 lith_cont = LithosphericPhases(Layers=[30 200 50], Phases=[3 4 2], Tlab=1250)
-add_box!(Phases, Temp, Grid; xlim=(400,1000), ylim=(-1000, 0.0), zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=150));
-add_box!(Phases, Temp, Grid; xlim=(200,1000), ylim=(-1000, 0.0), zlim=(-80.0, 0.0), phase = lith_cont,  T=HalfspaceCoolingTemp(Age=150));
+add_box!(Phases, Temp, Grid; xlim=(400.0,1000.0), ylim=(-1000.0, 0.0), zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=150));
+add_box!(Phases, Temp, Grid; xlim=(200.0,1000.0), ylim=(-1000.0, 0.0), zlim=(-80.0, 0.0), phase = lith_cont,  T=HalfspaceCoolingTemp(Age=150));
 
 lith_cont = LithosphericPhases(Layers=[30 200 10], Phases=[5 6 2], Tlab=1250)
-add_box!(Phases, Temp, Grid; xlim=(400,1000), ylim=(0, 1000),    zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
-add_box!(Phases, Temp, Grid; xlim=(200,1000), ylim=(0, 1000),    zlim=( -80.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
+add_box!(Phases, Temp, Grid; xlim=(400.0,1000.0), ylim=(0.0, 1000.0),    zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
+add_box!(Phases, Temp, Grid; xlim=(200.0,1000.0), ylim=(0.0, 1000.0),    zlim=( -80.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
 ```
 
 Define an oceanic plate with ridge
@@ -117,14 +117,14 @@ Define an oceanic plate with ridge
 ```julia
 v_spread_cm_yr = 3      #spreading velocity
 lith = LithosphericPhases(Layers=[15 45 10], Phases=[0 1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid; xlim=(-800 , 200), ylim=(-1000, -400.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
-add_box!(Phases, Temp, Grid; xlim=(-1000,-800), ylim=(-1000, -400.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
+add_box!(Phases, Temp, Grid; xlim=(-800.0 , 200.0), ylim=(-1000.0, -400.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
+add_box!(Phases, Temp, Grid; xlim=(-1000.0,-800.0), ylim=(-1000.0, -400.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
 
-add_box!(Phases, Temp, Grid; xlim=(-700,  200), ylim=(-400, 200.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
-add_box!(Phases, Temp, Grid; xlim=(-1000,-700), ylim=(-400, 200.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
+add_box!(Phases, Temp, Grid; xlim=(-700.0,  200.0), ylim=(-400.0, 200.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
+add_box!(Phases, Temp, Grid; xlim=(-1000.0,-700.0), ylim=(-400.0, 200.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
 
-add_box!(Phases, Temp, Grid; xlim=(-650,  200), ylim=(200, 1000.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
-add_box!(Phases, Temp, Grid; xlim=(-1000,-650), ylim=(200, 1000.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
+add_box!(Phases, Temp, Grid; xlim=(-650.0,  200.0), ylim=(200.0, 1000.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
+add_box!(Phases, Temp, Grid; xlim=(-1000.0,-650.0), ylim=(200.0, 1000.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
 ```
 
 Subducting parts of the oceanic plate
@@ -153,7 +153,7 @@ T_slab  = LinearWeightedTemperature( F1=HalfspaceCoolingTemp(Age=AgeTrench_Myrs)
 add_slab!(Phases, Temp, Grid, trench1, phase = lith, T=T_slab);
 ```
 
-Finally, it is often nice to see the deformation of the plate when it subducts. A simple way to do that is to put a `stripes` on top using `addStripes`, which has the same phase as the subducting crust.
+Finally, it is often nice to see the deformation of the plate when it subducts. A simple way to do that is to put a `stripes` on top using `add_stripes!`, which has the same phase as the subducting crust.
 
 ```julia
 add_stripes!(Phases, Grid; stripAxes = (1,1,0), phase = ConstantPhase(0), stripePhase = ConstantPhase(9), stripeWidth=50, stripeSpacing=200)
diff --git a/docs/src/man/tutorial_Polygon_structures.md b/docs/src/man/tutorial_Polygon_structures.md
index 1b4ec706..81bc6ff5 100644
--- a/docs/src/man/tutorial_Polygon_structures.md
+++ b/docs/src/man/tutorial_Polygon_structures.md
@@ -25,8 +25,8 @@ z        = LinRange(-660,50,nz)
 Cart     = CartData(xyz_grid(x, y, z))
 
 # initialize phase and temperature matrix
-Phase   = ones(Int32,size(X))
-Temp    = ones(Float64,size(X))*1350
+Phase = fill(1,nx,ny,nz);
+Temp = fill(1350.0, nx,ny,nz);
 
 # add different phases: crust->2, Mantle Lithosphere->3 Mantle->1
 add_box!(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))
@@ -47,13 +47,13 @@ To include the geological structures of a passive margin into the model, we use
 # unlimited number of points possible to create the polygon
 
 # add sediment basin 
-add_polygon!(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_polygon!(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 
-add_polygon!(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_polygon!(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 
-add_polygon!(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_polygon!(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 final model setup to a paraview file
@@ -61,8 +61,8 @@ For visualisation and comparison to actual measured data, the mode setup is save
 
 ```julia
 # # Save data to paraview:
-Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp)); 
-write_paraview(Data_Final, "Sedimentary_basin");
+Cart = addfield(Cart,(;Phase, Temp))
+write_paraview(Cart, "Sedimentary_basin");
 ```
 
 After importing and looking at the file to paraview, some unresolved areas might be visible as they are visible in this model. That is due to the resolution and shape of the polygon. To reduce those artefacts an increase in resolution or a change of the polygon angle might help.
diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index ea7e420c..64d736f0 100644
--- a/src/Setup_geometry.jl
+++ b/src/Setup_geometry.jl
@@ -130,11 +130,11 @@ end
 
 
 """
-    add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
+    add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::NTuple{2, _T} = (20,100), [ylim::NTuple{2, _T} = (1,10)], zlim::NTuple{2, _T} = (10,80),
             Origin=nothing, StrikeAngle=0, DipAngle=0,
             phase = ConstantPhase(1),
             T=nothing,
-            cell=false )
+            cell=false ) where _T
 
 Adds a box with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
 
@@ -189,11 +189,11 @@ julia> write_paraview(Grid,"LaMEM_ModelSetup")  # Save model to paraview
 ```
 """
 function add_box!(Phase, Temp, Grid::AbstractGeneralGrid;       # required input
-                xlim=Tuple{2}, ylim=nothing, zlim=Tuple{2},     # limits of the box
+                xlim::NTuple{2, _T} = (20,100), ylim=nothing, zlim::NTuple{2, _T} = (10,80),     # 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)
-                cell=false )                                    # if true, Phase and Temp are defined on cell centers 
+                cell=false ) where _T                           # 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)
@@ -247,9 +247,10 @@ end
 
 
 """
-    add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
+    add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::NTuple{2, _T} = (1,100), [ylim::NTuple{2, _T} = (0,20)], zlim::NTuple{2, _T} = (0,-100),
             phase = ConstantPhase(1),
-            T=nothing, cell=false )
+            T=nothing, cell=false ) where _T
+
 
 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.
 This simplifies creating model geometries in geodynamic models
@@ -348,9 +349,10 @@ end
 
 
 """
-    add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen=Tuple{3}, radius=Tuple{1},
+    add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen::NTuple{3, _T} = (0,0,-1), radius::Number,
             phase = ConstantPhase(1).
-            T=nothing, cell=false )
+            T=nothing, cell=false ) where _T
+
 
 Adds a sphere with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
 
@@ -390,9 +392,9 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 ```
 """
 function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid;    # required input
-    cen=Tuple{3}, radius=Tuple{1},                              # center and radius of the sphere
+    cen::NTuple{3, _T} = (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 )                                     # Sets the thermal structure (various functions are available)
+    T=nothing, cell=false )    where _T                         # Sets the thermal structure (various functions are available)
 
     # Retrieve 3D data arrays for the grid
     X,Y,Z = coordinate_grids(Grid, cell=cell)
@@ -412,10 +414,10 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid;    # required input
 end
 
 """
-    add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen=Tuple{3}, axes=Tuple{3},
+    add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen::NTuple{3, _T} = (-1,-1,-1), axes::NTuple{3, _T} = (0.2,0.1,0.5),
             Origin=nothing, StrikeAngle=0, DipAngle=0,
             phase = ConstantPhase(1).
-            T=nothing, cell=false )
+            T=nothing, cell=false ) where _T
 
 Adds an Ellipsoid with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
 
@@ -457,10 +459,10 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 ```
 """
 function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid;     # required input
-    cen=Tuple{3}, axes=Tuple{3},                                    # center and semi-axes of the ellpsoid
+    cen::NTuple{3, _T} = (-1,-1,-1), axes::NTuple{3, _T} = (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 )                                         # Sets the thermal structure (various functions are available)
+    T=nothing, cell=false )  where _T                               # Sets the thermal structure (various functions are available)
 
     if Origin==nothing
         Origin = cen  # center
@@ -496,9 +498,10 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid;     # required i
 end
 
 """
-    add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},
-            phase = ConstantPhase(1).
-            T=nothing, cell=false )
+    add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::NTuple{3, _T} = (-1,-1,-1.5), cap::NTuple{3, _T} = (-1,-1,-0.5), radius::Number,
+            phase = ConstantPhase(1),
+            T=nothing, cell=false ) where _T
+
 
 Adds a cylinder with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
 
@@ -539,9 +542,9 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 ```
 """
 function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid;  # required input
-    base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},               # center and radius of the sphere
+    base::NTuple{3, _T} = (-1,-1,-1.5), cap::NTuple{3, _T} = (-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)
+    T=nothing, cell=false )   where _T                          # Sets the thermal structure (various functions are available)
 
     # axis vector of cylinder
     axVec = cap .- base
@@ -593,7 +596,8 @@ end
 
 
 """
-        add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Vector(), ylim=Vector(2), zlim=Vector(), phase = ConstantPhase(1), T=nothing, cell=false )   
+        add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::NTuple{2, _T} = (0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false )   where _T
+
 
 Adds a polygon with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
 
@@ -625,7 +629,7 @@ LaMEM Grid:
   z           ϵ [-2.0 : 0.0]
 julia> Phases = zeros(Int32,   size(Grid.X));
 julia> Temp   = zeros(Float64, size(Grid.X));
-julia> add_polygon!(Phase, Temp, Cart; xlim=[0.0,0.0, 1.6, 2.0],ylim=[0.0,0.8], zlim=[0.0,-1.0,-2.0,0.0], phase = ConstantPhase(8), T=ConstantTemp(30))
+julia> add_polygon!(Phase, Temp, Cart; xlim=(0,0, 1.6, 2.0),ylim=(0,0.8), zlim=(0,-1,-2,0), phase = ConstantPhase(8), T=ConstantTemp(30))
 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}:
@@ -634,9 +638,15 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 
 """
 function add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid;   # required input
-    xlim::Vector=[], ylim::Vector=[], zlim::Vector=[],          # limits of the box
+    xlim=(), ylim::NTuple{2, _T} = (0,0.8), zlim=(),           # limits of the box
     phase = ConstantPhase(1),                                   # Sets the phase number(s) in the box
-    T=nothing, cell=false )                                     # Sets the thermal structure (various functions are available)
+    T=nothing, cell=false )  where _T                           # Sets the thermal structure (various functions are available)
+
+
+    xlim_ = Float64.(collect(xlim))
+    ylim_ = Float64.(collect(ylim))
+    zlim_ = Float64.(collect(zlim))
+
 
 # Retrieve 3D data arrays for the grid
 X,Y,Z = coordinate_grids(Grid, cell=cell)
@@ -646,8 +656,8 @@ 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] 
-        inpolygon!(ind_slice, 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,:]))
diff --git a/test/test_lamem.jl b/test/test_lamem.jl
index 10de4f1d..a022938d 100644
--- a/test/test_lamem.jl
+++ b/test/test_lamem.jl
@@ -93,7 +93,7 @@ rm("test_topo.dat")
 Grid    = read_LaMEM_inputfile("test_files/GeometricPrimitives.dat")
 Phases  = zeros(Int32,size(Grid.X));
 Temp    = zeros(Float64,size(Grid.X));
-add_sphere!(Phases,Temp,Grid, cen=(0,0,-6), radius=2, phase=ConstantPhase(1), T=ConstantTemp(800))
+add_sphere!(Phases,Temp,Grid, cen=(0,0,-6), radius=2.0, phase=ConstantPhase(1), T=ConstantTemp(800))
 @test Phases[55,55,55] == 1
 @test Phases[56,56,56] == 0
 @test Temp[44,52,21]   == 800.0
diff --git a/test/test_setup_geometry.jl b/test/test_setup_geometry.jl
index f91de7bf..7069e2ad 100644
--- a/test/test_setup_geometry.jl
+++ b/test/test_setup_geometry.jl
@@ -195,7 +195,7 @@ add_box!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0)
 
 # inclined slab
 Temp    = ones(Float64,size(Cart))*1350;
-add_box!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0,0),StrikeAngle=0, DipAngle=45, phase = ConstantPhase(5), T=TsMK);
+add_box!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0,0.0),StrikeAngle=0, DipAngle=45, phase = ConstantPhase(5), T=TsMK);
 @test sum(Temp)  ≈ 3.5125017626287365e8
 
 
@@ -230,7 +230,7 @@ Temp    = ones(Float64,(length(x),length(y),length(z)))*1350;
 add_box!(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 
-add_polygon!(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))
+add_polygon!(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
@@ -284,6 +284,7 @@ TsHC = HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=120, Adiabat=0.4)
 TsMK = McKenzie_subducting_slab(Tsurface = 20.0, Tmantle = 1350.0, v_cm_yr = 4.0, Adiabat = 0.0)
 T_slab = LinearWeightedTemperature(crit_dist=600, F1=TsHC, F2=TsMK);
 phase = LithosphericPhases(Layers=[5 7 88], Phases = [2 3 4], Tlab=nothing)
+
 t1 = Trench(Start = (400.0,400.0), End = (800.0,800.0),θ_max = 90.0, direction = 1.0, n_seg = 50, Length = 600.0, Thickness = 80.0, Lb = 500.0,d_decoupling = 100.0, type_bending =:Ribe,   WeakzoneThickness=10, WeakzonePhase=9)
 
 add_slab!(Phase,Temp,Cart, t1, phase=phase, T = T_slab)
@@ -298,7 +299,7 @@ z = range(-660,0,    nz);
 Grid2D = CartData(xyz_grid(x,0,z))
 Phases = zeros(Int64, nx, 1, nz);
 Temp = fill(1350.0, nx, 1, nz);
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1),  T=HalfspaceCoolingTemp(Age=40));    
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1),  T=HalfspaceCoolingTemp(Age=40));    
 
 trench = Trench(Start=(0.0,-100.0), End=(0.0,100.0), Thickness=80.0, θ_max=30.0, Length=300, Lb=150);
 add_slab!(Phases, Temp, Grid2D, trench, phase = ConstantPhase(2), T=HalfspaceCoolingTemp(Age=40));
@@ -324,13 +325,13 @@ Temp = fill(1350.0, nx, 1, nz);
 lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)
 
 # Lets add the overriding plate. Note that we add this twice with a different thickness to properly represent the transition around the trench
-add_box!(Phases, Temp, Grid2D; xlim=(200,1000), zlim=(-150.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
-add_box!(Phases, Temp, Grid2D; xlim=(0,200), zlim=(-60.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+add_box!(Phases, Temp, Grid2D; xlim=(200.0,1000.0), zlim=(-150.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+add_box!(Phases, Temp, Grid2D; xlim=(0.0,200.0), zlim=(-60.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
 
 # The horizontal part of the oceanic plate is as before
 v_spread_cm_yr = 3      #spreading velocity
 lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
-add_box!(Phases, Temp, Grid2D; xlim=(-800,0.0), zlim=(-150.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=v_spread_cm_yr));
+add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-150.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=v_spread_cm_yr));
 
 # Yet, now we add a trench as well. 
 AgeTrench_Myrs = 800e3/(v_spread_cm_yr/1e2)/1e6    #plate age @ trench
diff --git a/tutorials/Tutorial_polygon_geometry.jl b/tutorials/Tutorial_polygon_geometry.jl
index 6638727c..3036d5d6 100644
--- a/tutorials/Tutorial_polygon_geometry.jl
+++ b/tutorials/Tutorial_polygon_geometry.jl
@@ -15,7 +15,7 @@ Cart     = CartData(X,Y,Z, (Data=Z,))
 
 # initialize phase and temperature matrix
 Phase   = ones(Int32,size(X))
-Temp    = ones(Float64,size(X))*1350
+Temp    = ones(Float64,size(X)) * 1350
 
 # add different phases: crust->2, Mantle Lithosphere->3 Mantle->1
 add_box!(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)