diff --git a/NEWS.md b/NEWS.md
index c2a3145d..b109d014 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,13 +5,21 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.3.6] 2024-01-28 
+## [0.4.0] 2024-04-12
+
+### Changed
+
+- `DistributedCellField` now inherits from `CellField`. To accomodate the necessary API, we now save a pointer to the `DistributedTriangulation` where it is defined. This also requires `DistributedSingleFieldFESpace` to save the triangulation. Since PR[#141](https://github.com/gridap/GridapDistributed.jl/pull/141).
+- All the distributed `Multifield` cellfield types are now represented by a `DistributedMultiFieldCellField`. Both `DistributedMultiFieldFEFunction` and `DistributedMultiFieldFEBasis` structs have been removed and replaced with constant aliases, which makes it more consistent with single-field types. Since PR[#141](https://github.com/gridap/GridapDistributed.jl/pull/141).
+- Major refactor of ODE module. Implementation has been significantly simplified, while increasing the capability of the API. All `TransientDistributedObjects` structs have been removed, and replaced by `DistributedTransientObjects = DistributedObjects{TransientObject}`. Full support for EX/IM/IMEX methods. See Gridap's release for details. Since PR[#141](https://github.com/gridap/GridapDistributed.jl/pull/141).
+
+## [0.3.6] 2024-01-28
 
 ### Added
 
 - Added redistribution for MultiFieldFESpaces. Since PR [#140](https://github.com/gridap/GridapDistributed.jl/pull/140).
 
-### Fixed 
+### Fixed
 
 - Fixed issue [#142](https://github.com/gridap/GridapDistributed.jl/issues/142). Since PR [#142](https://github.com/gridap/GridapDistributed.jl/pull/142).
 
diff --git a/Project.toml b/Project.toml
index 511ee652..9c0da0a4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GridapDistributed"
 uuid = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 authors = ["S. Badia <santiago.badia@monash.edu>", "A. F. Martin <alberto.f.martin@anu.edu.au>", "F. Verdugo <f.verdugo.rojano@vu.nl>"]
-version = "0.3.6"
+version = "0.4.0"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -17,7 +17,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [compat]
 BlockArrays = "0.16.38"
 FillArrays = "0.8.4,1"
-Gridap = "0.17.23"
+Gridap = "0.18"
 LinearAlgebra = "1.3"
 MPI = "0.16, 0.17, 0.18, 0.19, 0.20"
 PartitionedArrays = "0.3.3"
diff --git a/src/Adaptivity.jl b/src/Adaptivity.jl
index 83313f86..d648260d 100644
--- a/src/Adaptivity.jl
+++ b/src/Adaptivity.jl
@@ -3,9 +3,9 @@
 
 const DistributedAdaptedDiscreteModel{Dc,Dp} = GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:AdaptedDiscreteModel{Dc,Dp}}}
 
-function DistributedAdaptedDiscreteModel(model  ::DistributedDiscreteModel,
-                                         parent ::DistributedDiscreteModel,
-                                         glue   ::AbstractArray{<:AdaptivityGlue})
+function DistributedAdaptedDiscreteModel(model  :: DistributedDiscreteModel,
+                                         parent :: DistributedDiscreteModel,
+                                         glue   :: AbstractArray{<:AdaptivityGlue})
   models = map(local_views(model),local_views(parent),glue) do model, parent, glue
     AdaptedDiscreteModel(model,parent,glue)
   end
diff --git a/src/Algebra.jl b/src/Algebra.jl
index e211bb83..cd1a1634 100644
--- a/src/Algebra.jl
+++ b/src/Algebra.jl
@@ -29,6 +29,57 @@ function Algebra.allocate_in_domain(matrix::BlockPMatrix)
   allocate_in_domain(BlockPVector{V},matrix)
 end
 
+# PSparseMatrix copy
+
+function Base.copy(a::PSparseMatrix)
+  mats = map(copy,partition(a))
+  cache = map(PartitionedArrays.copy_cache,a.cache)
+  return PSparseMatrix(mats,partition(axes(a,1)),partition(axes(a,2)),cache)
+end
+
+# PartitionedArrays extras
+
+function LinearAlgebra.axpy!(α,x::PVector,y::PVector)
+  @check partition(axes(x,1)) === partition(axes(y,1))
+  map(partition(x),partition(y)) do x,y
+    LinearAlgebra.axpy!(α,x,y)
+  end
+  consistent!(y) |> wait
+  return y
+end
+
+function LinearAlgebra.axpy!(α,x::BlockPVector,y::BlockPVector)
+  map(blocks(x),blocks(y)) do x,y
+    LinearAlgebra.axpy!(α,x,y)
+  end
+  return y
+end
+
+function Algebra.axpy_entries!(
+  α::Number, A::PSparseMatrix, B::PSparseMatrix;
+  check::Bool=true
+)
+# We should definitely check here that the index partitions are the same. 
+# However: Because the different matrices are assembled separately, the objects are not the 
+# same (i.e can't use ===). Checking the index partitions would then be costly...
+  @assert reduce(&,map(PartitionedArrays.matching_local_indices,partition(axes(A,1)),partition(axes(B,1))))
+  @assert reduce(&,map(PartitionedArrays.matching_local_indices,partition(axes(A,2)),partition(axes(B,2))))
+  map(partition(A),partition(B)) do A, B
+    Algebra.axpy_entries!(α,A,B;check)
+  end
+  return B
+end
+
+function Algebra.axpy_entries!(
+  α::Number, A::BlockPMatrix, B::BlockPMatrix;
+  check::Bool=true
+)
+  map(blocks(A),blocks(B)) do A, B
+    Algebra.axpy_entries!(α,A,B;check)
+  end
+  return B
+end
+
 # This might go to Gridap in the future. We keep it here for the moment.
 function change_axes(a::Algebra.ArrayCounter,axes)
   @notimplemented
@@ -95,6 +146,9 @@ function num_parts(comm::MPI.Comm)
   end
   nparts
 end
+@inline num_parts(comm::MPIArray) = num_parts(comm.comm)
+@inline num_parts(comm::DebugArray) = length(comm.items)
+@inline num_parts(comm::MPIVoidVector) = num_parts(comm.comm)
 
 function get_part_id(comm::MPI.Comm)
   if comm != MPI.COMM_NULL
@@ -104,23 +158,28 @@ function get_part_id(comm::MPI.Comm)
   end
   id
 end
+@inline get_part_id(comm::MPIArray) = get_part_id(comm.comm)
+@inline get_part_id(comm::MPIVoidVector) = get_part_id(comm.comm)
 
+"""
+    i_am_in(comm::MPIArray)
+    i_am_in(comm::DebugArray)
+  
+  Returns `true` if the processor is part of the subcommunicator `comm`.
+"""
 function i_am_in(comm::MPI.Comm)
   get_part_id(comm) >=0
 end
-
-function i_am_in(comm::MPIArray)
-  i_am_in(comm.comm)
-end
-
-function i_am_in(comm::MPIVoidVector)
-  i_am_in(comm.comm)
-end
+@inline i_am_in(comm::MPIArray) = i_am_in(comm.comm)
+@inline i_am_in(comm::MPIVoidVector) = i_am_in(comm.comm)
+@inline i_am_in(comm::DebugArray) = true
 
 function change_parts(x::Union{MPIArray,DebugArray,Nothing,MPIVoidVector}, new_parts; default=nothing)
-  x_new = map(new_parts) do _p
-    if isa(x,MPIArray) || isa(x,DebugArray)
+  x_new = map(new_parts) do p
+    if isa(x,MPIArray)
       PartitionedArrays.getany(x)
+    elseif isa(x,DebugArray) && (p <= length(x.items))
+      x.items[p]
     else
       default
     end
@@ -314,7 +373,7 @@ end
 
 function local_views(a::BlockPMatrix,new_rows::BlockPRange,new_cols::BlockPRange)
   vals = map(CartesianIndices(blocksize(a))) do I
-    local_views(a[Block(I)],new_rows[Block(I[1])],new_cols[Block(I[2])])
+    local_views(blocks(a)[I],blocks(new_rows)[I],blocks(new_cols)[I])
   end |> to_parray_of_arrays
   return map(mortar,vals)
 end
@@ -795,12 +854,10 @@ function local_views(a::PVectorAllocationTrackTouchedAndValues)
 end
 
 function _setup_prange_from_pvector_allocation(a::PVectorAllocationTrackTouchedAndValues)
-  test_dofs_prange = a.test_dofs_gids_prange # dof ids of the test space
-  ngrdofs = length(test_dofs_prange)
-  
+
   # Find the ghost rows
   allocations = local_views(a.allocations)
-  indices = partition(test_dofs_prange)
+  indices = partition(a.test_dofs_gids_prange)
   I_ghost_lids_to_dofs_ghost_lids = map(allocations, indices) do allocation, indices
     dofs_lids_touched = findall(allocation.touched)
     loc_to_gho = local_to_ghost(indices)
@@ -817,10 +874,11 @@ function _setup_prange_from_pvector_allocation(a::PVectorAllocationTrackTouchedA
     I_ghost_lids
   end
 
-  gids_ghost_to_global, gids_ghost_to_owner = map(
+  ghost_to_global, ghost_to_owner = map(
     find_gid_and_owner,I_ghost_lids_to_dofs_ghost_lids,indices) |> tuple_of_arrays
 
-  _setup_prange_impl_(ngrdofs,indices,gids_ghost_to_global,gids_ghost_to_owner)
+  ngids = length(a.test_dofs_gids_prange)
+  _setup_prange_impl_(ngids,indices,ghost_to_global,ghost_to_owner)
 end
 
 function Algebra.create_from_nz(a::PVectorAllocationTrackTouchedAndValues)
@@ -848,7 +906,7 @@ function find_gid_and_owner(ighost_to_jghost,jindices)
   jghost_to_local  = ghost_to_local(jindices)
   jlocal_to_global = local_to_global(jindices)
   jlocal_to_owner  = local_to_owner(jindices)
-  ighost_to_jlocal = view(jghost_to_local,ighost_to_jghost)
+  ighost_to_jlocal = sort(view(jghost_to_local,ighost_to_jghost))
 
   ighost_to_global = jlocal_to_global[ighost_to_jlocal]
   ighost_to_owner  = jlocal_to_owner[ighost_to_jlocal]
diff --git a/src/BlockPartitionedArrays.jl b/src/BlockPartitionedArrays.jl
index 5661432f..25e43466 100644
--- a/src/BlockPartitionedArrays.jl
+++ b/src/BlockPartitionedArrays.jl
@@ -145,17 +145,19 @@ end
 
 function Base.copyto!(y::BlockPVector,x::BlockPVector)
   @check blocklength(x) == blocklength(y)
-  for i in blockaxes(x,1)
-    copyto!(y[i],x[i])
+  yb, xb = blocks(y), blocks(x)
+  for i in 1:blocksize(x,1)
+    copyto!(yb[i],xb[i])
   end
   return y
 end
 
 function Base.copyto!(y::BlockPMatrix,x::BlockPMatrix)
   @check blocksize(x) == blocksize(y)
-  for i in blockaxes(x,1)
-    for j in blockaxes(x,2)
-      copyto!(y[i,j],x[i,j])
+  yb, xb = blocks(y), blocks(x)
+  for i in 1:blocksize(x,1)
+    for j in 1:blocksize(x,2)
+      copyto!(yb[i,j],xb[i,j])
     end
   end
   return y
@@ -169,6 +171,8 @@ function Base.fill!(a::BlockPVector,v)
 end
 
 function Base.sum(a::BlockPArray)
+  # TODO: This could use a single communication, instead of one for each block
+  # TODO: We could implement a generic reduce, that we apply to sum, all, any, etc..
   return sum(map(sum,blocks(a)))
 end
 
@@ -284,15 +288,47 @@ end
 
 # LinearAlgebra API
 
+function Base.:*(a::Number,b::BlockPArray)
+  mortar(map(bi -> a*bi,blocks(b)))
+end
+Base.:*(b::BlockPMatrix,a::Number) = a*b
+Base.:/(b::BlockPVector,a::Number) = (1/a)*b
+
+function Base.:*(a::BlockPMatrix,b::BlockPVector)
+  c = similar(b)
+  mul!(c,a,b)
+  return c
+end
+
+for op in (:+,:-)
+  @eval begin
+    function Base.$op(a::BlockPArray)
+      mortar(map($op,blocks(a)))
+    end
+    function Base.$op(a::BlockPArray,b::BlockPArray)
+      @assert blocksize(a) == blocksize(b)
+      mortar(map($op,blocks(a),blocks(b)))
+    end
+  end
+end
+
 function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector)
+  o = one(eltype(A))
+  mul!(y,A,x,o,o)
+end
+
+function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector,α::Number,β::Number)
+  yb, Ab, xb = blocks(y), blocks(A), blocks(x)
   z = zero(eltype(y))
   o = one(eltype(A))
-  for i in blockaxes(A,2)
-    fill!(y[i],z)
-    for j in blockaxes(A,2)
-      mul!(y[i],A[i,j],x[j],o,o)
+  for i in 1:blocksize(A,1)
+    fill!(yb[i],z)
+    for j in 1:blocksize(A,2)
+      mul!(yb[i],Ab[i,j],xb[j],α,o)
     end
+    rmul!(yb[i],β)
   end
+  return y
 end
 
 function LinearAlgebra.dot(x::BlockPVector,y::BlockPVector)
@@ -341,16 +377,19 @@ function Base.broadcasted(f, a::Union{BlockPArray,BlockPBroadcasted}, b::Number)
   return BlockPBroadcasted(blocks_out,blockaxes(a))
 end
 
-function Base.broadcasted(f,
-                        a::Union{BlockPArray,BlockPBroadcasted},
-                        b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}})
+function Base.broadcasted(
+  f,
+  a::Union{BlockPArray,BlockPBroadcasted},
+  b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}
+)
   Base.broadcasted(f,a,Base.materialize(b))
 end
 
 function Base.broadcasted(
   f,
   a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}},
-  b::Union{BlockPArray,BlockPBroadcasted})
+  b::Union{BlockPArray,BlockPBroadcasted}
+)
   Base.broadcasted(f,Base.materialize(a),b)
 end
 
diff --git a/src/CellData.jl b/src/CellData.jl
index e0f93c99..0814b398 100644
--- a/src/CellData.jl
+++ b/src/CellData.jl
@@ -1,32 +1,44 @@
 
-abstract type DistributedCellDatum <: GridapType end
-
 # DistributedCellPoint
 """
 """
-struct DistributedCellPoint{A<:AbstractArray{<:CellPoint}} <: DistributedCellDatum
+struct DistributedCellPoint{A<:AbstractArray{<:CellPoint},B<:DistributedTriangulation} <: CellDatum
   points::A
+  trian ::B
 end
 
 local_views(a::DistributedCellPoint) = a.points
+CellData.get_triangulation(a::DistributedCellPoint) = a.trian
+
+function CellData.DomainStyle(::Type{<:DistributedCellPoint{A}}) where A 
+  DomainStyle(eltype(A))
+end
 
 # DistributedCellField
 """
 """
-struct DistributedCellField{A,B} <: DistributedCellDatum
+struct DistributedCellField{A,B,C} <: CellField
   fields::A
-  metadata::B
+  trian ::B
+  metadata::C
   function DistributedCellField(
     fields::AbstractArray{<:CellField},
+    trian ::DistributedTriangulation,
     metadata=nothing)
 
     A = typeof(fields)
-    B = typeof(metadata)
-    new{A,B}(fields,metadata)
+    B = typeof(trian)
+    C = typeof(metadata)
+    new{A,B,C}(fields,trian,metadata)
   end
 end
 
 local_views(a::DistributedCellField) = a.fields
+CellData.get_triangulation(a::DistributedCellField) = a.trian
+
+function CellData.DomainStyle(::Type{<:DistributedCellField{A}}) where A 
+  DomainStyle(eltype(A))
+end
 
 # Constructors
 
@@ -34,14 +46,14 @@ function CellData.CellField(f::Function,trian::DistributedTriangulation)
   fields = map(trian.trians) do t
     CellField(f,t)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,trian)
 end
 
 function CellData.CellField(f::Number,trian::DistributedTriangulation)
   fields = map(trian.trians) do t
     CellField(f,t)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,trian)
 end
 
 function CellData.CellField(
@@ -49,7 +61,7 @@ function CellData.CellField(
   fields = map(f,trian.trians) do f,t
     CellField(f,t)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,trian)
 end
 
 # Evaluation
@@ -59,135 +71,128 @@ function (f::DistributedCellField)(x::DistributedCellPoint)
 end
 
 function Arrays.evaluate!(cache,f::DistributedCellField,x::DistributedCellPoint)
-  map(f.fields,x.points) do f,x
+  map(local_views(f),local_views(x)) do f,x
     evaluate!(nothing,f,x)
   end
 end
 
+# Given local CellFields and a set of original DistributedTriangulations, 
+# returns the DistributedTriangulation where the local CellFields are defined. 
+function _select_triangulation(fields,parents::DistributedCellField...)
+  trian_candidates = unique(objectid,map(get_triangulation,parents))
+  _select_triangulation(fields,trian_candidates...)
+end
+
+function _select_triangulation(fields,trian_candidates::DistributedTriangulation...)
+  if length(trian_candidates) == 1
+    return first(trian_candidates)
+  end
+
+  # Check if we can select one of the original triangulations
+  trians = map(local_views,trian_candidates)
+  t_id = map(fields,trians...) do f, trians...
+    f_id = objectid(get_triangulation(f))
+    return findfirst(tt -> objectid(tt) == f_id, trians)
+  end |> getany
+  if !isnothing(t_id)
+    return trian_candidates[t_id]
+  end
+
+  # If not, check if we can build a new DistributedTriangulation based on one of the original models. 
+  m_id = map(fields,trians...) do f, trians...
+    f_id = objectid(get_background_model(get_triangulation(f)))
+    return findfirst(tt -> objectid(get_background_model(tt)) == f_id, trians)
+  end |> getany
+  if !isnothing(m_id)
+    model = get_background_model(trian_candidates[m_id])
+    return DistributedTriangulation(map(get_triangulation,fields),model)
+  end
+
+  @error "Cannot select a triangulation for the operation"
+end
+
 # Operations
 
 function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField)
-  fields = map(a.fields) do f
+  fields = map(local_views(a)) do f
     evaluate!(nothing,k,f)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,get_triangulation(a))
 end
 
 function Arrays.evaluate!(
   cache,k::Operation,a::DistributedCellField,b::DistributedCellField)
-  fields = map(a.fields,b.fields) do f,g
+  fields = map(local_views(a),local_views(b)) do f,g
     evaluate!(nothing,k,f,g)
   end
-  DistributedCellField(fields)
+  trian = _select_triangulation(fields,a,b)
+  DistributedCellField(fields,trian)
 end
 
 function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::Number)
-  fields = map(a.fields) do f
+  fields = map(local_views(a)) do f
     evaluate!(nothing,k,f,b)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,get_triangulation(a))
 end
 
 function Arrays.evaluate!(cache,k::Operation,b::Number,a::DistributedCellField)
-  fields = map(a.fields) do f
+  fields = map(local_views(a)) do f
     evaluate!(nothing,k,b,f)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,get_triangulation(a))
 end
 
 function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::Function)
-  fields = map(a.fields) do f
+  fields = map(local_views(a)) do f
     evaluate!(nothing,k,f,b)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,get_triangulation(a))
 end
 
 function Arrays.evaluate!(cache,k::Operation,b::Function,a::DistributedCellField)
-  fields = map(a.fields) do f
+  fields = map(local_views(a)) do f
     evaluate!(nothing,k,b,f)
   end
-  DistributedCellField(fields)
+  DistributedCellField(fields,get_triangulation(a))
 end
 
 function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField...)
-  fields = map(map(i->i.fields,a)...) do f...
+  fields = map(map(local_views,a)...) do f...
     evaluate!(nothing,k,f...)
   end
-  DistributedCellField(fields)
+  trian = _select_triangulation(fields,a...)
+  DistributedCellField(fields,trian)
 end
 
-# Composition
-
-Base.:(∘)(f::Function,g::DistributedCellField) = Operation(f)(g)
-Base.:(∘)(f::Function,g::Tuple{DistributedCellField,DistributedCellField}) = Operation(f)(g[1],g[2])
-Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Number}) = Operation(f)(g[1],g[2])
-Base.:(∘)(f::Function,g::Tuple{Number,DistributedCellField}) = Operation(f)(g[1],g[2])
-Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Function}) = Operation(f)(g[1],g[2])
-Base.:(∘)(f::Function,g::Tuple{Function,DistributedCellField}) = Operation(f)(g[1],g[2])
-Base.:(∘)(f::Function,g::Tuple{Vararg{DistributedCellField}}) = Operation(f)(g...)
-
-# Define some of the well known arithmetic ops
-
-# Unary ops
-
-for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj)
-  @eval begin
-    ($op)(a::DistributedCellField) = Operation($op)(a)
-  end
-end
-
-# Binary ops
-
-for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/)
-  @eval begin
-    ($op)(a::DistributedCellField,b::DistributedCellField) = Operation($op)(a,b)
-    ($op)(a::DistributedCellField,b::Number) = Operation($op)(a,b)
-    ($op)(a::Number,b::DistributedCellField) = Operation($op)(a,b)
-    ($op)(a::DistributedCellField,b::Function) = Operation($op)(a,b)
-    ($op)(a::Function,b::DistributedCellField) = Operation($op)(a,b)
-  end
-end
-
-Base.broadcasted(f,a::DistributedCellField,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b)
-Base.broadcasted(f,a::Number,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b)
-Base.broadcasted(f,a::DistributedCellField,b::Number) = Operation((i,j)->f.(i,j))(a,b)
-Base.broadcasted(f,a::Function,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b)
-Base.broadcasted(f,a::DistributedCellField,b::Function) = Operation((i,j)->f.(i,j))(a,b)
-Base.broadcasted(::typeof(*),::typeof(∇),f::DistributedCellField) = Operation(Fields._extract_grad_diag)(∇(f))
-Base.broadcasted(::typeof(*),s::Fields.ShiftedNabla,f::DistributedCellField) = Operation(Fields._extract_grad_diag)(s(f))
-
-dot(::typeof(∇),f::DistributedCellField) = divergence(f)
-outer(::typeof(∇),f::DistributedCellField) = gradient(f)
-outer(f::DistributedCellField,::typeof(∇)) = transpose(gradient(f))
-cross(::typeof(∇),f::DistributedCellField) = curl(f)
-
 # Differential ops
 
 function Fields.gradient(a::DistributedCellField)
-  DistributedCellField(map(gradient,a.fields))
+  DistributedCellField(map(gradient,a.fields),get_triangulation(a))
 end
 
 function Fields.divergence(a::DistributedCellField)
-  DistributedCellField(map(divergence,a.fields))
+  DistributedCellField(map(divergence,a.fields),get_triangulation(a))
 end
 
 function Fields.DIV(a::DistributedCellField)
-  DistributedCellField(map(DIV,a.fields))
+  DistributedCellField(map(DIV,a.fields),get_triangulation(a))
 end
 
 function Fields.∇∇(a::DistributedCellField)
-  DistributedCellField(map(∇∇,a.fields))
+  DistributedCellField(map(∇∇,a.fields),get_triangulation(a))
 end
 
 function Fields.curl(a::DistributedCellField)
-  DistributedCellField(map(curl,a.fields))
+  DistributedCellField(map(curl,a.fields),get_triangulation(a))
 end
 
 # Integration related
 """
 """
-struct DistributedMeasure{A<:AbstractArray{<:Measure}} <: GridapType
+struct DistributedMeasure{A<:AbstractArray{<:Measure},B<:DistributedTriangulation} <: GridapType
   measures::A
+  trian::B
 end
 
 local_views(a::DistributedMeasure) = a.measures
@@ -196,18 +201,18 @@ function CellData.Measure(t::DistributedTriangulation,args...)
   measures = map(t.trians) do trian
     Measure(trian,args...)
   end
-  DistributedMeasure(measures)
+  DistributedMeasure(measures,t)
 end
 
 function CellData.Measure(tt::DistributedTriangulation{Dc,Dp},it::DistributedTriangulation{Dc,Dp},args...) where {Dc,Dp}
   measures = map(local_views(tt),local_views(it)) do ttrian, itrian
     Measure(ttrian,itrian,args...)
   end
-  return DistributedMeasure(measures)
+  return DistributedMeasure(measures,it)
 end
 
 function CellData.get_cell_points(a::DistributedMeasure)
-  DistributedCellPoint(map(get_cell_points,a.measures))
+  DistributedCellPoint(map(get_cell_points,a.measures),a.trian)
 end
 
 """
@@ -216,6 +221,8 @@ struct DistributedDomainContribution{A<:AbstractArray{<:DomainContribution}} <:
   contribs::A
 end
 
+CellData.num_domains(a::DistributedDomainContribution) = CellData.num_domains(getany(local_views(a)))
+
 local_views(a::DistributedDomainContribution) = a.contribs
 
 function Base.getindex(c::DistributedDomainContribution,t::DistributedTriangulation)
@@ -272,89 +279,77 @@ end
 
 (*)(a::DistributedDomainContribution,b::Number) = b*a
 
+# Jordi: This is ugly, but it is useful to re-use code from Gridap: 
+# A lot of the time, we create an empty DomainContribution and then add to it.
+# By dispatching here, this kind of code works verbatim for GridapDistributed. 
+# We could eventually replace this with an EmptyDomainContribution type.
+function (+)(a::CellData.DomainContribution,b::DistributedDomainContribution)
+  @assert iszero(CellData.num_domains(a))
+  return b
+end
+
 # Triangulation related
 
 function CellData.get_cell_points(a::DistributedTriangulation)
-  DistributedCellPoint(map(get_cell_points,a.trians))
+  DistributedCellPoint(map(get_cell_points,a.trians),a)
 end
 
 function CellData.get_normal_vector(a::DistributedTriangulation)
   fields = map(get_normal_vector,a.trians)
-  DistributedCellField(fields)
+  DistributedCellField(fields,a)
 end
 
 # Skeleton related
 
-function DistributedCellField(a::AbstractArray{<:SkeletonPair})
+function DistributedCellField(a::AbstractArray{<:SkeletonPair},trian::DistributedTriangulation)
   plus, minus = map(s->(s.plus,s.minus),a) |> tuple_of_arrays
-  dplus = DistributedCellField(plus)
-  dminus = DistributedCellField(minus)
+  tplus  = _select_triangulation(plus,trian)
+  tminus = _select_triangulation(minus,trian)
+  dplus  = DistributedCellField(plus,tplus)
+  dminus = DistributedCellField(minus,tminus)
   SkeletonPair(dplus,dminus)
 end
 
 function Base.getproperty(x::DistributedCellField, sym::Symbol)
   if sym in (:⁺,:plus)
-    DistributedCellField(map(i->i.plus,x.fields))
+    DistributedCellField(map(i->i.plus,local_views(x)),get_triangulation(x))
   elseif sym in (:⁻, :minus)
-    DistributedCellField(map(i->i.minus,x.fields))
+    DistributedCellField(map(i->i.minus,local_views(x)),get_triangulation(x))
   else
     getfield(x, sym)
   end
 end
 
-function Base.propertynames(x::DistributedCellField, private::Bool=false)
-  (fieldnames(typeof(x))...,:⁺,:plus,:⁻,:minus)
-end
-
-for op in (:outer,:*,:dot)
-  @eval begin
-    ($op)(a::DistributedCellField,b::SkeletonPair{<:DistributedCellField}) = Operation($op)(a,b)
-    ($op)(a::SkeletonPair{<:DistributedCellField},b::DistributedCellField) = Operation($op)(a,b)
-  end
-end
-
-function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::SkeletonPair{<:DistributedCellField})
-  plus = k(a.plus,b.plus)
-  minus = k(a.minus,b.minus)
-  SkeletonPair(plus,minus)
-end
-
-function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::DistributedCellField)
-  plus = k(a.plus,b.plus)
-  minus = k(a.minus,b.minus)
-  SkeletonPair(plus,minus)
-end
-
-CellData.jump(a::DistributedCellField) = DistributedCellField(map(jump,a.fields))
-CellData.jump(a::SkeletonPair{<:DistributedCellField}) = a.⁺ + a.⁻
-CellData.mean(a::DistributedCellField) = DistributedCellField(map(mean,a.fields))
-
+CellData.jump(a::DistributedCellField) = DistributedCellField(map(jump,a.fields),get_triangulation(a))
+CellData.mean(a::DistributedCellField) = DistributedCellField(map(mean,a.fields),get_triangulation(a))
 
 # DistributedCellDof
 
-struct DistributedCellDof{A} <: DistributedCellDatum
+struct DistributedCellDof{A<:AbstractArray{<:CellDof},B<:DistributedTriangulation} <: CellDatum
   dofs::A
-  function DistributedCellDof(dofs::AbstractArray{<:CellDof})
-      A = typeof(dofs)
-      new{A}(dofs)
-  end
+  trian::B
 end
 
 local_views(s::DistributedCellDof) = s.dofs
+CellData.get_triangulation(s::DistributedCellDof) = s.trian
+
+function CellData.DomainStyle(::Type{<:DistributedCellDof{A}}) where A 
+  DomainStyle(eltype(A))
+end
 
 (a::DistributedCellDof)(f) = evaluate(a,f)
 
 function Gridap.Arrays.evaluate!(cache,s::DistributedCellDof,f::DistributedCellField)
   map(local_views(s),local_views(f)) do s, f
-      evaluate!(nothing,s,f)
+    evaluate!(nothing,s,f)
   end
 end
 
 function Gridap.Arrays.evaluate!(cache, ::DistributedCellField, ::DistributedCellDof)
-@unreachable """\n
-CellField (f) objects cannot be evaluated at CellDof (s) objects.
-However, CellDofs objects can be evaluated at CellField objects.
-Did you mean evaluate(f,s) instead of evaluate(s,f), i.e.
-f(s) instead of s(f)?
-"""
+  @unreachable """\n
+  CellField (f) objects cannot be evaluated at CellDof (s) objects.
+  However, CellDofs objects can be evaluated at CellField objects.
+  Did you mean evaluate(f,s) instead of evaluate(s,f), i.e.
+  f(s) instead of s(f)?
+  """
 end
\ No newline at end of file
diff --git a/src/DivConformingFESpaces.jl b/src/DivConformingFESpaces.jl
index ac31e505..c814e75e 100644
--- a/src/DivConformingFESpaces.jl
+++ b/src/DivConformingFESpaces.jl
@@ -27,8 +27,9 @@ function _common_fe_space_constructor(model,cell_reffes;conformity,kwargs...)
      FESpace(m, cell_fe; kwargs...)
   end
   gids = generate_gids(model,spaces)
+  trian = DistributedTriangulation(map(get_triangulation,spaces),model)
   vector_type = _find_vector_type(spaces,gids)
-  DistributedSingleFieldFESpace(spaces,gids,vector_type)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
 end
 
 function _generate_sign_flips(model,cell_reffes)
diff --git a/src/FESpaces.jl b/src/FESpaces.jl
index 8972ff49..4480c2f4 100644
--- a/src/FESpaces.jl
+++ b/src/FESpaces.jl
@@ -130,8 +130,6 @@ function generate_gids(
   cell_to_ldofs::AbstractArray{<:AbstractArray},
   nldofs::AbstractArray{<:Integer})
 
-  ngcells = length(cell_range)
-
   # Find and count number owned dofs
   ldof_to_owner, nodofs = map(partition(cell_range),cell_to_ldofs,nldofs) do indices,cell_to_ldofs,nldofs
     ldof_to_owner = fill(Int32(0),nldofs)
@@ -260,7 +258,7 @@ end
 
 """
 """
-const DistributedSingleFieldFEFunction = DistributedCellField{A,<:DistributedFEFunctionData{T}} where {A,T}
+const DistributedSingleFieldFEFunction{A,B,T} = DistributedCellField{A,B,DistributedFEFunctionData{T}}
 
 function FESpaces.get_free_dof_values(uh::DistributedSingleFieldFEFunction)
   uh.metadata.free_values
@@ -269,21 +267,25 @@ end
 # Single field related
 """
 """
-struct DistributedSingleFieldFESpace{A,B,C} <: DistributedFESpace
+struct DistributedSingleFieldFESpace{A,B,C,D} <: DistributedFESpace
   spaces::A
   gids::B
-  vector_type::Type{C}
+  trian::C
+  vector_type::Type{D}
   function DistributedSingleFieldFESpace(
     spaces::AbstractArray{<:SingleFieldFESpace},
     gids::PRange,
-    vector_type::Type{C}) where C
+    trian::DistributedTriangulation,
+    vector_type::Type{D}) where D
     A = typeof(spaces)
     B = typeof(gids)
-    new{A,B,C}(spaces,gids,vector_type)
+    C = typeof(trian)
+    new{A,B,C,D}(spaces,gids,trian,vector_type)
   end
 end
 
 local_views(a::DistributedSingleFieldFESpace) = a.spaces
+CellData.get_triangulation(a::DistributedSingleFieldFESpace) = a.trian
 
 function FESpaces.get_vector_type(fs::DistributedSingleFieldFESpace)
   fs.vector_type
@@ -338,8 +340,9 @@ function _EvaluationFunction(func,
   f::DistributedSingleFieldFESpace,x::AbstractVector,isconsistent=false)
   free_values = change_ghost(x,f.gids,is_consistent=isconsistent,make_consistent=true)
   fields   = map(func,f.spaces,partition(free_values))
+  trian    = get_triangulation(f)
   metadata = DistributedFEFunctionData(free_values)
-  DistributedCellField(fields,metadata)
+  DistributedCellField(fields,trian,metadata)
 end
 
 function _EvaluationFunction(func,
@@ -347,56 +350,60 @@ function _EvaluationFunction(func,
   dirichlet_values::AbstractArray{<:AbstractVector},isconsistent=false)
   free_values = change_ghost(x,f.gids,is_consistent=isconsistent,make_consistent=true)
   fields   = map(func,f.spaces,partition(free_values),dirichlet_values)
+  trian    = get_triangulation(f)
   metadata = DistributedFEFunctionData(free_values)
-  DistributedCellField(fields,metadata)
+  DistributedCellField(fields,trian,metadata)
 end
 
 function FESpaces.get_fe_basis(f::DistributedSingleFieldFESpace)
   fields = map(get_fe_basis,f.spaces)
-  DistributedCellField(fields)
+  trian  = get_triangulation(f)
+  DistributedCellField(fields,trian)
 end
 
 function FESpaces.get_trial_fe_basis(f::DistributedSingleFieldFESpace)
   fields = map(get_trial_fe_basis,f.spaces)
-  DistributedCellField(fields)
+  trian  = get_triangulation(f)
+  DistributedCellField(fields,trian)
 end
 
 function FESpaces.get_fe_dof_basis(f::DistributedSingleFieldFESpace)
-  dofs = map(get_fe_dof_basis,local_views(f))
-  DistributedCellDof(dofs)
+  dofs  = map(get_fe_dof_basis,local_views(f))
+  trian = get_triangulation(f)
+  DistributedCellDof(dofs,trian)
 end
 
 function FESpaces.TrialFESpace(f::DistributedSingleFieldFESpace)
   spaces = map(TrialFESpace,f.spaces)
-  DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type)
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type)
 end
 
 function FESpaces.TrialFESpace(f::DistributedSingleFieldFESpace,fun)
   spaces = map(f.spaces) do s
     TrialFESpace(s,fun)
   end
-  DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type)
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type)
 end
 
 function FESpaces.TrialFESpace(fun,f::DistributedSingleFieldFESpace)
   spaces = map(f.spaces) do s
     TrialFESpace(fun,s)
   end
-  DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type)
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type)
 end
 
 function FESpaces.TrialFESpace!(f::DistributedSingleFieldFESpace,fun)
   spaces = map(f.spaces) do s
     TrialFESpace!(s,fun)
   end
-  DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type)
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type)
 end
 
 function FESpaces.HomogeneousTrialFESpace(f::DistributedSingleFieldFESpace)
   spaces = map(f.spaces) do s
     HomogeneousTrialFESpace(s)
   end
-  DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type)
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type)
 end
 
 function generate_gids(
@@ -478,8 +485,9 @@ function FESpaces.FESpace(model::DistributedDiscreteModel,reffe;split_own_and_gh
     FESpace(m,reffe;kwargs...)
   end
   gids =  generate_gids(model,spaces)
+  trian = DistributedTriangulation(map(get_triangulation,spaces),model)
   vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
-  DistributedSingleFieldFESpace(spaces,gids,vector_type)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
 end
 
 function FESpaces.FESpace(_trian::DistributedTriangulation,reffe;split_own_and_ghost=false,kwargs...)
@@ -492,7 +500,7 @@ function FESpaces.FESpace(_trian::DistributedTriangulation,reffe;split_own_and_g
   nldofs = map(num_free_dofs,spaces)
   gids = generate_gids(trian_gids,cell_to_ldofs,nldofs)
   vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
-  DistributedSingleFieldFESpace(spaces,gids,vector_type)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
 end
 
 function _find_vector_type(spaces,gids;split_own_and_ghost=false)
diff --git a/src/Geometry.jl b/src/Geometry.jl
index d80d9cb9..3a3b3568 100644
--- a/src/Geometry.jl
+++ b/src/Geometry.jl
@@ -398,6 +398,10 @@ function Geometry.get_background_model(a::DistributedTriangulation)
   a.model
 end
 
+function Geometry.num_cells(a::DistributedTriangulation)
+  sum(map(trian->num_cells(trian),local_views(a)))
+end
+
 # Triangulation constructors
 
 function Geometry.Triangulation(
diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl
index 27f4fc02..0472ef61 100644
--- a/src/GridapDistributed.jl
+++ b/src/GridapDistributed.jl
@@ -12,12 +12,12 @@ using Gridap.CellData
 using Gridap.Visualization
 using Gridap.FESpaces
 using Gridap.MultiField
-using Gridap.ODEs.TransientFETools
-using Gridap.ODEs.ODETools
+using Gridap.ODEs
 
 using MPI
 using PartitionedArrays
 const PArrays = PartitionedArrays
+using PartitionedArrays: getany
 
 using SparseArrays
 using WriteVTK
@@ -29,7 +29,6 @@ import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
 import LinearAlgebra: det, tr, cross, dot, ⋅, diag
 import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj, getproperty, propertynames
 import Gridap.Fields: grad2curl
-import Gridap.ODEs.ODETools: ∂t, ∂tt
 
 export FullyAssembledRows
 export SubAssembledRows
@@ -58,11 +57,7 @@ include("DivConformingFESpaces.jl")
 
 include("MultiField.jl")
 
-include("TransientDistributedCellField.jl")
-
-include("TransientMultiFieldDistributedCellField.jl")
-
-include("TransientFESpaces.jl")
+include("ODEs.jl")
 
 include("Adaptivity.jl")
 
diff --git a/src/MultiField.jl b/src/MultiField.jl
index 2d535ab4..e8927d92 100644
--- a/src/MultiField.jl
+++ b/src/MultiField.jl
@@ -1,31 +1,60 @@
 
-struct DistributedMultiFieldFEFunction{A,B,C} <: GridapType
+# DistributedMultiFieldCellField
+struct DistributedMultiFieldCellField{A,B,C} <: CellField
   field_fe_fun::A
   part_fe_fun::B
-  free_values::C
-  function DistributedMultiFieldFEFunction(
-    field_fe_fun::AbstractVector{<:DistributedSingleFieldFEFunction},
-    part_fe_fun::AbstractArray{<:MultiFieldFEFunction},
-    free_values::AbstractVector)
+  metadata::C
+  function DistributedMultiFieldCellField(
+    field_fe_fun::AbstractVector{<:DistributedCellField},
+    part_fe_fun::AbstractArray{<:CellField},
+    metadata=nothing)
     A = typeof(field_fe_fun)
     B = typeof(part_fe_fun)
-    C = typeof(free_values)
-    new{A,B,C}(field_fe_fun,part_fe_fun,free_values)
+    C = typeof(metadata)
+    new{A,B,C}(field_fe_fun,part_fe_fun,metadata)
   end
 end
 
+function CellData.get_triangulation(a::DistributedMultiFieldCellField)
+  trians = map(get_triangulation,a.field_fe_fun)
+  @check all(map(t -> t === first(trians), trians))
+  return first(trians)
+end
 
-function FESpaces.get_free_dof_values(uh::DistributedMultiFieldFEFunction)
-  uh.free_values
+function CellData.DomainStyle(::Type{<:DistributedMultiFieldCellField{A,B}}) where {A,B}
+  DomainStyle(eltype(B))
+end
+
+local_views(a::DistributedMultiFieldCellField) = a.part_fe_fun
+local_views(a::Vector{<:DistributedCellField}) = map(local_views,a)
+
+MultiField.num_fields(m::DistributedMultiFieldCellField) = length(m.field_fe_fun)
+Base.iterate(m::DistributedMultiFieldCellField) = iterate(m.field_fe_fun)
+Base.iterate(m::DistributedMultiFieldCellField,state) = iterate(m.field_fe_fun,state)
+Base.getindex(m::DistributedMultiFieldCellField,field_id::Integer) = m.field_fe_fun[field_id]
+
+function LinearAlgebra.dot(a::DistributedMultiFieldCellField,b::DistributedMultiFieldCellField)
+  @check num_fields(a) == num_fields(b)
+  return sum(map(dot,a.field_fe_fun,b.field_fe_fun))
 end
 
-local_views(a::DistributedMultiFieldFEFunction) = a.part_fe_fun
-MultiField.num_fields(m::DistributedMultiFieldFEFunction) = length(m.field_fe_fun)
-Base.iterate(m::DistributedMultiFieldFEFunction) = iterate(m.field_fe_fun)
-Base.iterate(m::DistributedMultiFieldFEFunction,state) = iterate(m.field_fe_fun,state)
-Base.getindex(m::DistributedMultiFieldFEFunction,field_id::Integer) = m.field_fe_fun[field_id]
+# DistributedMultiFieldFEFunction
+
+const DistributedMultiFieldFEFunction{A,B,T} = DistributedMultiFieldCellField{A,B,DistributedFEFunctionData{T}}
 
-local_views(a::Vector{<:DistributedCellField}) = [ai.fields for ai in a]
+function DistributedMultiFieldFEFunction(
+  field_fe_fun::AbstractVector{<:DistributedSingleFieldFEFunction},
+  part_fe_fun::AbstractArray{<:MultiFieldFEFunction},
+  free_values::AbstractVector)
+  metadata = DistributedFEFunctionData(free_values)
+  DistributedMultiFieldCellField(field_fe_fun,part_fe_fun,metadata)
+end
+
+function FESpaces.get_free_dof_values(uh::DistributedMultiFieldFEFunction)
+  uh.metadata.free_values
+end
+
+# DistributedMultiFieldFESpace
 
 """
 """
@@ -47,6 +76,15 @@ struct DistributedMultiFieldFESpace{MS,A,B,C,D} <: DistributedFESpace
   end
 end
 
+function CellData.get_triangulation(a::DistributedMultiFieldFESpace)
+  trians = map(get_triangulation,a.field_fe_space)
+  @check all(map(t -> t === first(trians), trians))
+  return first(trians)
+end
+
+MultiField.MultiFieldStyle(::Type{<:DistributedMultiFieldFESpace{MS}}) where MS = MS()
+MultiField.MultiFieldStyle(a::DistributedMultiFieldFESpace) = MultiField.MultiFieldStyle(typeof(a))
+
 local_views(a::DistributedMultiFieldFESpace) = a.part_fe_space
 MultiField.num_fields(m::DistributedMultiFieldFESpace) = length(m.field_fe_space)
 Base.iterate(m::DistributedMultiFieldFESpace) = iterate(m.field_fe_space)
@@ -63,8 +101,9 @@ function FESpaces.get_free_dof_ids(fs::DistributedMultiFieldFESpace)
 end
 
 function MultiField.restrict_to_field(
-  f::DistributedMultiFieldFESpace,free_values::AbstractVector,field::Integer)
-  values = map(f.part_fe_space,partition(free_values)) do u,x
+  f::DistributedMultiFieldFESpace,free_values::AbstractVector,field::Integer
+)
+  values = map(local_views(f),partition(free_values)) do u,x
     restrict_to_field(u,x,field)
   end
   gids = f.field_fe_space[field].gids
@@ -75,8 +114,11 @@ function FESpaces.zero_dirichlet_values(f::DistributedMultiFieldFESpace)
   map(zero_dirichlet_values,f.field_fe_space)
 end
 
-function FESpaces.FEFunction(
-  f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false)
+function FESpaces.get_dirichlet_dof_values(f::DistributedMultiFieldFESpace)
+  return map(get_dirichlet_dof_values,f.field_fe_space)
+end
+
+function FESpaces.FEFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false)
   free_values  = change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true)
   part_fe_fun  = map(FEFunction,f.part_fe_space,partition(free_values))
   field_fe_fun = DistributedSingleFieldFEFunction[]
@@ -89,8 +131,25 @@ function FESpaces.FEFunction(
   DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
 end
 
-function FESpaces.EvaluationFunction(
-  f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false)
+function FESpaces.FEFunction(
+  f::DistributedMultiFieldFESpace,x::AbstractVector,
+  dirichlet_values::AbstractArray{<:AbstractVector},isconsistent=false
+)
+  free_values  = GridapDistributed.change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true)
+  part_dirvals = to_parray_of_arrays(dirichlet_values)
+  part_fe_fun  = map(FEFunction,f.part_fe_space,partition(free_values),part_dirvals)
+  field_fe_fun = DistributedSingleFieldFEFunction[]
+  for i in 1:num_fields(f)
+    free_values_i = restrict_to_field(f,free_values,i)
+    dirichlet_values_i = dirichlet_values[i]
+    fe_space_i = f.field_fe_space[i]
+    fe_fun_i = FEFunction(fe_space_i,free_values_i,dirichlet_values_i,true)
+    push!(field_fe_fun,fe_fun_i)
+  end
+  DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
+end
+
+function FESpaces.EvaluationFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false)
   free_values  = change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true)
   part_fe_fun  = map(EvaluationFunction,f.part_fe_space,partition(free_values))
   field_fe_fun = DistributedSingleFieldFEFunction[]
@@ -122,6 +181,24 @@ function FESpaces.interpolate!(objects,free_values::AbstractVector,fe::Distribut
   DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
 end
 
+function Gridap.FESpaces.interpolate!(
+  objects::Union{<:DistributedMultiFieldCellField,<:DistributedCellField},
+  free_values::AbstractVector,
+  fe::DistributedMultiFieldFESpace
+)
+  part_fe_fun = map(local_views(objects),partition(free_values),local_views(fe)) do objects,x,f
+    interpolate!(objects,x,f)
+  end
+  field_fe_fun = GridapDistributed.DistributedSingleFieldFEFunction[]
+  for i in 1:num_fields(fe)
+    free_values_i = Gridap.MultiField.restrict_to_field(fe,free_values,i)
+    fe_space_i = fe.field_fe_space[i]
+    fe_fun_i = FEFunction(fe_space_i,free_values_i)
+    push!(field_fe_fun,fe_fun_i)
+  end
+  GridapDistributed.DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
+end
+
 function FESpaces.interpolate_everywhere(objects,fe::DistributedMultiFieldFESpace)
   free_values = zero_free_values(fe)
   part_fe_fun = map(partition(free_values),local_views(fe)) do x,f
@@ -178,49 +255,6 @@ function FESpaces.interpolate_everywhere(
   DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
 end
 
-
-"""
-"""
-struct DistributedMultiFieldFEBasis{A,B} <: GridapType
-  field_fe_basis::A
-  part_fe_basis::B
-  function DistributedMultiFieldFEBasis(
-    field_fe_basis::AbstractVector{<:DistributedCellField},
-    part_fe_basis::AbstractArray{<:MultiFieldCellField})
-    A = typeof(field_fe_basis)
-    B = typeof(part_fe_basis)
-    new{A,B}(field_fe_basis,part_fe_basis)
-  end
-end
-
-local_views(a::DistributedMultiFieldFEBasis) = a.part_fe_basis
-MultiField.num_fields(m::DistributedMultiFieldFEBasis) = length(m.field_fe_basis)
-Base.iterate(m::DistributedMultiFieldFEBasis) = iterate(m.field_fe_basis)
-Base.iterate(m::DistributedMultiFieldFEBasis,state) = iterate(m.field_fe_basis,state)
-Base.getindex(m::DistributedMultiFieldFEBasis,field_id::Integer) = m.field_fe_basis[field_id]
-
-function FESpaces.get_fe_basis(f::DistributedMultiFieldFESpace)
-  part_mbasis = map(get_fe_basis,f.part_fe_space)
-  field_fe_basis = DistributedCellField[]
-  for i in 1:num_fields(f)
-    basis_i = map(b->b[i],part_mbasis)
-    bi = DistributedCellField(basis_i)
-    push!(field_fe_basis,bi)
-  end
-  DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis)
-end
-
-function FESpaces.get_trial_fe_basis(f::DistributedMultiFieldFESpace)
-  part_mbasis = map(get_trial_fe_basis,f.part_fe_space)
-  field_fe_basis = DistributedCellField[]
-  for i in 1:num_fields(f)
-    basis_i = map(b->b[i],part_mbasis)
-    bi = DistributedCellField(basis_i)
-    push!(field_fe_basis,bi)
-  end
-  DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis)
-end
-
 function FESpaces.TrialFESpace(objects,a::DistributedMultiFieldFESpace)
   TrialFESpace(a,objects)
 end
@@ -237,6 +271,30 @@ function FESpaces.TrialFESpace(a::DistributedMultiFieldFESpace{MS},objects) wher
   DistributedMultiFieldFESpace(f_dspace,p_mspace,gids,vector_type)
 end
 
+# DistributedMultiFieldFEBasis
+
+const DistributedMultiFieldFEBasis{A} = DistributedMultiFieldCellField{A,<:AbstractArray{<:FEBasis}}
+
+function FESpaces.get_fe_basis(f::DistributedMultiFieldFESpace)
+  part_mbasis = map(get_fe_basis,f.part_fe_space)
+  field_fe_basis = map(1:num_fields(f)) do i
+    space_i  = f.field_fe_space[i]
+    basis_i = map(b->b[i],part_mbasis)
+    DistributedCellField(basis_i,get_triangulation(space_i))
+  end
+  DistributedMultiFieldCellField(field_fe_basis,part_mbasis)
+end
+
+function FESpaces.get_trial_fe_basis(f::DistributedMultiFieldFESpace)
+  part_mbasis = map(get_trial_fe_basis,f.part_fe_space)
+  field_fe_basis = map(1:num_fields(f)) do i
+    space_i  = f.field_fe_space[i]
+    basis_i = map(b->b[i],part_mbasis)
+    DistributedCellField(basis_i,get_triangulation(space_i))
+  end
+  DistributedMultiFieldCellField(field_fe_basis,part_mbasis)
+end
+
 # Factory
 
 function MultiField.MultiFieldFESpace(
diff --git a/src/ODEs.jl b/src/ODEs.jl
new file mode 100644
index 00000000..1073ea4b
--- /dev/null
+++ b/src/ODEs.jl
@@ -0,0 +1,163 @@
+
+# Distributed FESpace commons
+
+function Arrays.evaluate!(transient_space::DistributedFESpace, space::DistributedFESpace, t::Real)
+  map(local_views(transient_space),local_views(space)) do transient_space, space
+    Arrays.evaluate!(transient_space,space,t)
+  end
+  return transient_space
+end
+
+# SingleField FESpace
+
+const DistributedTransientTrialFESpace = DistributedSingleFieldFESpace{<:AbstractArray{<:ODEs.TransientTrialFESpace}}
+
+function ODEs.TransientTrialFESpace(space::DistributedSingleFieldFESpace,transient_dirichlet::Union{Function, AbstractVector{<:Function}})
+  spaces = map(local_views(space)) do space
+    ODEs.TransientTrialFESpace(space,transient_dirichlet)
+  end
+  gids  = get_free_dof_ids(space)
+  trian = get_triangulation(space)
+  vector_type = get_vector_type(space)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+end
+
+function ODEs.TransientTrialFESpace(space::DistributedSingleFieldFESpace)
+  spaces = map(local_views(space)) do space
+    ODEs.TransientTrialFESpace(space)
+  end
+  gids  = get_free_dof_ids(space)
+  trian = get_triangulation(space)
+  vector_type = get_vector_type(space)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+end
+
+function ODEs.allocate_space(space::DistributedTransientTrialFESpace)
+  spaces = map(local_views(space)) do space
+    ODEs.allocate_space(space)
+  end
+  gids  = get_free_dof_ids(space)
+  trian = get_triangulation(space)
+  vector_type = get_vector_type(space)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+end
+
+function ODEs.time_derivative(space::DistributedSingleFieldFESpace)
+  spaces = map(ODEs.time_derivative,local_views(space))
+  gids   = get_free_dof_ids(space)
+  trian  = get_triangulation(space)
+  vector_type = get_vector_type(space)
+  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+end
+
+for T in [:Real,:Nothing]
+  @eval begin
+    function Arrays.evaluate(space::DistributedTransientTrialFESpace, t::$T)
+      spaces = map(local_views(space)) do space
+        Arrays.evaluate(space,t)
+      end
+      gids  = get_free_dof_ids(space)
+      trian = get_triangulation(space)
+      vector_type = get_vector_type(space)
+      DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+    end
+  end
+end
+
+# SingleField CellField
+
+const DistributedTransientSingleFieldCellField = DistributedCellField{<:AbstractArray{<:ODEs.TransientSingleFieldCellField}}
+
+function ODEs.TransientCellField(f::DistributedCellField,derivatives::Tuple)
+  fields = map(local_views(f),map(local_views,derivatives)...) do f, derivatives...
+    ODEs.TransientCellField(f,Tuple(derivatives))
+  end
+  DistributedCellField(fields,get_triangulation(f),f.metadata)
+end
+
+function ODEs.time_derivative(f::DistributedTransientSingleFieldCellField)
+  fields = map(local_views(f)) do field
+    ODEs.time_derivative(field)
+  end
+  DistributedCellField(fields,get_triangulation(f))
+end
+
+# MultiField FESpace
+
+function ODEs.has_transient(space::DistributedMultiFieldFESpace)
+  getany(map(ODEs.has_transient,local_views(space)))
+end
+
+function ODEs.allocate_space(space::DistributedMultiFieldFESpace)
+  if !ODEs.has_transient(space)
+    return space
+  end
+  field_fe_space = map(ODEs.allocate_space,space.field_fe_space)
+  style = MultiFieldStyle(space)
+  spaces = to_parray_of_arrays(map(local_views,field_fe_space))
+  part_fe_spaces = map(s -> MultiFieldFESpace(s;style),spaces)
+  gids   = get_free_dof_ids(space)
+  vector_type = get_vector_type(space)
+  DistributedMultiFieldFESpace(field_fe_space,part_fe_spaces,gids,vector_type)
+end
+
+function ODEs.time_derivative(space::DistributedMultiFieldFESpace)
+  if !ODEs.has_transient(space)
+    return space
+  end
+  field_fe_space = map(ODEs.time_derivative,space.field_fe_space)
+  style = MultiFieldStyle(space)
+  spaces = to_parray_of_arrays(map(local_views,field_fe_space))
+  part_fe_spaces = map(s -> MultiFieldFESpace(s;style),spaces)
+  gids   = get_free_dof_ids(space)
+  vector_type = get_vector_type(space)
+  DistributedMultiFieldFESpace(field_fe_space,part_fe_spaces,gids,vector_type)
+end
+
+for T in [:Real,:Nothing]
+  @eval begin
+    function Arrays.evaluate(space::DistributedMultiFieldFESpace, t::$T)
+      if !ODEs.has_transient(space)
+        return space
+      end
+      field_fe_space = map(s->Arrays.evaluate(s,t),space.field_fe_space)
+      style = MultiFieldStyle(space)
+      spaces = to_parray_of_arrays(map(local_views,field_fe_space))
+      part_fe_spaces = map(s -> MultiFieldFESpace(s;style),spaces)
+      gids = get_free_dof_ids(space)
+      vector_type = get_vector_type(space)
+      DistributedMultiFieldFESpace(field_fe_space,part_fe_spaces,gids,vector_type)
+    end
+  end
+end
+
+# MultiField CellField
+
+const DistributedTransientMultiFieldCellField{A} = 
+  DistributedMultiFieldCellField{A,<:AbstractArray{<:ODEs.TransientMultiFieldCellField}}
+
+function ODEs.TransientCellField(f::DistributedMultiFieldCellField,derivatives::Tuple)
+  field_fe_fun = map(1:num_fields(f)) do i
+    f_i = f[i]
+    df_i = Tuple(map(df -> df[i],derivatives))
+    ODEs.TransientCellField(f_i,df_i)
+  end
+  fields = to_parray_of_arrays(map(local_views,field_fe_fun))
+  part_fe_fun = map(ODEs.TransientMultiFieldCellField,fields)
+  DistributedMultiFieldCellField(field_fe_fun,part_fe_fun,f.metadata)
+end
+
+function ODEs.TransientMultiFieldCellField(fields::AbstractVector{<:ODEs.TransientSingleFieldCellField})
+  cellfield = MultiFieldCellField(map(f -> f.cellfield,fields))
+  n_derivatives = length(first(fields).derivatives)
+  @check all(map(f -> length(f.derivatives) == n_derivatives,fields))
+  derivatives = Tuple(map(i -> MultiFieldCellField(map(f -> f.derivatives[i],fields)),1:n_derivatives))
+  TransientMultiFieldCellField(cellfield,derivatives,fields)
+end
+
+function ODEs.time_derivative(f::DistributedTransientMultiFieldCellField)
+  field_fe_fun = map(ODEs.time_derivative,f.field_fe_fun)
+  fields = to_parray_of_arrays(map(local_views,field_fe_fun))
+  part_fe_fun = map(ODEs.TransientMultiFieldCellField,fields)
+  DistributedMultiFieldCellField(field_fe_fun,part_fe_fun)
+end
diff --git a/src/TransientDistributedCellField.jl b/src/TransientDistributedCellField.jl
deleted file mode 100644
index 8934f40f..00000000
--- a/src/TransientDistributedCellField.jl
+++ /dev/null
@@ -1,107 +0,0 @@
-# Transient Distributed CellField
-abstract type TransientDistributedCellField <: DistributedCellDatum  end
-
-# Transient SingleField
-struct TransientSingleFieldDistributedCellField{A} <: TransientDistributedCellField
-  cellfield::A
-  derivatives::Tuple
-end
-
-local_views(f::TransientSingleFieldDistributedCellField) = local_views(f.cellfield)
-
-# Constructors
-function TransientFETools.TransientCellField(single_field::DistributedSingleFieldFEFunction,derivatives::Tuple)
-  TransientSingleFieldDistributedCellField(single_field,derivatives)
-end
-
-function TransientFETools.TransientCellField(single_field::DistributedCellField,derivatives::Tuple)
-TransientSingleFieldDistributedCellField(single_field,derivatives)
-end
-
-# Time derivative
-function ∂t(f::TransientDistributedCellField)
-  cellfield, derivatives = first_and_tail(f.derivatives)
-  TransientCellField(cellfield,derivatives)
-end
-
-∂tt(f::TransientDistributedCellField) = ∂t(∂t(f))
-
-# Integration related
-function Fields.integrate(f::TransientDistributedCellField,b::DistributedMeasure)
-  integrate(f.cellfield,b)
-end
-
-# Differential Operations
-Fields.gradient(f::TransientDistributedCellField) = gradient(f.cellfield)
-Fields.∇∇(f::TransientDistributedCellField) = ∇∇(f.cellfield)
-
-# Unary ops
-for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj)
-  @eval begin
-    ($op)(a::TransientDistributedCellField) = ($op)(a.cellfield)
-  end
-end
-
-# Binary ops
-for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/)
-  @eval begin
-    ($op)(a::TransientDistributedCellField,b::TransientDistributedCellField) = ($op)(a.cellfield,b.cellfield)
-    ($op)(a::TransientDistributedCellField,b::DistributedCellField) = ($op)(a.cellfield,b)
-    ($op)(a::DistributedCellField,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
-    ($op)(a::TransientDistributedCellField,b::Number) = ($op)(a.cellfield,b)
-    ($op)(a::Number,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
-    ($op)(a::TransientDistributedCellField,b::Function) = ($op)(a.cellfield,b)
-    ($op)(a::Function,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
-  end
-end
-
-Base.broadcasted(f,a::TransientDistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a.cellfield,b.cellfield)
-Base.broadcasted(f,a::TransientDistributedCellField,b::DistributedCellField) = broadcasted(f,a.cellfield,b)
-Base.broadcasted(f,a::DistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
-Base.broadcasted(f,a::Number,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
-Base.broadcasted(f,a::TransientDistributedCellField,b::Number) = broadcasted(f,a.cellfield,b)
-Base.broadcasted(f,a::Function,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
-Base.broadcasted(f,a::TransientDistributedCellField,b::Function) = broadcasted(f,a.cellfield,b)
-Base.broadcasted(a::typeof(*),b::typeof(∇),f::TransientDistributedCellField) = broadcasted(a,b,f.cellfield)
-Base.broadcasted(a::typeof(*),s::Fields.ShiftedNabla,f::TransientDistributedCellField) = broadcasted(a,s,f.cellfield)
-
-dot(::typeof(∇),f::TransientDistributedCellField) = dot(∇,f.cellfield)
-outer(::typeof(∇),f::TransientDistributedCellField) = outer(∇,f.cellfield)
-outer(f::TransientDistributedCellField,::typeof(∇)) = outer(f.cellfield,∇)
-cross(::typeof(∇),f::TransientDistributedCellField) = cross(∇,f.cellfield)
-
-Gridap.Arrays.evaluate!(cache,k::Operation,a::TransientDistributedCellField,b::DistributedCellField) = evaluate!(cache,k,a.cellfield,b)
-Gridap.Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::TransientDistributedCellField) = evaluate!(cache,k,a,b.cellfield)
-
-Base.:(∘)(f::Function,g::Tuple{TransientDistributedCellField,DistributedCellField}) = Operation(f)(g[1],g[2])
-Base.:(∘)(f::Function,g::Tuple{DistributedCellField,TransientDistributedCellField}) = Operation(f)(g[1],g[2])
-
-# Skeleton related
-function Base.getproperty(f::TransientDistributedCellField, sym::Symbol)
-  if sym in (:⁺,:plus,:⁻, :minus)
-    derivatives = ()
-    cellfield = DistributedCellField(f.cellfield,sym)
-    for iderivative in f.derivatives
-      derivatives = (derivatives...,DistributedCellField(iderivative))
-    end
-    return TransientSingleFieldCellField(cellfield,derivatives)
-  else
-    return getfield(f, sym)
-  end
-end
-
-Base.propertynames(x::TransientDistributedCellField, private::Bool=false) = propertynames(x.cellfield, private)
-
-for op in (:outer,:*,:dot)
-  @eval begin
-    ($op)(a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = ($op)(a.cellfield,b)
-    ($op)(a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = ($op)(a,b.cellfield)
-  end
-end
-
-Arrays.evaluate!(cache,k::Operation,a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = evaluate!(cache,k,a.cellfield,b)
-
-Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = evaluate!(cache,k,a,b.cellfield)
-
-CellData.jump(a::TransientDistributedCellField) = jump(a.cellfield)
-CellData.mean(a::TransientDistributedCellField) = mean(a.cellfield)
diff --git a/src/TransientFESpaces.jl b/src/TransientFESpaces.jl
deleted file mode 100644
index 938d86c8..00000000
--- a/src/TransientFESpaces.jl
+++ /dev/null
@@ -1,58 +0,0 @@
-# Functions for transient FE spaces
-
-Fields.evaluate(U::DistributedSingleFieldFESpace,t::Nothing) = U
-
-(U::DistributedSingleFieldFESpace)(t) = U
-
-∂t(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U)
-∂t(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂t.(U.field_fe_space))
-∂tt(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U)
-∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_fe_spaces))
-
-function TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace})
-  MultiFieldFESpace(spaces)
-end
-
-# Functions for transient FE Functions
-
-function ODETools.allocate_jacobian(
-  op::TransientFETools.TransientFEOperatorFromWeakForm,
-  t0::Real,
-  duh::Union{DistributedCellField,DistributedMultiFieldFEFunction},
-  cache)
-  _matdata_jacobians = TransientFETools.fill_initial_jacobians(op,t0,duh)
-  matdata = _vcat_distributed_matdata(_matdata_jacobians)
-  allocate_matrix(op.assem_t,matdata)
-end
-
-function ODETools.jacobians!(
-  A::AbstractMatrix,
-  op::TransientFETools.TransientFEOperatorFromWeakForm,
-  t::Real,
-  xh::TransientDistributedCellField,
-  γ::Tuple{Vararg{Real}},
-  cache)
-  _matdata_jacobians = TransientFETools.fill_jacobians(op,t,xh,γ)
-  matdata = _vcat_distributed_matdata(_matdata_jacobians)
-  assemble_matrix_add!(A,op.assem_t, matdata)
-  A
-end
-
-function _vcat_distributed_matdata(_matdata)
-  term_to_cellmat = map(a->a[1],local_views(_matdata[1]))
-  term_to_cellidsrows = map(a->a[2],local_views(_matdata[1]))
-  term_to_cellidscols = map(a->a[3],local_views(_matdata[1]))
-  for j in 2:length(_matdata)
-    term_to_cellmat_j = map(a->a[1],local_views(_matdata[j]))
-    term_to_cellidsrows_j = map(a->a[2],local_views(_matdata[j]))
-    term_to_cellidscols_j = map(a->a[3],local_views(_matdata[j]))
-    term_to_cellmat = map((a,b)->vcat(a,b),local_views(term_to_cellmat),local_views(term_to_cellmat_j))
-    term_to_cellidsrows = map((a,b)->vcat(a,b),local_views(term_to_cellidsrows),local_views(term_to_cellidsrows_j))
-    term_to_cellidscols = map((a,b)->vcat(a,b),local_views(term_to_cellidscols),local_views(term_to_cellidscols_j))
-  end
-  map( (a,b,c) -> (a,b,c),
-    local_views(term_to_cellmat),
-    local_views(term_to_cellidsrows),
-    local_views(term_to_cellidscols)
-  )
-end
diff --git a/src/TransientMultiFieldDistributedCellField.jl b/src/TransientMultiFieldDistributedCellField.jl
deleted file mode 100644
index 4bd280db..00000000
--- a/src/TransientMultiFieldDistributedCellField.jl
+++ /dev/null
@@ -1,63 +0,0 @@
-# Transient MultiField
-struct TransientMultiFieldDistributedCellField{A} <: TransientDistributedCellField
-  cellfield::A
-  derivatives::Tuple
-  transient_single_fields::Vector{<:TransientDistributedCellField} # used to iterate
-end
-
-local_views(f::TransientMultiFieldDistributedCellField) = local_views(f.cellfield)
-
-# Constructors
-function TransientFETools.TransientCellField(multi_field::DistributedMultiFieldFEFunction,derivatives::Tuple)
-  transient_single_fields = _to_transient_single_distributed_fields(multi_field,derivatives)
-  TransientMultiFieldDistributedCellField(multi_field,derivatives,transient_single_fields)
-end
-
-# Get single index
-function Base.getindex(f::TransientMultiFieldDistributedCellField,ifield::Integer)
-  single_field = f.cellfield[ifield]
-  single_derivatives = ()
-  for ifield_derivatives in f.derivatives
-    single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
-  end
-  TransientSingleFieldDistributedCellField(single_field,single_derivatives)
-end
-
-# Get multiple indices
-function Base.getindex(f::TransientMultiFieldDistributedCellField,indices::Vector{<:Int})
-  cellfield = DistributedMultiFieldCellField(f.cellfield[indices],DomainStyle(f.cellfield))
-  derivatives = ()
-  for derivative in f.derivatives
-    derivatives = (derivatives...,DistributedMultiFieldCellField(derivative[indices],DomainStyle(derivative)))
-  end
-  transient_single_fields = _to_transient_single_distributed_fields(cellfield,derivatives)
-  TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_fields)
-end
-
-function _to_transient_single_distributed_fields(multi_field,derivatives)
-  transient_single_fields = TransientDistributedCellField[]
-  for ifield in 1:num_fields(multi_field)
-    single_field = multi_field[ifield]
-    single_derivatives = ()
-    for ifield_derivatives in derivatives
-      single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
-    end
-    transient_single_field = TransientSingleFieldDistributedCellField(single_field,single_derivatives)
-    push!(transient_single_fields,transient_single_field)
-  end
-  transient_single_fields
-end
-
-# Iterate functions
-Base.iterate(f::TransientMultiFieldDistributedCellField)  = iterate(f.transient_single_fields)
-Base.iterate(f::TransientMultiFieldDistributedCellField,state)  = iterate(f.transient_single_fields,state)
-
-# Time derivative
-function ∂t(f::TransientMultiFieldDistributedCellField)
-  cellfield, derivatives = first_and_tail(f.derivatives)
-  transient_single_field_derivatives = TransientDistributedCellField[]
-  for transient_single_field in f.transient_single_fields
-    push!(transient_single_field_derivatives,∂t(transient_single_field))
-  end
-  TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_field_derivatives)
-end
diff --git a/test/HeatEquationTests.jl b/test/HeatEquationTests.jl
index 403d8ead..fa96dc11 100644
--- a/test/HeatEquationTests.jl
+++ b/test/HeatEquationTests.jl
@@ -1,7 +1,7 @@
 module HeatEquationTests
 
 using Gridap
-using Gridap.ODEs
+using Gridap.ODEs, Gridap.Algebra
 using GridapDistributed
 using PartitionedArrays
 using Test
@@ -10,65 +10,62 @@ function main(distribute,parts)
   ranks = distribute(LinearIndices((prod(parts),)))
   θ = 0.2
 
-  u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
-  u(t::Real) = x -> u(x,t)
-  f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
+  ut(t) = x -> (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
+  u = TimeSpaceFunction(ut)
+  ft(t) = x -> ∂t(u)(t,x) - Δ(u)(t,x)
+  f = TimeSpaceFunction(ft)
 
   domain = (0,1,0,1)
-  partition = (4,4)
-  model = CartesianDiscreteModel(ranks,parts,domain,partition)
+  model = CartesianDiscreteModel(ranks,parts,domain,(4,4))
 
   order = 2
-
   reffe = ReferenceFE(lagrangian,Float64,order)
-  V0 = FESpace(
-    model,
-    reffe,
-    conformity=:H1,
-    dirichlet_tags="boundary"
-  )
+  V0 = FESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
   U = TransientTrialFESpace(V0,u)
 
   Ω = Triangulation(model)
   degree = 2*order
   dΩ = Measure(Ω,degree)
 
-  #
-  m(u,v) = ∫(u*v)dΩ
-  a(u,v) = ∫(∇(v)⋅∇(u))dΩ
+  m(t,u,v) = ∫(u*v)dΩ
+  a(t,u,v) = ∫(∇(v)⋅∇(u))dΩ
   b(t,v) = ∫(v*f(t))dΩ
 
-  res(t,u,v) = a(u,v) + m(∂t(u),v) - b(t,v)
-  jac(t,u,du,v) = a(du,v)
-  jac_t(t,u,dut,v) = m(dut,v)
+  res(t,u,v) = a(t,u,v) + m(t,∂t(u),v) - b(t,v)
+  jac(t,u,du,v) = a(t,du,v)
+  jac_t(t,u,dut,v) = m(t,dut,v)
 
   op = TransientFEOperator(res,jac,jac_t,U,V0)
-  op_constant = TransientConstantMatrixFEOperator(m,a,b,U,V0)
+
+  assembler = SparseMatrixAssembler(U,V0,SubAssembledRows())
+  op_constant = TransientLinearFEOperator(
+    (a,m),b,U,V0,constant_forms=(true,true),assembler=assembler
+  )
 
   t0 = 0.0
   tF = 1.0
   dt = 0.1
 
   U0 = U(0.0)
-  uh0 = interpolate_everywhere(u(0.0),U0)
+  uh0 = interpolate_everywhere(u(0.0),U0);
 
   ls = LUSolver()
-  ode_solver = ThetaMethod(ls,dt,θ)
+  linear_ode_solver = ThetaMethod(ls,dt,θ)
+  sol_t_const = solve(linear_ode_solver,op_constant,t0,tF,uh0)
 
-  sol_t = solve(ode_solver,op,uh0,t0,tF)
-  sol_t_const = solve(ode_solver,op_constant,uh0,t0,tF)
+  nls = NewtonRaphsonSolver(ls,1.0e-6,10)
+  nonlinear_ode_solver = ThetaMethod(nls,dt,θ)
+  sol_t = solve(nonlinear_ode_solver,op,t0,tF,uh0)
 
   l2(w) = w*w
-
   tol = 1.0e-6
-
-  for (uh_tn, tn) in sol_t
+  for (tn, uh_tn) in sol_t
     e = u(tn) - uh_tn
     el2 = sqrt(sum( ∫(l2(e))dΩ ))
     @test el2 < tol
   end
 
-  for (uh_tn, tn) in sol_t_const
+  for (tn, uh_tn) in sol_t_const
     e = u(tn) - uh_tn
     el2 = sqrt(sum( ∫(l2(e))dΩ ))
     @test el2 < tol
diff --git a/test/MultiFieldTests.jl b/test/MultiFieldTests.jl
index cb53e322..9a0475ab 100644
--- a/test/MultiFieldTests.jl
+++ b/test/MultiFieldTests.jl
@@ -1,8 +1,7 @@
 module MultiFieldTests
 
 using Gridap
-using Gridap.FESpaces
-using Gridap.MultiField
+using Gridap.FESpaces, Gridap.MultiField, Gridap.Algebra
 using GridapDistributed
 using PartitionedArrays
 using Test
@@ -68,6 +67,15 @@ function main(distribute, parts, mfs)
   @test l2_error(p,ph1,dΩ) < 1.0e-9
   @test l2_error(u,uh2,dΩ) < 1.0e-9
   @test l2_error(p,ph2,dΩ) < 1.0e-9
+
+  a1(x,y) = ∫(x⋅y)dΩ
+  a2((u,p),(v,q)) = ∫(u⋅v + p⋅q)dΩ
+  A1 = assemble_matrix(a1,UxP,UxP)
+  A2 = assemble_matrix(a2,UxP,UxP)
+
+  x1 = allocate_in_domain(A1); fill!(x1,1.0)
+  x2 = allocate_in_domain(A2); fill!(x2,1.0)
+  @test norm(A1*x1-A2*x2) < 1.0e-9
 end
 
 function main(distribute, parts)
diff --git a/test/PLaplacianTests.jl b/test/PLaplacianTests.jl
index 63a2d447..2347dfc4 100644
--- a/test/PLaplacianTests.jl
+++ b/test/PLaplacianTests.jl
@@ -7,13 +7,6 @@ using PartitionedArrays
 using Test
 using SparseArrays
 
-# Overload required for the tests below
-function Base.copy(a::PSparseMatrix)
-  a_matrix_partition = similar(a.matrix_partition)
-  copy!(a_matrix_partition, a.matrix_partition)
-  PSparseMatrix(a_matrix_partition,a.row_partition,a.col_partition)
-end
-
 function main(distribute,parts)
   main(distribute,parts,FullyAssembledRows(),SparseMatrixCSR{0,Float64,Int})
   main(distribute,parts,SubAssembledRows(),SparseMatrixCSC{Float64,Int})
diff --git a/test/StokesOpenBoundaryTests.jl b/test/StokesOpenBoundaryTests.jl
index 8431e89b..3af86558 100644
--- a/test/StokesOpenBoundaryTests.jl
+++ b/test/StokesOpenBoundaryTests.jl
@@ -3,7 +3,7 @@ module StokesOpenBoundaryTests
 using Gridap
 using LinearAlgebra
 using Test
-using Gridap.ODEs
+using Gridap.ODEs, Gridap.Algebra
 using GridapDistributed
 using PartitionedArrays
 
@@ -12,16 +12,19 @@ function main(distribute,parts)
 
   θ = 0.5
 
-  u(x,t) = VectorValue(x[1],x[2])*t
-  u(t::Real) = x -> u(x,t)
+  ut(t) = x -> VectorValue(x[1],x[2])*t
+  u = TimeSpaceFunction(ut)
 
-  p(x,t) = (x[1]-x[2])*t
-  p(t::Real) = x -> p(x,t)
-  q(x) = t -> p(x,t)
+  pt(t) = x -> (x[1]-x[2])*t
+  p = TimeSpaceFunction(pt)
+  q(x) = t -> p(t,x)
 
-  f(t) = x -> ∂t(u)(t)(x)-Δ(u(t))(x)+ ∇(p(t))(x)
-  g(t) = x -> (∇⋅u(t))(x)
-  h(t) = x -> ∇(u(t))(x)⋅VectorValue(0.0,1.0) - p(t)(x)*VectorValue(0.0,1.0)
+  ft(t) = x -> ∂t(u)(t,x) - Δ(u)(t,x) + ∇(p)(t,x)
+  gt(t) = x -> (∇⋅u)(t,x)
+  ht(t) = x -> ∇(u)(t,x)⋅VectorValue(0.0,1.0) - p(t,x)*VectorValue(0.0,1.0)
+  f = TimeSpaceFunction(ft)
+  g = TimeSpaceFunction(gt)
+  h = TimeSpaceFunction(ht)
 
   domain = (0,1,0,1)
   partition = (4,4)
@@ -59,20 +62,16 @@ function main(distribute,parts)
   dΓ = Measure(Γ,degree)
 
   #
-  a(u,v) = ∫(∇(u)⊙∇(v))dΩ
-  b((v,q),t) = ∫(v⋅f(t))dΩ + ∫(q*g(t))dΩ + ∫(v⋅h(t))dΓ
-  m(ut,v) = ∫(ut⋅v)dΩ
+  a(t,u,v) = ∫(∇(u)⊙∇(v))dΩ
+  b(t,(v,q)) = ∫(v⋅f(t))dΩ + ∫(q*g(t))dΩ + ∫(v⋅h(t))dΓ
+  m(t,ut,v) = ∫(ut⋅v)dΩ
 
   X = TransientMultiFieldFESpace([U,P])
   Y = MultiFieldFESpace([V0,Q])
 
-  res(t,(u,p),(v,q)) = a(u,v) + m(∂t(u),v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t)
-  jac(t,(u,p),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ
-  jac_t(t,(u,p),(dut,dpt),(v,q)) = m(dut,v)
-
-  b((v,q)) = b((v,q),0.0)
-
-  mat((du1,du2),(v1,v2)) = a(du1,v1)+a(du2,v2)
+  res(t,(u,p),(v,q)) = a(t,u,v) + m(t,∂t(u),v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b(t,(v,q))
+  jac(t,(u,p),(du,dp),(v,q)) = a(t,du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ
+  jac_t(t,(u,p),(dut,dpt),(v,q)) = m(t,dut,v)
 
   U0 = U(0.0)
   P0 = P(0.0)
@@ -81,33 +80,28 @@ function main(distribute,parts)
   ph0 = interpolate_everywhere(p(0.0),P0)
   xh0 = interpolate_everywhere([uh0,ph0],X0)
 
-  op = TransientFEOperator(res,jac,jac_t,X,Y)
+  op = TransientFEOperator(res,(jac,jac_t),X,Y)
 
   t0 = 0.0
   tF = 1.0
   dt = 0.1
 
-  ls = LUSolver()
-  ode_solver = ThetaMethod(ls,dt,θ)
+  ls  = LUSolver()
+  nls = NewtonRaphsonSolver(ls,1.0e-6,10)
+  ode_solver = ThetaMethod(nls,dt,θ)
 
-  sol_t = solve(ode_solver,op,xh0,t0,tF)
+  sol_t = solve(ode_solver,op,t0,tF,xh0)
 
   l2(w) = w⋅w
-
-
   tol = 1.0e-6
-  _t_n = t0
-
-  result = Base.iterate(sol_t)
-
-  for (xh_tn, tn) in sol_t
+  for (tn, xh_tn) in sol_t
     uh_tn = xh_tn[1]
     ph_tn = xh_tn[2]
-    writevtk(Ω,"output/tmp_stokes_OB_sol_$tn.vtu",cellfields=["u"=>uh_tn,"p"=>ph_tn])
+    #writevtk(Ω,"output/tmp_stokes_OB_sol_$tn.vtu",cellfields=["u"=>uh_tn,"p"=>ph_tn])
     e = u(tn) - uh_tn
-    el2 = sqrt(sum( ∫(l2(e))dΩ ))
+    el2 = sqrt(sum(∫(l2(e))dΩ))
     e = p(tn) - ph_tn
-    el2 = sqrt(sum( ∫(l2(e))dΩ ))
+    el2 = sqrt(sum(∫(l2(e))dΩ))
     @test el2 < tol
   end
 end
diff --git a/test/TransientDistributedCellFieldTests.jl b/test/TransientDistributedCellFieldTests.jl
index e356abe4..b6328677 100644
--- a/test/TransientDistributedCellFieldTests.jl
+++ b/test/TransientDistributedCellFieldTests.jl
@@ -2,8 +2,8 @@ module TransientDistributedCellFieldTests
 
 using Gridap
 using GridapDistributed
-using Gridap.ODEs.ODETools: ∂t, ∂tt
-using Gridap.ODEs.TransientFETools: TransientCellField
+using Gridap.ODEs: ∂t, ∂tt
+using Gridap.ODEs: TransientCellField
 using PartitionedArrays
 using Test
 
@@ -27,16 +27,13 @@ function main(distribute,parts)
   @test isa(dda(0),GridapDistributed.DistributedCellField)
 
   b(t) = TransientCellField(a(t),(da(t),dda(t)))
-  @test isa(b(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(b(0),GridapDistributed.TransientSingleFieldDistributedCellField)
+  @test isa(b(0),GridapDistributed.DistributedTransientSingleFieldCellField)
 
   db(t) = ∂t(b(t))
-  @test isa(db(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(db(0),GridapDistributed.TransientSingleFieldDistributedCellField)
+  @test isa(db(0),GridapDistributed.DistributedTransientSingleFieldCellField)
 
   ddb(t) = ∂t(db(t))
-  @test isa(ddb(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(ddb(0),GridapDistributed.TransientSingleFieldDistributedCellField)
+  @test isa(ddb(0),GridapDistributed.DistributedTransientSingleFieldCellField)
 
   @test (∑(∫(a(0.5))dΩ)) ≈ 0.25
   @test (∑(∫(da(0.5))dΩ)) ≈ 1.0
diff --git a/test/TransientMultiFieldDistributedCellFieldTests.jl b/test/TransientMultiFieldDistributedCellFieldTests.jl
index bdb6eca6..ca83333e 100644
--- a/test/TransientMultiFieldDistributedCellFieldTests.jl
+++ b/test/TransientMultiFieldDistributedCellFieldTests.jl
@@ -2,9 +2,9 @@ module TransientMultiFieldDistributedCellFieldTests
 
 using Gridap
 using GridapDistributed
-using Gridap.ODEs.ODETools: ∂t, ∂tt
-using Gridap.ODEs.TransientFETools: TransientCellField
-using Gridap.ODEs.TransientFETools: TransientTrialFESpace, TransientMultiFieldFESpace
+using Gridap.ODEs: ∂t, ∂tt
+using Gridap.ODEs: TransientCellField
+using Gridap.ODEs: TransientTrialFESpace, TransientMultiFieldFESpace
 using PartitionedArrays
 using Test
 
@@ -36,28 +36,22 @@ function main(distribute,parts)
   @test isa(dda(0),GridapDistributed.DistributedMultiFieldFEFunction)
 
   b(t) = TransientCellField(a(t),(da(t),dda(t)))
-  @test isa(b(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(b(0),GridapDistributed.TransientMultiFieldDistributedCellField)
+  @test isa(b(0),GridapDistributed.DistributedTransientMultiFieldCellField)
 
   db(t) = ∂t(b(t))
-  @test isa(db(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(db(0),GridapDistributed.TransientMultiFieldDistributedCellField)
+  @test isa(db(0),GridapDistributed.DistributedTransientMultiFieldCellField)
 
   ddb(t) = ∂t(db(t))
-  @test isa(ddb(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(ddb(0),GridapDistributed.TransientMultiFieldDistributedCellField)
+  @test isa(ddb(0),GridapDistributed.DistributedTransientMultiFieldCellField)
 
   b1(t) = b(t)[1]
-  @test isa(b1(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(b1(0),GridapDistributed.TransientSingleFieldDistributedCellField)
+  @test isa(b1(0),GridapDistributed.DistributedTransientSingleFieldCellField)
 
   db1(t) = ∂t(b1(t))
-  @test isa(db1(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(db1(0),GridapDistributed.TransientSingleFieldDistributedCellField)
+  @test isa(db1(0),GridapDistributed.DistributedTransientSingleFieldCellField)
 
   ddb1(t) = ∂t(db1(t))
-  @test isa(ddb1(0),GridapDistributed.TransientDistributedCellField)
-  @test isa(ddb1(0),GridapDistributed.TransientSingleFieldDistributedCellField)
+  @test isa(ddb1(0),GridapDistributed.DistributedTransientSingleFieldCellField)
 
   @test (∑(∫(b(0.5)[1])dΩ)) == (∑(∫(b1(0.5))dΩ))
   @test (∑(∫(db(0.5)[1])dΩ)) == (∑(∫(db1(0.5))dΩ))