diff --git a/NEWS.md b/NEWS.md
index 24ebcf1..65ba24d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,12 @@ 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).
 
+## Unreleased
+
+### Fixed
+
+- Fixed distributed interpolators for Vector-Valued FESpaces. Since PR[#152](https://github.com/gridap/GridapDistributed.jl/pull/152).
+
 ## [0.4.3] 2024-07-18
 
 ### Added
diff --git a/src/CellData.jl b/src/CellData.jl
index 95518ce..6befc88 100644
--- a/src/CellData.jl
+++ b/src/CellData.jl
@@ -355,73 +355,99 @@ f(s) instead of s(f)?
 end
 
 # Interpolation at arbitrary points (returns -Inf if the point is not found)
-Arrays.evaluate!(cache,f::DistributedCellField,x::Point) = evaluate!(cache,Interpolable(f),x)
-Arrays.evaluate!(cache,f::DistributedCellField,x::AbstractVector{<:Point}) = evaluate!(cache,Interpolable(f),x)
+Arrays.evaluate!(cache,f::DistributedCellField,x::Point) = evaluate(Interpolable(f),x)
+Arrays.evaluate!(cache,f::DistributedCellField,x::AbstractVector{<:Point}) = evaluate(Interpolable(f),x)
 
-struct DistributedInterpolable{A} <: Function
+struct DistributedInterpolable{Tx,Ty,A} <: Function
   interps::A
+  function DistributedInterpolable(interps::AbstractArray{<:Interpolable})
+    Tx,Ty = map(interps) do I
+      fi = I.uh
+      trian = get_triangulation(fi)
+      x = mean(testitem(get_cell_coordinates(trian)))
+      return typeof(x), return_type(fi,x)
+    end |> tuple_of_arrays
+    Tx = getany(Tx)
+    Ty = getany(Ty)
+    A  = typeof(interps)
+    new{Tx,Ty,A}(interps)
+  end
 end
+
 local_views(a::DistributedInterpolable) = a.interps
 
 function Interpolable(f::DistributedCellField;kwargs...) 
-  interps = map(local_views(f)) do fun
-      Interpolable(fun,kwargs...)
+  interps = map(local_views(f)) do f
+    Interpolable(f,kwargs...)
   end
   DistributedInterpolable(interps)
 end
 
 (a::DistributedInterpolable)(x) = evaluate(a,x)
 
-function Gridap.Arrays.return_cache(I::DistributedInterpolable,x::Point)
-  caches = map(local_views(I)) do fi
-      trian = get_triangulation(fi.uh)
-      y=mean(testitem(get_cell_coordinates(trian)))
-      @check typeof(testitem(x)) == typeof(y) "Can only evaluate DistributedInterpolable at physical points of the same dimension of the underlying triangulation"
-      return_cache(fi,y)
+Arrays.return_cache(f::DistributedInterpolable,x::Point) = return_cache(f,[x])
+Arrays.evaluate!(caches,I::DistributedInterpolable,x::Point) = first(evaluate!(caches,I,[x]))
+
+function Arrays.return_cache(I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty}
+  msg = "Can only evaluate DistributedInterpolable at physical points of the same dimension of the underlying triangulation"
+  @check Tx == eltype(x) msg
+  caches = map(local_views(I)) do I
+    trian = get_triangulation(I.uh)
+    y = mean(testitem(get_cell_coordinates(trian)))
+    return_cache(I,y)
   end
-  caches 
+  caches
 end
-Gridap.Arrays.return_cache(f::DistributedInterpolable,x::AbstractVector{<:Point}) = Gridap.Arrays.return_cache(f,testitem(x))
 
-function Gridap.Arrays.evaluate!(caches,I::DistributedInterpolable,x::Point)
-  y=map(local_views(I),local_views(caches)) do fi,cache
+function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty}
+  _allgather(x) = PartitionedArrays.getdata(getany(gather(x;destination=:all)))
+
+  # Evaluate in local portions of the domain. Only keep points inside the domain.
+  nx = length(x)
+  my_ids, my_vals = map(local_views(I),local_views(cache)) do I, cache
+    ids  = Vector{Int}(undef,nx)
+    vals = Vector{Ty}(undef,nx)
+    k = 1
+    yi = zero(Ty)
+    for (i,xi) in enumerate(x)
+      inside = true
       try
-          evaluate!(cache,fi,x)
+        yi = evaluate!(cache,I,xi)
       catch
-          -Inf
+        inside = false
       end
-  end
-  # reduce(max,y)
-  z=gather(y)
-  map_main(local_views(z)) do zi
-      reduce(max,zi)
-    end
-end
-
-function Gridap.Arrays.evaluate!(caches,I::DistributedInterpolable,v::AbstractVector{<:Point})
-  n=length(local_views(I))
-  m=length(v)
-  y=map(local_views(I),local_views(caches)) do fi,cache
-      w=Vector{Float64}(undef,m)
-      for (i,x) in enumerate(v)
-          try
-              w[i]=evaluate!(cache,fi,x)
-          catch
-              w[i]=-Inf
-          end
+      if inside
+        ids[k] = i
+        vals[k] = yi
+        k += 1
       end
-      return w
     end
-  # z=gather(y,destination=:all)
-  z=gather(y)
-  map_main(local_views(z)) do zi
-      w=Vector{Float64}(undef,m)
-      for i=0:m-1
-        w[i+1]=reduce(max,zi.data[zi.ptrs[1:n].+i])
-      end
-      return w
+    resize!(ids,k-1)
+    resize!(vals,k-1)
+    return ids, vals
+  end |> tuple_of_arrays
+
+  # Communicate results, so that every (id,value) pair is known by every process
+  if Ty <: VectorValue
+    D = num_components(Ty)
+    vals_d = Vector{Vector{eltype(Ty)}}(undef,D)
+    for d in 1:D
+      my_vals_d = map(y_p -> map(y_p_i -> y_p_i[d],y_p),my_vals)
+      vals_d[d] = _allgather(my_vals_d)
     end
-  # reduce((v,w)->broadcast(max,v,w),y)
+    vals = map(VectorValue,vals_d...)
+  else
+    vals = _allgather(my_vals)
+  end
+  ids = _allgather(my_ids)
+
+  # Combine results
+  w = Vector{Ty}(undef,nx)
+  for (i,v) in zip(ids,vals)
+    w[i] = v
+  end
+
+  return w
 end
 
 # Support for distributed Dirac deltas
diff --git a/src/Geometry.jl b/src/Geometry.jl
index 3b9b4dc..baf1967 100644
--- a/src/Geometry.jl
+++ b/src/Geometry.jl
@@ -237,13 +237,13 @@ function emit_cartesian_descriptor(
   return DistributedCartesianDescriptor(new_ranks,new_mesh_partition,new_desc)
 end
 
-const DistributedCartesianDiscreteModel{Dc,Dp,A,B,C} = 
+const DistributedCartesianDiscreteModel{Dc,Dp,A,B,C} =
   GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:CartesianDiscreteModel},B,<:DistributedCartesianDescriptor}
 
 function Geometry.CartesianDiscreteModel(
   ranks::AbstractArray{<:Integer}, # Distributed array with the rank IDs
   parts::NTuple{N,<:Integer},      # Number of ranks (parts) in each direction
-  args...;isperiodic=map(i->false,parts),kwargs...) where N 
+  args...;isperiodic=map(i->false,parts),kwargs...) where N
 
   desc = CartesianDescriptor(args...;isperiodic=isperiodic,kwargs...)
   nc = desc.partition
@@ -262,7 +262,7 @@ function Geometry.CartesianDiscreteModel(
     models = map(ranks,upartition) do rank, upartition
       cmin = gcids[first(upartition)]
       cmax = gcids[last(upartition)]
-      CartesianDiscreteModel(desc,cmin,cmax)  
+      CartesianDiscreteModel(desc,cmin,cmax)
     end
     gids = PRange(upartition)
     metadata = DistributedCartesianDescriptor(ranks,parts,desc)
@@ -271,10 +271,10 @@ function Geometry.CartesianDiscreteModel(
 end
 
 function _cartesian_model_with_periodic_bcs(ranks,parts,desc)
-  # We create and extended CartesianDescriptor for the local models: 
-  # If a direction is periodic and partitioned: 
+  # We create and extended CartesianDescriptor for the local models:
+  # If a direction is periodic and partitioned:
   #   - we add a ghost cell at either side, which will be made periodic by the index partition.
-  #   - We move the origin to accomodate the new cells. 
+  #   - We move the origin to accomodate the new cells.
   #   - We turn OFF the periodicity in the local model, since periodicity will be taken care of
   #     by the global index partition.
   _map = desc.map #! Important: the map should be periodic if you want to integrate on the ghost cells.
@@ -288,7 +288,7 @@ function _cartesian_model_with_periodic_bcs(ranks,parts,desc)
   end |> tuple_of_arrays
   _desc = CartesianDescriptor(Point(_origin),_sizes,_partition;map=_map,isperiodic=_isperiodic)
 
-  # We create the global index partition, which has the original number of cells per direction. 
+  # We create the global index partition, which has the original number of cells per direction.
   # Globally, the periodicity is turned ON in the directions which are periodic and partitioned
   # (if a direction is not partitioned, the periodicity is handled locally).
   ghost = map(i->true,parts)
@@ -399,7 +399,7 @@ function Geometry.DiscreteModel(
           cell_to_mask[jcell] = true
         end
       end
-    end 
+    end
     lcell_to_cell = findall(cell_to_mask)
     lcell_to_part = zeros(Int32,length(lcell_to_cell))
     lcell_to_part .= cell_to_part[lcell_to_cell]
@@ -408,11 +408,11 @@ function Geometry.DiscreteModel(
 
   partition = map(parts,lcell_to_cell,lcell_to_part) do part, lcell_to_cell, lcell_to_part
     LocalIndices(ncells, part, lcell_to_cell, lcell_to_part)
-  end 
+  end
 
-  # This is required to provide the hint that the communication 
-  # pattern underlying partition is symmetric, so that we do not have 
-  # to execute the algorithm the reconstructs the reciprocal in the 
+  # This is required to provide the hint that the communication
+  # pattern underlying partition is symmetric, so that we do not have
+  # to execute the algorithm the reconstructs the reciprocal in the
   # communication graph
   assembly_neighbors(partition;symmetric=true)
 
@@ -427,7 +427,7 @@ end
 
 # UnstructuredDiscreteModel
 
-const DistributedUnstructuredDiscreteModel{Dc,Dp,A,B,C} = 
+const DistributedUnstructuredDiscreteModel{Dc,Dp,A,B,C} =
   GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:UnstructuredDiscreteModel},B,C}
 
 function Geometry.UnstructuredDiscreteModel(model::GenericDistributedDiscreteModel)
@@ -439,9 +439,9 @@ end
 
 # Simplexify
 
-function Geometry.simplexify(model::DistributedDiscreteModel)
+function Geometry.simplexify(model::DistributedDiscreteModel;kwargs...)
   _model = UnstructuredDiscreteModel(model)
-  ref_model = refine(_model, refinement_method = "simplexify")
+  ref_model = refine(_model; refinement_method = "simplexify", kwargs...)
   return UnstructuredDiscreteModel(Adaptivity.get_model(ref_model))
 end
 
@@ -696,7 +696,7 @@ function add_ghost_cells(dmodel::DistributedDiscreteModel{Dm},
 
     cache=fetch_vector_ghost_values_cache(mcell_intrian,partition(gids))
     fetch_vector_ghost_values!(mcell_intrian,cache) |> wait
-    
+
     dreffes=map(local_views(dmodel)) do model
       ReferenceFE{Dt}
     end
@@ -721,7 +721,7 @@ function generate_cell_gids(dmodel::DistributedDiscreteModel{Dm},
     # count number owned cells
     notcells, tcell_to_mcell = map(
       local_views(dmodel),local_views(dtrian),PArrays.partition(mgids)) do model,trian,partition
-      lid_to_owner = local_to_owner(partition)  
+      lid_to_owner = local_to_owner(partition)
       part = part_id(partition)
       glue = get_glue(trian,Val(Dt))
       @assert isa(glue,FaceToFaceGlue)
@@ -756,7 +756,7 @@ function generate_cell_gids(dmodel::DistributedDiscreteModel{Dm},
 
     # Prepare new partition
     ngtcells = reduction(+,notcells,destination=:all,init=zero(eltype(notcells)))
-    partition = map(ngtcells, 
+    partition = map(ngtcells,
                     mcell_to_gtcell,
                     tcell_to_mcell,
                     PArrays.partition(mgids)) do ngtcells,mcell_to_gtcell,tcell_to_mcell,partition
diff --git a/test/AdaptivityUnstructuredTests.jl b/test/AdaptivityUnstructuredTests.jl
index fdf1939..588f15a 100644
--- a/test/AdaptivityUnstructuredTests.jl
+++ b/test/AdaptivityUnstructuredTests.jl
@@ -74,7 +74,7 @@ function main(distribute,parts,ncells)
   child2 = refine(parent1, refinement_method = "simplexify" )
   test_adaptivity(ranks,parent1,child2)
 
-  parent2 = simplexify(parent1)
+  parent2 = simplexify(parent1,positive=true)
 
   if Dc == 2
     i_am_main(ranks) && println("UnstructuredAdaptivityTests: nvb")
@@ -92,4 +92,4 @@ function main(distribute)
   main(distribute,(2,2,1),(4,4,4))
 end
 
-end # module AdaptivityUnstructuredTests
\ No newline at end of file
+end # module AdaptivityUnstructuredTests
diff --git a/test/CellDataTests.jl b/test/CellDataTests.jl
index 2485374..3baa2b8 100644
--- a/test/CellDataTests.jl
+++ b/test/CellDataTests.jl
@@ -58,7 +58,7 @@ function main(distribute,parts)
   u3 = CellField(2.0,Ω)
   u = _my_op∘(u1,u2,u3)
 
- order = 1
+  order = 1
   reffe = ReferenceFE(lagrangian,Float64,order)
   V = TestFESpace(model,reffe)
   uh = interpolate_everywhere(x->x[1]+x[2],V)
@@ -67,15 +67,21 @@ function main(distribute,parts)
   x3 = Point(0.9,0.9)
   v = [x1,x2,x3]
 
-  u1 = uh(x1)
-  u2 = uh(x2)
-  uv = uh(v)
+  @test uh(x1) == 0.2
+  @test uh(x2) == 1.0
+  @test uh(v) == [0.2,1.0,1.8]
 
-  map_main(u1,u2,uv) do u1,u2,v
-    @test u1 == 0.2
-    @test u2 == 1.0
-    @test v ==[0.2,1.0,1.8]
-  end
+  reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+  V = TestFESpace(model,reffe)
+  uh = interpolate_everywhere(x->x,V)
+  x1 = Point(0.1,0.1)
+  x2 = Point(0.1,0.9)
+  x3 = Point(0.9,0.9)
+  v = [x1,x2,x3]
+
+  @test uh(x1) == x1
+  @test uh(x2) == x2
+  @test uh(v)  == v
 
   # Point δ
   δ=DiracDelta{0}(model;tags=2)
@@ -89,7 +95,7 @@ function main(distribute,parts)
   @test sum(δ(f)) ≈ 8.0
   @test sum(δ(3.0)) ≈ 12.0
   @test sum(δ(x->2*x)) ≈ VectorValue(16.0,0.0)
- 
+
 end
 
 end # module