diff --git a/src/Arrays/AppendedArrays.jl b/src/Arrays/AppendedArrays.jl
new file mode 100644
index 000000000..9b86f5d77
--- /dev/null
+++ b/src/Arrays/AppendedArrays.jl
@@ -0,0 +1,118 @@
+
+"""
+"""
+function lazy_append(a::AbstractArray,b::AbstractArray)
+  AppendedArray(a,b)
+end
+
+"""
+"""
+function lazy_split(f::AbstractArray,n::Integer)
+  _lazy_split(f,n)
+end
+
+function lazy_split(f::CompressedArray,n::Integer)
+  _a , _b =_lazy_split(f,n)
+  a = _compact_values_ptrs(_a)
+  b = _compact_values_ptrs(_b)
+  a,b
+end
+
+function _lazy_split(f,n)
+  l = length(f)
+  @assert n <= l
+  ids_a = collect(1:n)
+  ids_b = collect((n+1):l)
+  a = reindex(f,ids_a)
+  b = reindex(f,ids_b)
+  a,b
+end
+
+function _compact_values_ptrs(a)
+  v = a.values
+  p = a.ptrs
+  touched = fill(false,length(v))
+  for i in p
+    touched[i] = true
+  end
+  if all(touched) || length(p) == 0
+    return a
+  else
+    j_to_i = findall(touched)
+    i_to_j = zeros(Int,length(v))
+    i_to_j[j_to_i] = 1:length(j_to_i)
+    _v = v[j_to_i]
+    _p = i_to_j[p]
+    return CompressedArray(_v, _p)
+  end
+end
+
+struct AppendedArray{T,A,B} <: AbstractVector{T}
+  a::A
+  b::B
+  function AppendedArray(a::AbstractArray,b::AbstractArray)
+    A = typeof(a)
+    B = typeof(b)
+    ai = testitem(a)
+    bi = testitem(b)
+    T = eltype(collect((ai,bi)))
+    new{T,A,B}(a,b)
+  end
+end
+
+Base.IndexStyle(::Type{<:AppendedArray}) = IndexLinear()
+
+function Base.getindex(v::AppendedArray,i::Integer)
+  l = length(v.a)
+  if i > l
+    v.b[i-l]
+  else
+    v.a[i]
+  end
+end
+
+Base.size(v::AppendedArray) = (length(v.a)+length(v.b),)
+
+function array_cache(v::AppendedArray)
+  ca = array_cache(v.a)
+  cb = array_cache(v.b)
+  (ca,cb)
+end
+
+function getindex!(cache,v::AppendedArray,i::Integer)
+  ca, cb = cache
+  l = length(v.a)
+  if i > l
+    getindex!(cb,v.b,i-l)
+  else
+    getindex!(ca,v.a,i)
+  end
+end
+
+function apply(f,a::AppendedArray...)
+  la = map(ai->length(ai.a),a)
+  lb = map(ai->length(ai.b),a)
+  if all(la .== first(la)) && all(lb .== first(lb))
+    c_a = apply(f,map(ai->ai.a,a)...)
+    c_b = apply(f,map(ai->ai.b,a)...)
+    lazy_append(c_a,c_b)
+  else
+    s = common_size(a...)
+    apply(Fill(f,s...),a...)
+  end
+end
+
+function apply(::Type{T},f,a::AppendedArray...) where T
+  la = map(ai->length(ai.a),a)
+  lb = map(ai->length(ai.b),a)
+  if all(la .== first(la)) && all(lb .== first(lb))
+    c_a = apply(T,f,map(ai->ai.a,a)...)
+    c_b = apply(T,f,map(ai->ai.b,a)...)
+    lazy_append(c_a,c_b)
+  else
+    s = common_size(a...)
+    apply(T,Fill(f,s...),a...)
+  end
+end
+
+
diff --git a/src/Arrays/Arrays.jl b/src/Arrays/Arrays.jl
index ef83816a7..8f18113f9 100644
--- a/src/Arrays/Arrays.jl
+++ b/src/Arrays/Arrays.jl
@@ -86,6 +86,10 @@ export unpair_arrays
 
 export matvec_muladd!
 
+export lazy_append
+export lazy_split
+export AppendedArray
+
 import Base: size
 import Base: getindex, setindex!
 import Base: similar
@@ -118,4 +122,6 @@ include("SubVectors.jl")
 
 include("ArrayPairs.jl")
 
+include("AppendedArrays.jl")
+
 end # module
diff --git a/src/Arrays/Interface.jl b/src/Arrays/Interface.jl
index e75203120..64ac7ab76 100644
--- a/src/Arrays/Interface.jl
+++ b/src/Arrays/Interface.jl
@@ -188,8 +188,8 @@ function test_array(
   t = true
   for i in eachindex(a)
     ai = getindex!(cache,a,i)
-    t = t && (typeof(ai) == eltype(a))
-    t = t && (typeof(ai) == T)
+    t = t && (typeof(ai) <: eltype(a))
+    t = t && (typeof(ai) <: T)
   end
   @test t
   @test IndexStyle(a) == IndexStyle(b)
diff --git a/src/FESpaces/ExtendedFESpaces.jl b/src/FESpaces/ExtendedFESpaces.jl
index ef1865570..f9eede12c 100644
--- a/src/FESpaces/ExtendedFESpaces.jl
+++ b/src/FESpaces/ExtendedFESpaces.jl
@@ -77,11 +77,15 @@ function reindex(a::ExtendedVector,ptrs::IdentityVector)
   a
 end
 
-function reindex(a::ExtendedVector,trian::Triangulation)
-  ptrs = get_cell_id(trian)
-  _extended_reindex(a,ptrs)
+function reindex(a::ExtendedVector,ptrs::SkeletonPair)
+  _extended_reindex(a,ptrs::SkeletonPair)
 end
 
+#function reindex(a::ExtendedVector,trian::Triangulation)
+#  ptrs = get_cell_id(trian)
+#  _extended_reindex(a,ptrs)
+#end
+
 function _extended_reindex(a,ptrs::SkeletonPair)
   left= _extended_reindex(a,ptrs.left)
   right = _extended_reindex(a,ptrs.right)
diff --git a/src/Fields/FieldArrays.jl b/src/Fields/FieldArrays.jl
index a928b6d2e..21dd0d4ac 100644
--- a/src/Fields/FieldArrays.jl
+++ b/src/Fields/FieldArrays.jl
@@ -284,3 +284,26 @@ for op in (:+,:-)
   end
 end
 
+# More optimizations
+
+function evaluate_field_array(a::AppendedArray,b::AppendedArray)
+  if (length(a.a) == length(b.a)) && (length(a.b) == length(b.b))
+    c_a = evaluate_field_array(a.a,b.a)
+    c_b = evaluate_field_array(a.b,b.b)
+    lazy_append(c_a,c_b)
+  else
+    _evaluate_field_array(a,b)
+  end
+end
+
+function evaluate_field_array(a::AppendedArray,b::AbstractArray)
+  n = length(a.a)
+  _b = lazy_append(lazy_split(b,n)...)
+  evaluate_field_array(a,_b)
+end
+
+function field_array_gradient(a::AppendedArray)
+  c_a = field_array_gradient(a.a)
+  c_b = field_array_gradient(a.b)
+  lazy_append(c_a,c_b)
+end
diff --git a/src/Geometry/AppendedTriangulations.jl b/src/Geometry/AppendedTriangulations.jl
new file mode 100644
index 000000000..7fe20414b
--- /dev/null
+++ b/src/Geometry/AppendedTriangulations.jl
@@ -0,0 +1,79 @@
+
+function lazy_append(a::Triangulation,b::Triangulation)
+  AppendedTriangulation(a,b)
+end
+
+struct AppendedTriangulation{Dc,Dp} <: Triangulation{Dc,Dp}
+  a::Triangulation{Dc,Dp}
+  b::Triangulation{Dc,Dp}
+end
+
+function get_cell_coordinates(trian::AppendedTriangulation)
+  a = get_cell_coordinates(trian.a)
+  b = get_cell_coordinates(trian.b)
+  lazy_append(a,b)
+end
+
+function get_reffes(trian::AppendedTriangulation)
+  vcat(get_reffes(trian.a),get_reffes(trian.b))
+end
+
+function get_cell_type(trian::AppendedTriangulation)
+  a = get_cell_type(trian.a)
+  b = get_cell_type(trian.b) .+ Int8(length(get_reffes(trian.a)))
+  lazy_append(a,b)
+end
+
+function reindex(f::AbstractArray,trian::AppendedTriangulation)
+  a = reindex(f,trian.a)
+  b = reindex(f,trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_id(trian::AppendedTriangulation)
+  a = get_cell_id(trian.a)
+  b = get_cell_id(trian.b)
+  lazy_append(a,b)
+end
+
+function restrict(f::AbstractArray,trian::AppendedTriangulation)
+  a = restrict(f,trian.a)
+  b = restrict(f,trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_reffes(trian::AppendedTriangulation)
+  a = get_cell_reffes(trian.a)
+  b = get_cell_reffes(trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_shapefuns(trian::AppendedTriangulation)
+  a = get_cell_shapefuns(trian.a)
+  b = get_cell_shapefuns(trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_map(trian::AppendedTriangulation)
+  a = get_cell_map(trian.a)
+  b = get_cell_map(trian.b)
+  lazy_append(a,b)
+end
+
+function get_normal_vector(trian::AppendedTriangulation)
+  cm = get_cell_map(trian)
+  a = get_array(get_normal_vector(trian.a))
+  b = get_array(get_normal_vector(trian.b))
+  GenericCellField(lazy_append(a,b),cm)
+end
+
+function CellQuadrature(trian::AppendedTriangulation,degree1::Integer,degree2::Integer)
+  quad1 = CellQuadrature(trian.a,degree1)
+  quad2 = CellQuadrature(trian.b,degree2)
+  lazy_append(quad1,quad2)
+end
+
+function CellQuadrature(trian::AppendedTriangulation,degree::Integer)
+  CellQuadrature(trian,degree,degree)
+end
+
diff --git a/src/Geometry/CartesianGrids.jl b/src/Geometry/CartesianGrids.jl
index e49699533..3675d477e 100644
--- a/src/Geometry/CartesianGrids.jl
+++ b/src/Geometry/CartesianGrids.jl
@@ -215,6 +215,11 @@ function field_array_gradient(a::CartesianMap)
   Fill(j,length(a))
 end
 
+function field_array_gradient(a::Reindexed{T,N,A}) where {T,N,A<:CartesianMap}
+  g = field_array_gradient(a.i_to_v)
+  reindex(g,a.j_to_i)
+end
+
 function get_cell_map(grid::CartesianGrid{D,T,typeof(identity)} where {D,T})
   CartesianMap(grid.node_coords.data)
 end
diff --git a/src/Geometry/CellQuadratures.jl b/src/Geometry/CellQuadratures.jl
index 907d6be20..770d8f176 100644
--- a/src/Geometry/CellQuadratures.jl
+++ b/src/Geometry/CellQuadratures.jl
@@ -77,6 +77,12 @@ function _get_coordinates(q::Fill{<:Quadrature})
   Fill(coords,length(q))
 end
 
+function _get_coordinates(q::AppendedArray)
+  a = _get_coordinates(q.a)
+  b = _get_coordinates(q.b)
+  lazy_append(a,b)
+end
+
 function _get_weights(q::AbstractArray{<:Quadrature})
   @notimplemented "Not implemented, since we dont need it"
 end
@@ -91,6 +97,12 @@ function _get_weights(q::Fill{<:Quadrature})
   Fill(w,length(q))
 end
 
+function _get_weights(q::AppendedArray)
+  a = _get_weights(q.a)
+  b = _get_weights(q.b)
+  lazy_append(a,b)
+end
+
 """
     integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
 
@@ -108,3 +120,8 @@ function integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
   integrate(get_array(f),q,w,j)
 end
 
+function lazy_append(quad1::CellQuadrature,quad2::CellQuadrature)
+  array = lazy_append(quad1.array,quad2.array)
+  CellQuadrature(array)
+end
+
diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl
index 859d152ce..9922c7ed2 100644
--- a/src/Geometry/Geometry.jl
+++ b/src/Geometry/Geometry.jl
@@ -27,10 +27,13 @@ using Gridap.ReferenceFEs: _get_offsets
 using Gridap.ReferenceFEs: _get_offset
 using Gridap.ReferenceFEs: _find_unique_with_indices
 
+using Gridap.Arrays: Reindexed
+
 import Gridap.Arrays: array_cache
 import Gridap.Arrays: getindex!
 import Gridap.Arrays: reindex
 import Gridap.Arrays: get_array
+import Gridap.Arrays: lazy_append
 
 import Gridap.Fields: field_cache
 import Gridap.Fields: evaluate_field!
@@ -202,6 +205,8 @@ export QPointCellField
 
 export RestrictedDiscreteModel
 
+export AppendedTriangulation
+
 include("GridTopologies.jl")
 
 include("GridTopologyMocks.jl")
@@ -252,4 +257,6 @@ include("CellQuadratures.jl")
 
 include("QPointCellFields.jl")
 
+include("AppendedTriangulations.jl")
+
 end # module
diff --git a/src/Geometry/RestrictedTriangulations.jl b/src/Geometry/RestrictedTriangulations.jl
index eb888e3d6..bcc491e45 100644
--- a/src/Geometry/RestrictedTriangulations.jl
+++ b/src/Geometry/RestrictedTriangulations.jl
@@ -54,3 +54,9 @@ end
 function get_cell_id(trian::RestrictedTriangulation)
   trian.cell_to_oldcell
 end
+
+function get_cell_map(trian::RestrictedTriangulation)
+  cell_map = get_cell_map(trian.oldtrian)
+  reindex(cell_map,trian.cell_to_oldcell)
+end
+
diff --git a/src/Geometry/TriangulationPortions.jl b/src/Geometry/TriangulationPortions.jl
index 0fff4675b..ffdf766ca 100644
--- a/src/Geometry/TriangulationPortions.jl
+++ b/src/Geometry/TriangulationPortions.jl
@@ -23,3 +23,7 @@ function get_cell_coordinates(trian::TriangulationPortion)
   reindex(get_cell_coordinates(trian.oldtrian),trian.cell_to_oldcell)
 end
 
+function get_cell_map(trian::TriangulationPortion)
+  cell_map = get_cell_map(trian.oldtrian)
+  reindex(cell_map,trian.cell_to_oldcell)
+end
diff --git a/test/ArraysTests/AppendedArraysTests.jl b/test/ArraysTests/AppendedArraysTests.jl
new file mode 100644
index 000000000..9069ab0b7
--- /dev/null
+++ b/test/ArraysTests/AppendedArraysTests.jl
@@ -0,0 +1,38 @@
+module AppendedArraysTests
+
+using Test
+using Gridap.Arrays
+
+a = collect(11:20)
+b = collect(Float64,101:120)
+c = lazy_append(a,b)
+r = vcat(a,b)
+test_array(c,r)
+
+a = collect(Float64,11:20)
+b = collect(Float64,101:120)
+c = lazy_append(a,b)
+r = vcat(a,b)
+test_array(c,r)
+
+d = apply(-,c)
+r = -c
+test_array(d,r)
+@test isa(d,AppendedArray)
+
+e = apply(-,c,d)
+r = c-d
+test_array(e,r)
+@test isa(e,AppendedArray)
+
+d = apply(Float64,-,c)
+r = -c
+test_array(d,r)
+@test isa(d,AppendedArray)
+
+e = apply(Float64,-,c,d)
+r = c-d
+test_array(e,r)
+@test isa(e,AppendedArray)
+
+end # module
diff --git a/test/ArraysTests/runtests.jl b/test/ArraysTests/runtests.jl
index 76eb61a76..8158eefab 100644
--- a/test/ArraysTests/runtests.jl
+++ b/test/ArraysTests/runtests.jl
@@ -26,4 +26,6 @@ using Test
 
 @testset "ArrayPairs" begin include("ArrayPairsTests.jl") end
 
+@testset "AppendedArrays" begin include("AppendedArraysTests.jl") end
+
 end # module
diff --git a/test/FieldsTests/FieldArraysTests.jl b/test/FieldsTests/FieldArraysTests.jl
index 99485fa33..6d4ba9e07 100644
--- a/test/FieldsTests/FieldArraysTests.jl
+++ b/test/FieldsTests/FieldArraysTests.jl
@@ -32,7 +32,6 @@ agx = fill(gx,l)
 a∇gx = fill(∇gx,l)
 test_array_of_fields(ag,ax,agx,grad=a∇gx)
 
-
 struct FieldPlaceHolder <: Field end
 
 ag = apply_to_field_array(FieldPlaceHolder,bcast(+),af,af)
@@ -91,4 +90,49 @@ agx = fill(r,l)
 a∇gx = fill(∇r,l)
 test_array_of_fields(ag,ax,agx,grad=a∇gx)
 
+# lazy_append
+
+np = 4
+p = Point(1,2)
+x = fill(p,np)
+
+v = 3.0
+d = 2
+f = MockField{d}(v)
+fx = evaluate(f,x)
+∇fx = evaluate(∇(f),x)
+
+l = 10
+af = Fill(f,l)
+ax = fill(x,l)
+afx = fill(fx,l)
+a∇fx = fill(∇fx,l)
+
+np = 5
+p = Point(4,3)
+x = fill(p,np)
+
+v = 5.0
+d = 2
+f = MockField{d}(v)
+fx = evaluate(f,x)
+∇fx = evaluate(∇(f),x)
+
+l = 15
+bf = Fill(f,l)
+bx = fill(x,l)
+bfx = fill(fx,l)
+b∇fx = fill(∇fx,l)
+
+cf = lazy_append(af,bf)
+cx = lazy_append(ax,bx)
+
+cfx = evaluate(cf,cx)
+∇cfx = evaluate(∇(cf),cx)
+rfx = vcat(afx,bfx)
+∇rfx = vcat(a∇fx,b∇fx)
+test_array_of_fields(cf,cx,rfx,grad=∇rfx)
+@test isa(cfx, AppendedArray)
+@test isa(∇cfx, AppendedArray)
+
 end # module
diff --git a/test/GeometryTests/AppendedTriangulationsTests.jl b/test/GeometryTests/AppendedTriangulationsTests.jl
new file mode 100644
index 000000000..aca65b4b4
--- /dev/null
+++ b/test/GeometryTests/AppendedTriangulationsTests.jl
@@ -0,0 +1,80 @@
+module AppendedTriangulationsTests
+
+using Test
+using Gridap.ReferenceFEs
+using Gridap.Arrays
+using Gridap.Geometry
+using Gridap.Visualization
+using Gridap.FESpaces
+using Gridap.Fields
+using Gridap.Integration
+
+domain = (0,1,0,1)
+partition = (10,10)
+model = CartesianDiscreteModel(domain,partition)
+
+ncells = num_cells(model)
+nin = ceil(Int,2*ncells/3)
+cell_to_mask = fill(false,ncells)
+cell_to_mask[1:nin] .= true
+
+grid = get_grid(model)
+
+trian_in = RestrictedTriangulation(grid,cell_to_mask)
+trian_out = RestrictedTriangulation(grid,collect(Bool, .! cell_to_mask))
+
+trian = lazy_append(trian_out,trian_in)
+test_triangulation(trian)
+
+order = 1
+quad = CellQuadrature(trian,2*order)
+quad_in = CellQuadrature(trian_in,2*order)
+quad_out = CellQuadrature(trian_out,2*order)
+
+q = get_coordinates(quad)
+w = get_weights(quad)
+@test isa(q,AppendedArray)
+@test isa(w,AppendedArray)
+
+V = TestFESpace(model=model,valuetype=Float64,order=order,reffe=:Lagrangian,conformity=:H1)
+
+u(x) = x[1]+x[2]
+
+_v = interpolate(V,u)
+v = restrict(_v,trian)
+
+e = u - v
+el2 = sqrt(sum(integrate(e*e,trian,quad)))
+@test el2 < 1.0e-8
+
+_dv = get_cell_basis(V)
+dv = restrict(_dv,trian)
+
+cellmat =  integrate(∇(dv)*∇(dv),trian,quad)
+@test isa(cellmat,AppendedArray)
+@test isa(cellmat.a,CompressedArray)
+@test isa(cellmat.b,CompressedArray)
+
+#writevtk(trian_in,"trian_in")
+#writevtk(trian_out,"trian_out")
+#writevtk(trian,"trian",cellfields=["v"=>v])
+
+# Append triangulations of different cell type
+
+domain = (0,1,0,1)
+partition = (10,10)
+grid1 = CartesianGrid(domain,partition)
+
+domain = (1,2,0,1)
+partition = (10,10)
+grid2 = simplexify(CartesianGrid(domain,partition))
+
+trian = lazy_append(grid1,grid2)
+test_triangulation(trian)
+
+d = mktempdir()
+f = joinpath(d,"trian")
+writevtk(trian,f)
+rm(d,recursive=true)
+
+end # module
diff --git a/test/GeometryTests/runtests.jl b/test/GeometryTests/runtests.jl
index d3c6f50cb..448856d2e 100644
--- a/test/GeometryTests/runtests.jl
+++ b/test/GeometryTests/runtests.jl
@@ -46,4 +46,6 @@ using Test
 
 @testset "CellQuadratures" begin include("CellQuadraturesTests.jl") end
 
+@testset "AppendedTriangulations" begin include("AppendedTriangulationsTests.jl") end
+
 end # module