diff --git a/.travis.yml b/.travis.yml
index c6d9b6106..72f4f34ec 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -20,6 +20,8 @@ jobs:
       julia: 1.1
       os: linux
       script:
+        - julia --project=docs/ -e 'using Pkg;
+                                    Pkg.add(PackageSpec(name="UnstructuredGrids", rev="master"))'
         - julia --project=docs/ -e 'using Pkg;
                                     Pkg.develop(PackageSpec(path=pwd()));
                                     Pkg.instantiate()'
diff --git a/Manifest.toml b/Manifest.toml
index d099b77d6..c6f2ae51f 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -197,7 +197,9 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [[UnstructuredGrids]]
 deps = ["WriteVTK"]
-git-tree-sha1 = "8d383ac05dcbbaa3b1222b441370e485e2abc831"
+git-tree-sha1 = "e780e3e9a2895f689cd68465221b866a6820112e"
+repo-rev = "master"
+repo-url = "https://github.com/gridap/UnstructuredGrids.jl.git"
 uuid = "9cd2e232-b78e-5341-8154-33bf4eafe965"
 version = "0.2.0"
 
diff --git a/src/CellValues/CellValuesGallery.jl b/src/CellValues/CellValuesGallery.jl
index d573e53ce..560f1b3f7 100644
--- a/src/CellValues/CellValuesGallery.jl
+++ b/src/CellValues/CellValuesGallery.jl
@@ -20,6 +20,7 @@ import Base: eltype
 import Base: size
 import Base: getindex
 import Base: IndexStyle
+import Gridap: compress
 
 struct CellValueFromArray{T,N,V<:AbstractArray{T,N}} <: IndexCellValue{T,N}
   v::V
@@ -58,6 +59,10 @@ size(self::CellVectorFromDataAndPtrs) = (length(self.ptrs)-1,)
 
 IndexStyle(::Type{CellVectorFromDataAndPtrs{T,V,P}}) where {T,V,P} = IndexLinear()
 
+function compress(cv::CellVectorFromDataAndPtrs)
+  (cv.data, cv.ptrs)
+end
+
 struct CellVectorFromDataAndStride{T,V} <: IndexCellArray{T,1,CachedSubVector{T,V},1}
   data::V
   stride::Int
diff --git a/src/CellValues/CellValuesOperations.jl b/src/CellValues/CellValuesOperations.jl
index 8b00fb854..7c804dc50 100644
--- a/src/CellValues/CellValuesOperations.jl
+++ b/src/CellValues/CellValuesOperations.jl
@@ -1,6 +1,7 @@
 module CellValuesOperations
 
 using Gridap
+using Gridap.Helpers
 
 using Gridap.Kernels: CellSumKernel
 using Gridap.Kernels: CellNewAxisKernel
@@ -8,6 +9,7 @@ using Gridap.Kernels: CellNewAxisKernel
 export cellsum
 export cellmean
 export cellnewaxis
+export compress
 import Base: +,-,*,/,\, ==, ≈
 import LinearAlgebra: inv, det
 import TensorValues: inner, outer, meas
@@ -111,4 +113,8 @@ function cellnewaxis(ca::CellMap;dim::Int)
   apply(k,ca)
 end
 
+function compress(cv::CellVector)
+  @notimplemented
+end
+
 end # module
diff --git a/src/CellValues/files.jl b/src/CellValues/files.jl
index 0d333fcc7..50301653b 100644
--- a/src/CellValues/files.jl
+++ b/src/CellValues/files.jl
@@ -29,12 +29,12 @@ include("MapApply.jl")
 include("CellMapApply.jl")
 @reexport using Gridap.CellMapApply
 
-include("CellValuesGallery.jl")
-#@reexport using Gridap.CellValuesGallery
-
 include("CellValuesOperations.jl")
 @reexport using Gridap.CellValuesOperations
 
+include("CellValuesGallery.jl")
+#@reexport using Gridap.CellValuesGallery
+
 include("CellValuesAppend.jl")
 @reexport using Gridap.CellValuesAppend
 
diff --git a/src/Geometry/GeometryWrappers.jl b/src/Geometry/GeometryWrappers.jl
index d91d82d2a..4a1f2eaca 100644
--- a/src/Geometry/GeometryWrappers.jl
+++ b/src/Geometry/GeometryWrappers.jl
@@ -23,9 +23,41 @@ RefCell(polytope::Polytope) = _ref_cell_from_polytope(polytope)
 Construct a Grid from a polytope
 """
 function Grid(polytope::Polytope{D},dim::Int) where D
+  if dim < D
+    _grid_d(polytope,dim)
+  elseif dim==D
+    _grid_D(polytope)
+  else
+    @unreachable
+  end
+end
+
+# Helpers
+
+function _grid_D(polytope::Polytope{D}) where D
+
+  order = 1
+  reffe = LagrangianRefFE{D,Float64}(polytope,order)
+  points = reffe.dofbasis.nodes
+
+  cell_to_nodes = [ [i for i in 1:length(points)], ]
+
+  cells_data,cells_ptrs = generate_data_and_ptrs(cell_to_nodes)
+  extrusion = polytope.extrusion.array.data
+
+  l = 1
+  co = ConstantCellValue(order,l)
+  ct = ConstantCellValue(extrusion,l)
+
+  UnstructuredGrid(points,cells_data,cells_ptrs,ct,co)
+
+end
+
+function _grid_d(polytope::Polytope{D},dim::Int) where D
   @assert dim < D
 
-  reffe = LagrangianRefFE{D,Float64}(polytope, 1)
+  order = 1
+  reffe = LagrangianRefFE{D,Float64}(polytope,order)
   points = reffe.dofbasis.nodes
 
   dim_to_jface_to_vertices, dim_to_jface_to_code = _faces(polytope)
@@ -49,7 +81,7 @@ function Grid(polytope::Polytope{D},dim::Int) where D
 
 end
 
-# Helpers
+
 
 function _ref_cell_from_polytope(polytope::Polytope{D}) where D
 
diff --git a/src/Geometry/Simplexify.jl b/src/Geometry/Simplexify.jl
new file mode 100644
index 000000000..75f4749f0
--- /dev/null
+++ b/src/Geometry/Simplexify.jl
@@ -0,0 +1,254 @@
+module Simplexify
+
+using Gridap
+using Gridap.Helpers
+using UnstructuredGrids.Kernels: refine_grid_connectivity
+using UnstructuredGrids.Kernels: generate_tface_to_face
+using Gridap.DiscreteModels: DiscreteModelFromData
+
+export simplexify
+
+function simplexify(grid::Grid{D,Z}) where {D,Z}
+  ugrid = UnstructuredGrid(grid)
+
+  ltcell_to_lpoints = _generate_ltcell_to_lpoints(celltypes(ugrid))
+
+  tcells_data, tcells_ptrs = refine_grid_connectivity(
+    ugrid.cells_data, ugrid.cells_ptrs, ltcell_to_lpoints)
+
+  order = 1
+  _check_order(order,cellorders(ugrid))
+
+  ntcells = length(tcells_ptrs) - 1
+  ex = tuple(fill(TET_AXIS,Z)...)
+  _ct = ConstantCellValue(ex,ntcells)
+  _co = ConstantCellValue(order,ntcells)
+
+  UnstructuredGrid(
+    points(ugrid),
+    tcells_data,
+    tcells_ptrs,
+    _ct,
+    _co)
+
+end
+
+function simplexify(model::DiscreteModel{D}) where D
+
+  grid = Grid(model)
+
+  graph = FullGridGraph(model)
+
+  facelabels = FaceLabels(model)
+
+  tgrid = simplexify(grid)
+
+  tgraph = FullGridGraph(tgrid)
+
+  tgrids = _generate_tgrids(tgrid,tgraph)
+
+  ltcell_to_lnodes = _generate_ltcell_to_lpoints(celltypes(grid))
+
+  tfacelabels = _generate_tfacelabels(
+    tgrids,grid,facelabels,tgraph,graph,ltcell_to_lnodes)
+
+  DiscreteModelFromData(tgrids,tgraph,tfacelabels)
+
+end
+
+# Helpers
+
+const UNSET_ID = 0
+
+function _generate_tfacelabels(
+  tgrids,grid,facelabels,tgraph,graph,ltcell_to_lnodes)
+
+  dim_to_tface_to_label = [ fill(UNSET_ID,ncells(g)) for g in tgrids ]
+
+  dim_to_face_to_label = facelabels.dim_to_nface_to_label
+
+  tgrid = tgrids[end]
+  tpolytope = CellPolytopes(tgrid).value
+  polytope = CellPolytopes(grid).value
+
+  _fill_dim_to_tface_to_label!(
+    dim_to_tface_to_label,
+    dim_to_face_to_label,
+    tpolytope,
+    polytope,
+    tgraph,
+    graph,
+    ltcell_to_lnodes)
+
+  FaceLabels(
+    dim_to_tface_to_label, facelabels.tag_to_labels, facelabels.tag_to_name)
+
+end
+
+function _fill_dim_to_tface_to_label!(
+  dim_to_tface_to_label,
+  dim_to_face_to_label,
+  tpolytope,
+  polytope,
+  tgraph,
+  graph,
+  ltcell_to_lnodes)
+
+  D = length(dim_to_face_to_label)-1
+  d = 0
+  dim_to_tface_to_label[d+1] = dim_to_face_to_label[d+1]
+
+  trefcell = RefCell(tpolytope)
+  refcell = RefCell(polytope)
+
+  for d in 1:(D-1)
+
+    cell_to_faces = connections(graph,D,d)
+    cell_to_faces_data, cell_to_faces_ptrs = compress(cell_to_faces)
+
+    tcell_to_tfaces = connections(tgraph,D,d)
+
+    ntfaces = maximum(tcell_to_tfaces.data)
+
+    ltface_to_ltnodes = trefcell.faces[d+1]
+    lface_to_lnodes = refcell.faces[d+1]
+
+    tface_to_face = generate_tface_to_face(
+      cell_to_faces_data,
+      cell_to_faces_ptrs,
+      tcell_to_tfaces.data,
+      tcell_to_tfaces.ptrs,
+      ltcell_to_lnodes,
+      ltface_to_ltnodes,
+      lface_to_lnodes,
+      ntfaces)
+
+    _update_labels!(
+      dim_to_tface_to_label[d+1],dim_to_face_to_label[d+1],tface_to_face)
+
+  end
+
+  d = D
+  ncells = length(dim_to_face_to_label[d+1])
+  nltcells = length(ltcell_to_lnodes)
+  tcell_to_cell = _generate_tcell_to_cell(ncells,nltcells)
+  _update_labels!(
+    dim_to_tface_to_label[d+1],dim_to_face_to_label[d+1],tcell_to_cell)
+
+  for d = 1:(D-1)
+
+    for j in (d+1):D
+
+      dface_to_jfaces = connections(tgraph,d,j)
+      dface_to_label = dim_to_tface_to_label[d+1]
+      jface_to_label = dim_to_tface_to_label[j+1]
+      _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)
+
+    end
+
+  end
+
+end
+
+function _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)
+
+  ndfaces = length(dface_to_label)
+  @assert ndfaces == length(dface_to_jfaces)
+
+  for dface in 1:ndfaces
+
+    dlabel = dface_to_label[dface]
+    if dlabel != UNSET_ID
+      continue
+    end
+
+    jfaces = dface_to_jfaces[dface]
+    for jface in jfaces
+      jlabel = jface_to_label[jface]
+      if jlabel != UNSET_ID
+        dface_to_label[dface] = jlabel
+        break
+      end
+    end
+
+  end
+
+end
+
+function _update_labels!(tface_to_label,face_to_label,tface_to_face)
+  for tface in 1:length(tface_to_label)
+    face = tface_to_face[tface]
+    if face != 0
+      tface_to_label[tface] = face_to_label[face]
+    end
+  end
+end
+
+function _generate_tcell_to_cell(ncells,nltcells)
+  ntcells = ncells * nltcells
+  tcell_to_cell = Vector{Int}(undef,ntcells)
+  tcell = 1
+  for cell in 1:ncells
+    for ltcell in 1:nltcells
+      tcell_to_cell[tcell] = cell
+      tcell += 1
+    end
+  end
+  tcell_to_cell
+end
+
+function _generate_tgrids(tgrid::Grid{D},tgraph) where D
+  tgrids = Grid[]
+  for d in 0:(D-1)
+    tface_to_vertices = connections(tgraph,d,0)
+    tgrid_d = _grid_of_dim(points(tgrid),tface_to_vertices,d)
+    push!(tgrids,tgrid_d)
+  end
+  push!(tgrids,tgrid)
+  tgrids
+end
+
+function _grid_of_dim(coords,face_to_vertices,Z)
+
+  nfaces = length(face_to_vertices)
+  fcode = tuple([TET_AXIS for i in 1:Z]...)
+
+  order = 1
+
+  _points = coords
+  _cells_data = face_to_vertices.data
+  _cells_ptrs = face_to_vertices.ptrs
+  _ctypes = ConstantCellValue(fcode,nfaces)
+  _corders = ConstantCellValue(order,nfaces)
+
+  UnstructuredGrid(
+    _points, _cells_data, _cells_ptrs, _ctypes, _corders)
+
+end
+
+function _check_order(order,co)
+  @notimplemented
+end
+
+function _check_order(order,co::ConstantCellValue)
+  @notimplementedif co.value != order
+end
+
+function _generate_ltcell_to_lpoints(ct)
+  @notimplemented
+end
+
+function _generate_ltcell_to_lpoints(ct::ConstantCellValue)
+  extrusion = ct.value
+  if extrusion == (HEX_AXIS, HEX_AXIS)
+    ltcell_to_lpoints = [[1,2,3],[4,3,2]]
+    return ltcell_to_lpoints
+  elseif extrusion == (HEX_AXIS, HEX_AXIS, HEX_AXIS)
+    ltcell_to_lpoints = [
+      [7,3,2,1], [7,5,2,1], [7,4,3,2], [7,4,8,2], [7,6,5,2], [7,6,8,2]]
+  else
+    @notimplemented
+  end
+end
+
+end # module
diff --git a/src/Geometry/UnstructuredGeometry.jl b/src/Geometry/UnstructuredGeometry.jl
index 61f47371a..be41531ef 100644
--- a/src/Geometry/UnstructuredGeometry.jl
+++ b/src/Geometry/UnstructuredGeometry.jl
@@ -55,6 +55,12 @@ function UnstructuredGrid(
     corders)
 end
 
+function UnstructuredGrid(grid::Grid)
+  @notimplemented
+end
+
+UnstructuredGrid(grid::UnstructuredGrid) = grid
+
 points(self::UnstructuredGrid) = CellValueFromArray(self.points)
 
 cells(self::UnstructuredGrid) = CellVectorFromDataAndPtrs(self.cells_data,self.cells_ptrs)
diff --git a/src/Geometry/files.jl b/src/Geometry/files.jl
index df5310a6d..261518269 100644
--- a/src/Geometry/files.jl
+++ b/src/Geometry/files.jl
@@ -28,3 +28,6 @@ include("BoundaryGrids.jl")
 
 include("SkeletonGrids.jl")
 @reexport using Gridap.SkeletonGrids
+
+include("Simplexify.jl")
+@reexport using Gridap.Simplexify
diff --git a/src/RefFEs/RefFEs.jl b/src/RefFEs/RefFEs.jl
index 8d761a1f8..da965c747 100644
--- a/src/RefFEs/RefFEs.jl
+++ b/src/RefFEs/RefFEs.jl
@@ -61,10 +61,26 @@ struct LagrangianRefFE{D,T} <: RefFE{D,T}
   # this type is unstable
 end
 
+#@fverdugo. This is only a temporary hack to make work the interface of
+#anisotropic order also for n-simplices. Not that an @notimplemented error
+# will be raised for simplices if all orders are not the same.
+function LagrangianRefFE{D,T}(polytope::Polytope{D}, orders::Vector{Int64}) where {D,T}
+  @assert length(orders) == D
+  @assert D > 0
+  if all(polytope.extrusion.array .== HEX_AXIS)
+    return _LagrangianRefFE(T, polytope, orders)
+  elseif all( orders .== orders[1] )
+    order = orders[1]
+    return LagrangianRefFE{D,T}(polytope,order)
+  else
+    @notimplemented
+  end
+end
+
 # @santiagobadia: Temporary constructor that uses the old NodeArray.
 # I keep it because I am not considering anisotropic order
 # in the new method. Future development.
-function LagrangianRefFE{D,T}(polytope::Polytope{D},
+function _LagrangianRefFE(::Type{T},polytope::Polytope{D},
   orders::Vector{Int64}) where {D,T}
   nodes=NodesArray(polytope,orders)
   dofsb = LagrangianDOFBasis{D,T}(nodes.coordinates)
diff --git a/src/Visualization/Vtkio.jl b/src/Visualization/Vtkio.jl
index b3f33f024..711f8a184 100644
--- a/src/Visualization/Vtkio.jl
+++ b/src/Visualization/Vtkio.jl
@@ -340,6 +340,12 @@ end
 
 function _refgrid end
 
+function _refgrid(::Val{V},nref::Int) where V
+  @notimplementedif nref != 0
+  polytope = Polytope(V)
+  Grid(polytope,length(V))
+end
+
 function _refgrid(::Val{(HEX_AXIS,)},nref::Int)
   n = 2^nref
   CartesianGrid(domain=(0.0,1.0),partition=(n,))
diff --git a/test/CellValuesTests/CellValuesGalleryTests.jl b/test/CellValuesTests/CellValuesGalleryTests.jl
index 4c14207f9..60f0facce 100644
--- a/test/CellValuesTests/CellValuesGalleryTests.jl
+++ b/test/CellValuesTests/CellValuesGalleryTests.jl
@@ -32,6 +32,10 @@ ca = CellVectorFromDataAndPtrs(data,ptrs)
 a = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
 test_index_cell_array( ca, a )
 
+d,p = compress(ca)
+@test d === ca.data
+@test p === ca.ptrs
+
 a = [ data[3*(i-1)+1:3*i] for i in 1:4]
 ca = CellVectorFromDataAndStride(data,3)
 test_index_cell_array( ca, a )
diff --git a/test/GeometryTests/GeometryWrappersTests.jl b/test/GeometryTests/GeometryWrappersTests.jl
index 4972599a2..a2be3b59a 100644
--- a/test/GeometryTests/GeometryWrappersTests.jl
+++ b/test/GeometryTests/GeometryWrappersTests.jl
@@ -6,10 +6,20 @@ using Test
 
 # Polytope to UnstructuredGrid
 
+t = TET_AXIS
+polytope = Polytope((t,t,t))
+
+grid = Grid(polytope,1)
+grid = Grid(polytope,2)
+grid = Grid(polytope,3)
+#writevtk(grid,"grid")
+
 t = HEX_AXIS
 polytope = Polytope((t,t,t))
 
+grid = Grid(polytope,1)
 grid = Grid(polytope,2)
+grid = Grid(polytope,3)
 #writevtk(grid,"grid")
 
 trian = Triangulation(grid)
diff --git a/test/GeometryTests/SimplexifyTests.jl b/test/GeometryTests/SimplexifyTests.jl
new file mode 100644
index 000000000..249b44c13
--- /dev/null
+++ b/test/GeometryTests/SimplexifyTests.jl
@@ -0,0 +1,23 @@
+module SimplexifyTests
+
+using Gridap
+
+grid = CartesianGrid(domain=(0.0,1.0,-1.0,2.0),partition=(3,4))
+tgrid = simplexify(grid)
+test_grid(tgrid,20,24)
+#writevtk(tgrid,"tgrid")
+
+grid = CartesianGrid(partition=(3,4,3))
+tgrid = simplexify(grid)
+test_grid(tgrid,80,216)
+#writevtk(tgrid,"tgrid")
+
+model = CartesianDiscreteModel(partition=(3,4,2))
+tmodel = simplexify(model)
+
+model = CartesianDiscreteModel(partition=(3,4))
+tmodel = simplexify(model)
+#writevtk(tmodel,"tmodel")
+
+
+end # module
diff --git a/test/GeometryTests/runtests.jl b/test/GeometryTests/runtests.jl
index a3a769d8c..75b9e3c89 100644
--- a/test/GeometryTests/runtests.jl
+++ b/test/GeometryTests/runtests.jl
@@ -22,4 +22,6 @@ using Test
 
 @testset "SkeletonGrids" begin include("SkeletonGridsTests.jl") end
 
+@testset "SimplexifyTests" begin include("SimplexifyTests.jl") end
+
 end # module
diff --git a/test/VisualizationTests/VtkioTests.jl b/test/VisualizationTests/VtkioTests.jl
index 2d6d01506..0f277b873 100644
--- a/test/VisualizationTests/VtkioTests.jl
+++ b/test/VisualizationTests/VtkioTests.jl
@@ -135,6 +135,11 @@ writevtk(trian,f,nref=2,celldata=["r"=>rand(9)])
 writevtk(trian,f,nref=2,cellfields=["u"=>u,"v"=>v])
 writevtk(trian,f,nref=2,celldata=["r"=>rand(9)],cellfields=["u"=>u,"v"=>v])
 
+grid = CartesianGrid(partition=(3,3))
+grid = simplexify(grid)
+trian = Triangulation(grid)
+writevtk(trian,f)
+
 rm(d,recursive=true)
 
 end