diff --git a/src/LaMEM.jl b/src/LaMEM.jl
index 08304f93..949af667 100644
--- a/src/LaMEM.jl
+++ b/src/LaMEM.jl
@@ -28,13 +28,13 @@ using .LaMEM_Model
 export  LaMEM_Model, Model, Write_LaMEM_InputFile, create_initialsetup,
         Scaling, Grid, Time, FreeSurface, BoundaryConditions, VelocityBox, SolutionParams,
         Solver, ModelSetup, 
-        geom_Sphere,
+        geom_Sphere, geom_Ellipsoid, geom_Box, geom_RidgeSeg, geom_Hex, geom_Layer, geom_Cylinder,
         Output,
         Multigrid, print_short, 
         Materials, Phase, Softening, PhaseTransition, PhaseAggregate, Dike, PassiveTracers,
         add_vbox!, rm_vbox!, rm_last_vbox!, 
         add_phase!, rm_phase!, rm_last_phase!, replace_phase!, add_petsc!, add_softening!, add_phaseaggregate!,
-        add_phasetransition!, add_dike!, add_geom!, set_air, copy_phase,
+        add_phasetransition!, add_dike!, add_geom!, rm_geom!, set_air, copy_phase,
         add_topography!, AboveSurface!, BelowSurface!,
         prepare_lamem, isdefault, hasplasticity,
         add_geoparams_rheologies,
diff --git a/src/LaMEM_ModelGeneration/Model.jl b/src/LaMEM_ModelGeneration/Model.jl
index f906765c..4a9b079d 100644
--- a/src/LaMEM_ModelGeneration/Model.jl
+++ b/src/LaMEM_ModelGeneration/Model.jl
@@ -216,7 +216,7 @@ end
 Prepares a LaMEM run for the parameters that are specified in `model`, without running the simulation
 
     1) Create the `*.dat` file
-    2) Write markers to disk
+    2) Write markers to disk in case we use a "files" setup
 
 This is useful if you want to prepare a model on one machine but run it on another one (e.g. a cluster)
 
@@ -265,6 +265,9 @@ end
 
 Creates the initial model setup of LaMEM from `model`, which includes:
 - Writing the LaMEM (*.dat) input file
+
+and in case we do not employt geometric primitives to create the setup:
+
 - Write the VTK file (if requested when `model.Output.write_VTK_setup=true`)
 - Write the marker files to disk (if `model.ModelSetup.msetup="files"`)
 
@@ -280,6 +283,7 @@ function create_initialsetup(model::Model, cores::Int64=1, args::String=""; verb
 
     Write_LaMEM_InputFile(model, model.Output.param_file_name)
     
+    
     if !isnothing(model.FreeSurface.Topography)
         Save_LaMEMTopography(model.FreeSurface.Topography, model.FreeSurface.surf_topo_file)
     end
diff --git a/src/LaMEM_ModelGeneration/ModelSetup.jl b/src/LaMEM_ModelGeneration/ModelSetup.jl
index 5ba538a9..3a05984f 100644
--- a/src/LaMEM_ModelGeneration/ModelSetup.jl
+++ b/src/LaMEM_ModelGeneration/ModelSetup.jl
@@ -1,6 +1,8 @@
 # Model Setup
 
-export ModelSetup, Write_LaMEM_InputFile, geom_Sphere, set_geom!
+export ModelSetup, Write_LaMEM_InputFile, 
+        geom_Sphere, geom_Ellipsoid, geom_Box, geom_RidgeSeg, geom_Hex, geom_Layer, geom_Cylinder,
+        set_geom!
 
 """
     Structure that contains the LaMEM Model Setup and Advection options
@@ -8,7 +10,7 @@ export ModelSetup, Write_LaMEM_InputFile, geom_Sphere, set_geom!
     $(TYPEDFIELDS)
 """
 Base.@kwdef mutable struct ModelSetup
-    "Setup type - can be `geom` (phases are assigned from geometric primitives), `files` (from julia input), `polygons` (from geomIO input, which requires `poly_file` to be specified) "
+    "Setup type - can be `geom` (phases are assigned from geometric primitives, using `add_geom!(model, ...)`), `files` (from julia input), `polygons` (from geomIO input, which requires `poly_file` to be specified) "
     msetup::String          = "files"              
      
     "add random noise to the particle location"
@@ -91,13 +93,20 @@ end
 
 
 function show_short(io::IO, d::ModelSetup)
-    println(io,"|-- Model setup options :  Type=$(d.msetup); ")
+    str_geom = ""
+    if d.msetup == "geom"
+        str_geom = "$(length(d.geom_primitives)) geometric primitive objects"
+    end
+
+    println(io,"|-- Model setup options :  Type=$(d.msetup); $(str_geom)")
 
     return nothing
 end
 
 
 # Geometric primitives======
+
+# -------
 """
     LaMEM geometric primitive `sphere` object 
     
@@ -122,7 +131,7 @@ end
 
 
 function show(io::IO, d::geom_Sphere)
-    println(io, "Sphere(ph=$(d.phase), r=$(d.radius), cen=$(d.center), T=$(d.Temperature)=$(d.cstTemp))")
+    println(io, "Sphere(ph=$(d.phase), radius=$(d.radius), center=$(d.center), Temperature=$(d.Temperature),  cstTemp = $(d.cstTemp))")
     return nothing
 end
 
@@ -144,6 +153,328 @@ function Write_LaMEM_InputFile(io, d::geom_Sphere)
     println(io, "    <SphereEnd>")
     return nothing
 end
+# -------
+
+# -------
+"""
+    LaMEM geometric primitive `Ellipsoid` object 
+    
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef struct geom_Ellipsoid
+    "phase"
+    phase::Int64          = 1          
+     
+    "semi-axes of ellipsoid in `x`, `y` and `z` "
+    axes::Vector{Float64}      = [2.0, 1.5, 1.0]
+      
+    "center of sphere"
+    center::Vector{Float64} = [1.0, 2.0, 3.0]
+        
+    "optional: Temperature of the sphere. possibilities: [constant, or nothing]"
+    Temperature::Union{String,Nothing} = nothing
+        
+    "required in case of [constant]: temperature value [in Celcius in case of GEO units]"
+    cstTemp::Union{Float64,Nothing}     = nothing   
+end
+
+
+function show(io::IO, d::geom_Ellipsoid)
+    println(io, "Ellipsoid(ph=$(d.phase), axes=$(d.axes), cen=$(d.center), T=$(d.Temperature)=$(d.cstTemp))")
+    return nothing
+end
+
+"""
+    Write_LaMEM_InputFile(io, d::geom_Ellipsoid)
+"""
+function Write_LaMEM_InputFile(io, d::geom_Ellipsoid)
+    fields    = fieldnames(typeof(d))
+    println(io, "    <EllipsoidStart>")
+    for f in fields
+        if !isnothing(getfield(d,f))
+            name = rpad(String(f),15)
+            comment = get_doc(geom_Ellipsoid, f)
+            data = getfield(d,f) 
+            println(io,"        $name  = $(write_vec(data))     # $(comment)")
+        end
+    end
+    println(io, "    <EllipsoidEnd>")
+    return nothing
+end
+# -------
+
+# -------
+"""
+    LaMEM geometric primitive `Box` object 
+    
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef struct geom_Box
+    "phase"
+    phase::Int64          = 1          
+     
+    "box bound coordinates: `left`, `right`, `front`, `back`, `bottom`, `top` "
+    bounds::Vector{Float64}      = [1.0, 2.0, 1.0, 2.0, 1.0, 2.0]
+        
+    "optional: Temperature structure. possibilities: [constant, linear, halfspace]"
+    Temperature::Union{String,Nothing} = nothing
+        
+    "required in case of [`constant`]: temperature value [in Celcius in case of GEO units]"
+    cstTemp::Union{Float64,Nothing}     = nothing   
+
+    "required in case of [`linear,halfspace`]: temperature @ top [in Celcius in case of GEO units]"
+    topTemp::Union{Float64,Nothing}     = nothing   
+    
+    "required in case of [`linear,halfspace`]: temperature @ top [in Celcius in case of GEO units]"
+    botTemp::Union{Float64,Nothing}     = nothing   
+    
+    "required in case of [`halfspace`]: thermal age of lithosphere [in Myrs if GEO units are used]"
+    thermalAge::Union{Float64,Nothing}     = nothing   
+end
+
+function show(io::IO, d::geom_Box)
+    println(io, "Box(ph=$(d.phase), bounds=$(d.bounds), T=$(d.Temperature)=$(d.cstTemp), [top,box]=[$(d.topTemp), $(d.botTemp)], thermalAge=$(d.thermalAge))")
+    return nothing
+end
+
+"""
+    Write_LaMEM_InputFile(io, d::geom_Box)
+"""
+function Write_LaMEM_InputFile(io, d::geom_Box)
+    fields    = fieldnames(typeof(d))
+    println(io, "    <BoxStart>")
+    for f in fields
+        if !isnothing(getfield(d,f))
+            name = rpad(String(f),15)
+            comment = get_doc(geom_Box, f)
+            data = getfield(d,f) 
+            println(io,"        $name  = $(write_vec(data))     # $(comment)")
+        end
+    end
+    println(io, "    <BoxEnd>")
+    return nothing
+end
+# -------
+
+# -------
+"""
+    LaMEM geometric primitive `RidgeSeg` object 
+    
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef struct geom_RidgeSeg
+    "phase"
+    phase::Int64          = 1          
+     
+    "box bound coordinates: `left`, `right`, `front`, `back`, `bottom`, `top` "
+    bounds::Vector{Float64}      = [1.0, 2.0, 1.0, 2.0, 1.0, 2.0]
+        
+    "coordinate order: left, right [can be different for oblique ridge]"
+    ridgeseg_x::Vector{Float64}      = [1.5, 1.5]
+    
+    "coordinate order: front, back [can be different for oblique ridge]"
+    ridgeseg_y::Vector{Float64}      = [1.0, 2.0]
+    
+    "initial temperature structure [ridge must be set to `halfspace_age` --> setTemp=4]"
+    Temperature::String = "halfspace_age"
+
+    "required in case of [`linear,halfspace`]: temperature @ top [in Celcius in case of GEO units]"
+    topTemp::Float64     = 0.0   
+    
+    "required in case of [`linear,halfspace`]: temperature @ top [in Celcius in case of GEO units]"
+    botTemp::Float64    = 1300.0   
+    
+    "minimum age of seafloor at ridge [in `Myr` in case of GEO units]"
+    age0::Float64     = 0.001   
+
+    "[optional] parameter that indicates the maximum thermal age of a plate "
+    maxAge::Union{Float64,Nothing}      = nothing
+
+    "[optional] parameter that indicates the spreading velocity of the plate; if not defined it uses bvel_velin specified elsewhere"
+    v_spread::Union{Float64,Nothing}    = nothing
+end
+
+function show(io::IO, d::geom_RidgeSeg)
+    println(io, "RidgeSeg(ph=$(d.phase), bounds=$(d.bounds), ridgeseg_x=$(d.ridgeseg_x),ridgeseg_y=$(d.ridgeseg_y),  T=$(d.Temperature) = [top,box]=[$(d.topTemp), $(d.botTemp)], age0=$(d.age0)), maxAge=$(d.maxAge), v_spread=$(d.v_spread))")
+    return nothing
+end
+
+"""
+    Write_LaMEM_InputFile(io, d::geom_RidgeSeg)
+"""
+function Write_LaMEM_InputFile(io, d::geom_RidgeSeg)
+    fields    = fieldnames(typeof(d))
+    println(io, "    <RidgeSegStart>")
+    for f in fields
+        if !isnothing(getfield(d,f))
+            name = rpad(String(f),15)
+            comment = get_doc(geom_RidgeSeg, f)
+            data = getfield(d,f) 
+            println(io,"        $name  = $(write_vec(data))     # $(comment)")
+        end
+    end
+    println(io, "    <RidgeSegEnd>")
+    return nothing
+end
+# -------
+
+# -------
+"""
+    LaMEM geometric primitive `Hex` object to define hexahedral elements
+    
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef struct geom_Hex
+    "phase"
+    phase::Int64          = 1          
+     
+    """
+    `x`-`y`-`z` coordinates for each of 8 nodes (24 parameters)
+    (counter)-clockwise for an arbitrary face, followed by the opposite face
+    """
+    coord::Vector{Float64}      = [0.25, 0.25, 0.25,   0.5, 0.2, 0.2,   0.6, 0.7, 0.25,   0.3, 0.5, 0.3,   0.2, 0.3, 0.75,   0.6, 0.15, 0.75,   0.5, 0.6, 0.80,   0.2, 0.4, 0.75]
+end
+
+function show(io::IO, d::geom_Hex)
+    println(io, "Hex(ph=$(d.phase), coord=$(d.coord))")
+    return nothing
+end
+
+"""
+    Write_LaMEM_InputFile(io, d::geom_Hex)
+"""
+function Write_LaMEM_InputFile(io, d::geom_Hex)
+    fields    = fieldnames(typeof(d))
+    println(io, "    <HexStart>")
+    for f in fields
+        if !isnothing(getfield(d,f))
+            name = rpad(String(f),15)
+            comment = get_doc(geom_Hex, f)
+            data = getfield(d,f) 
+            println(io,"        $name  = $(write_vec(data))     # $(comment)")
+        end
+    end
+    println(io, "    <HexEnd>")
+    return nothing
+end
+# -------
+
+# -------
+"""
+    LaMEM geometric primitive `Layer` object 
+    
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef struct geom_Layer
+    "phase"
+    phase::Int64          = 1          
+     
+    "top of layer"
+    top::Float64         = 5.0
+    
+    "bottom of layer"
+    bottom::Float64      = 3.0
+    
+    "optional: add a cosine perturbation on top of the interface (if 1)"
+    cosine::Union{Int64,Nothing} = nothing
+    
+    "required if cosine: wavelength in x-direction"
+    wavelength::Union{Float64,Nothing} = nothing
+
+    "required if cosine: amplitude of perturbation"
+    amplitude::Union{Float64,Nothing} = nothing
+
+    "optional: Temperature structure. possibilities: [constant, linear, halfspace]"
+    Temperature::Union{String,Nothing} = nothing
+
+    "required in case of [`constant`]: temperature value [in Celcius in case of GEO units]"
+    cstTemp::Union{Float64,Nothing}     = nothing   
+
+    "required in case of [`linear,halfspace`]: temperature @ top [in Celcius in case of GEO units]"
+    topTemp::Union{Float64,Nothing}     = nothing   
+    
+    "required in case of [`linear,halfspace`]: temperature @ top [in Celcius in case of GEO units]"
+    botTemp::Union{Float64,Nothing}     = nothing   
+    
+    "required in case of [`halfspace`]: thermal age of lithosphere [in Myrs if GEO units are used]"
+    thermalAge::Union{Float64,Nothing}  = nothing   
+end
+
+function show(io::IO, d::geom_Layer)
+    println(io, "Layer(ph=$(d.phase), bot/top=[$(d.bot),$(d.top)], T=$(d.Temperature)=$(d.cstTemp), [Ttop,Tbot]=[$(d.topTemp), $(d.botTemp)], thermalAge=$(d.thermalAge), cosine perturbation=$(d.cosine), wavelength=$(d.wavelength), amplitude=$(d.amplitude))")
+    return nothing
+end
+
+"""
+    Write_LaMEM_InputFile(io, d::geom_Layer)
+"""
+function Write_LaMEM_InputFile(io, d::geom_Layer)
+    fields    = fieldnames(typeof(d))
+    println(io, "    <LayerStart>")
+    for f in fields
+        if !isnothing(getfield(d,f))
+            name = rpad(String(f),15)
+            comment = get_doc(geom_Layer, f)
+            data = getfield(d,f) 
+            println(io,"        $name  = $(write_vec(data))     # $(comment)")
+        end
+    end
+    println(io, "    <LayerEnd>")
+    return nothing
+end
+# -------
+
+# -------
+"""
+    LaMEM geometric primitive `Cylinder` object 
+    
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef struct geom_Cylinder
+    "phase"
+    phase::Int64          = 1          
+     
+    "radius of cylinder "
+    radius::Float64      = 1.5
+    
+    "center of base of cylinder"
+    base::Vector{Float64}      = [1.0, 2.0, 3.0]
+    
+    "center of cap of cylinder"
+    cap::Vector{Float64}      = [3.0, 5.0, 7.0]
+    
+    "optional: Temperature structure. possibilities: [constant]"
+    Temperature::Union{String,Nothing} = nothing
+        
+    "required in case of [`constant`]: temperature value [in Celcius in case of GEO units]"
+    cstTemp::Union{Float64,Nothing}     = nothing   
+
+end
+
+function show(io::IO, d::geom_Cylinder)
+    println(io, "Cylinder(ph=$(d.phase), radius=$(d.radius) base=$(d.base), cap=$(d.cap), T=$(d.Temperature)=$(d.cstTemp) )")
+    return nothing
+end
+
+"""
+    Write_LaMEM_InputFile(io, d::geom_Cylinder)
+"""
+function Write_LaMEM_InputFile(io, d::geom_Cylinder)
+    fields    = fieldnames(typeof(d))
+    println(io, "    <CylinderStart>")
+    for f in fields
+        if !isnothing(getfield(d,f))
+            name = rpad(String(f),15)
+            comment = get_doc(geom_Cylinder, f)
+            data = getfield(d,f) 
+            println(io,"        $name  = $(write_vec(data))     # $(comment)")
+        end
+    end
+    println(io, "    <CylinderEnd>")
+    return nothing
+end
+# -------
+
 # ==========================
 
 
@@ -192,8 +523,6 @@ function Write_LaMEM_InputFile(io, d::ModelSetup)
         for object in d.geom_primitives
             Write_LaMEM_InputFile(io, object)
         end
-
-        println("WARNING! Still need to implement geometric primitives!!")
     end
 
 
diff --git a/src/LaMEM_ModelGeneration/Utils.jl b/src/LaMEM_ModelGeneration/Utils.jl
index ac072fa2..6b5f26e7 100644
--- a/src/LaMEM_ModelGeneration/Utils.jl
+++ b/src/LaMEM_ModelGeneration/Utils.jl
@@ -4,7 +4,7 @@ import LaMEM.IO_functions: Read_LaMEM_simulation, Read_LaMEM_timestep
 export  add_phase!, rm_phase!, rm_last_phase!, replace_phase!,
         add_vbox!, rm_last_vbox!, rm_vbox!,
         add_softening!, add_phasetransition!, add_phaseaggregate!,
-        add_dike!, add_geom! , cross_section,
+        add_dike!, add_geom!, rm_geom!, cross_section,
         set_air, copy_phase,
         isdefault, hasplasticity,
         stress_strainrate_0D
@@ -73,7 +73,6 @@ function add_phase!(model::Model, phases...)
 end
 
 
-
 """
     rm_last_phase!(model::Model, phase::Phase)
 This removes the last added `phase` from `model`
@@ -196,14 +195,43 @@ This adds an internal geometric primitive object `geom_object` to the LaMEM Mode
 
 Currently available primitive geom objects are:
 - `geom_Sphere`
+- `geom_Ellipsoid`
+- `geom_Box`
+- `geom_Layer`
+- `geom_Cylinder`
+- `geom_RidgeSeg`
+- `geom_Hex`
+
 """
 function add_geom!(model::Model, geom_object)
     push!(model.ModelSetup.geom_primitives, geom_object);
-    set_geom!(model, geom_object)
+    model.ModelSetup.msetup = "geom";
+    
+    #set_geom!(model, geom_object)
+
     return nothing
 end
 
 
+"""
+    add_geom!(model::Model, geom_object)
+Add several geometric objects @ once.
+"""
+function add_geom!(model::Model, geom_objects...) 
+    for geom_object in geom_objects
+        add_geom!(model, geom_object)
+    end
+end
+
+"""
+    rm_geom!(model::Model)
+This removes all existing geometric objects from `model`
+"""
+function rm_geom!(model::Model) 
+    model.ModelSetup.geom_primitives = []
+    return nothing
+end
+
 """
 
 This sets the geometry 
diff --git a/test/test_julia_setups.jl b/test/test_julia_setups.jl
index 375c0756..1689d866 100644
--- a/test/test_julia_setups.jl
+++ b/test/test_julia_setups.jl
@@ -161,7 +161,43 @@ using GeophysicalModelGenerator
 end
 
 
-@testset "various julia setups" begin
+@testset "phase diagrams" begin
     # This tests requires an update of LaMEM, to allow phase diagram names that are longer
     #include("test_julia_setup_phase_diagrams.jl")
+end
+
+
+
+@testset "build-in geometries" begin
+    # This tests the ability to write output files when we have build-in geometrical objects,
+    # such as spheres, boxes, etc. as opoosed to doing this through the GeophysicalModelGenerator
+    model  = Model(Grid(nel=(16,16), x=[-2,2], z=[-1,1]),
+                    Time(nstep_max=2, dt=1, dt_max=10), 
+                    Solver(SolverType="direct"),
+                    Output(out_dir="example_1"))
+
+    # Specify material properties
+    matrix = Phase(ID=0,Name="matrix",eta=1e20,rho=3000)
+    sphere = Phase(ID=1,Name="sphere",eta=1e23,rho=3200)
+    add_phase!(model, sphere, matrix)
+
+    geom_sphere = geom_Sphere();
+    rm_geom!(model)
+    add_geom!(model, geom_sphere)
+    out = run_lamem(model)
+    @test isnothing(out)
+
+    geom_ellipsoid = geom_Ellipsoid();
+    geom_box = geom_Box();
+    geom_layer = geom_Layer();
+    geom_cylinder = geom_Cylinder();
+    geom_ridge = geom_RidgeSeg();
+    geom_hex = geom_Hex();
+
+    add_geom!(model, geom_ellipsoid, geom_box, geom_layer, geom_cylinder, geom_ridge, geom_hex)
+    out = run_lamem(model)
+    @test isnothing(out)
+
+
+    
 end
\ No newline at end of file