diff --git a/Project.toml b/Project.toml
index 6fe569e..1b060eb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,7 +8,6 @@ AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
-GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -33,7 +32,6 @@ BlockArrays = "1.9"
 ForwardDiff = "1.3"
 Gridap = "0.19.8"
 GridapDistributed = "0.4"
-GridapGmsh = "0.7"
 GridapPETSc = "0.5"
 GridapSolvers = "0.6"
 IterativeSolvers = "0.9"
diff --git a/src/ComputationalModels/ComputationalModels.jl b/src/ComputationalModels/ComputationalModels.jl
index dc73ac4..e44badf 100644
--- a/src/ComputationalModels/ComputationalModels.jl
+++ b/src/ComputationalModels/ComputationalModels.jl
@@ -19,9 +19,6 @@ using GridapSolvers.SolverInterfaces: finished, print_message, converged
 using LinearAlgebra
 using WriteVTK
 
-using GridapGmsh 
-using GridapGmsh: GmshDiscreteModel
-
 import Base.getindex
 
 include("BoundaryConditions.jl")
@@ -48,7 +45,6 @@ export TestFESpace
 export TestFESpace!
 
 include("GridapExtras.jl")
-export GmshDiscreteModel
 export evaluate!
 
 include("Drivers.jl")
diff --git a/src/Exports.jl b/src/Exports.jl
index 3818ab2..6cb9034 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -131,7 +131,6 @@ end
 @publish ComputationalModels  StaticLinearModel
 @publish ComputationalModels  solve!
 @publish ComputationalModels  dirichlet_preconditioning!
-@publish ComputationalModels  GmshDiscreteModel
 @publish ComputationalModels  updateBC!
 @publish ComputationalModels  PostProcessor
 @publish ComputationalModels  vtk_save
diff --git a/test/Project.toml b/test/Project.toml
index d898967..dba39ca 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/test/SimulationsTests/BoundaryConditions.jl b/test/SimulationsTests/BoundaryConditions.jl
index 0149727..5215687 100644
--- a/test/SimulationsTests/BoundaryConditions.jl
+++ b/test/SimulationsTests/BoundaryConditions.jl
@@ -1,4 +1,4 @@
-using Gridap
+using Gridap, GridapGmsh
 using Gridap.FESpaces
 
 @testset "Dirichlet BC testing analytical mapping" begin
diff --git a/test/TestWeakForms/WeakForms.jl b/test/TestWeakForms/WeakForms.jl
index 79c16ca..7709cca 100644
--- a/test/TestWeakForms/WeakForms.jl
+++ b/test/TestWeakForms/WeakForms.jl
@@ -129,76 +129,3 @@ end
 
 
 end
-
-# @testset "Assembly Jacobian Mechanics" begin
-
-#     pname = "test"
-#     meshfile = "test_cubicsphere.msh"
-#     simdir = datadir("sims", pname)
-#     setupfolder(simdir)
-  
-#     geomodel = GmshDiscreteModel(datadir("models", meshfile))
-  
-#     # Constitutive models
-#     ∇umacro = TensorValue(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0) * 1e-2
-  
-#     Kmodel = EvolutiveKinematics(Mechano; F=(t) -> ((∇u) -> ∇u + one(∇u) + t * ∇umacro))
-#     physmodel_matrix = MooneyRivlin3D(λ=10.0, μ1=1.0, μ2=1.0, Kinematic=Kmodel)
-#     physmodel_inclussion = MooneyRivlin3D(λ=50.0, μ1=5.0, μ2=5.0, Kinematic=Kmodel)
-  
-#     # Setup integration
-#     order = 1
-#     degree = 2 * order
-#     Ω = Triangulation(geomodel)
-#     dΩ = Measure(Ω, degree)
-#     Ω₁ = Interior(geomodel, tags="Phase0")
-#     dΩ₁ = Measure(Ω₁, degree)
-#     Ω₂ = Interior(geomodel, tags="Inclusions")
-#     dΩ₂ = Measure(Ω₂, degree)
-  
-#     # Dirichlet conditions 
-#     evolu(Λ) = 1.0
-#     dir_u_tags = ["Corners"]
-#     dir_u_values = [[0.0, 0.0, 0.0]]
-#     dir_u_timesteps = [evolu]
-#     D_bc = DirichletBC(dir_u_tags, dir_u_values, dir_u_timesteps)
-  
-#     #  FE spaces
-#     reffeu = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
-#     V = TestFESpace(geomodel, reffeu, D_bc, :H1)
-#     U = TrialFESpace(V, D_bc, 1.0)
-  
-#     #  residual and jacobian function of load factor
-#     res(Λ) = (u, v) -> residual(physmodel_matrix, u, v, dΩ₁, Λ) + residual(physmodel_inclussion, u, v, dΩ₂, Λ)
-#     jac(Λ) = (u, du, v) -> jacobian(physmodel_matrix, u, du, v, dΩ₁, Λ) + jacobian(physmodel_inclussion, u, du, v, dΩ₂, Λ)
-  
-#     uh = FEFunction(V, zeros(Float64, num_free_dofs(V)))
-
-#     function jach(uh, φh)
-#         jac((du, dφ), (v, vφ)) = jacobian(modelelectro, (uh, φh), (du, dφ), (v, vφ), dΩ)
-#     end
-#     function jac_mech(uh, φh)
-#         jac(du, v) = jacobian(modelelectro, Mechano, (uh, φh), du, v, dΩ)
-#     end
-#     function jac_elech(uh, φh)
-#         jac(dφ, vφ) = jacobian(modelelectro, Electro, (uh, φh), dφ, vφ, dΩ)
-#     end
-
-#     jac_ = assemble_matrix(jach(uh, φh), V, V)
-#     jac_m = assemble_matrix(jac_mech(uh, φh), Vu, Vu)
-#     jac_e = assemble_matrix(jac_elech(uh, φh), Vφ, Vφ)
-#     ≈
-#     @test norm(jac_) ≈ 18.934585248125135
-#     @test jac_[1] ≈ 0.7777777777777775
-#     @test jac_[end] ≈ -1.3333333333333326
-
-#     @test norm(jac_m) ≈ 16.420402846143244
-#     @test jac_m[1] ≈ 0.7777777777777775
-#     @test jac_m[end] ≈ 2.1111111111111094
-
-#     @test norm(jac_e) ≈ 4.216370213557837
-#     @test jac_e[1] ≈ -1.3333333333333321
-#     @test jac_e[end] ≈ -1.3333333333333326
-
-
-# end
\ No newline at end of file
diff --git a/test/data/Stokes.jl b/test/data/Stokes.jl
index 8a1ac2d..0f096f7 100644
--- a/test/data/Stokes.jl
+++ b/test/data/Stokes.jl
@@ -1,4 +1,4 @@
-using Gridap, GridapSolvers
+using Gridap, GridapSolvers, GridapGmsh
 using GridapSolvers.LinearSolvers
 using TimerOutputs
 using Gridap.FESpaces