From d5bc4c962dd7874949c5be6a968206a3e987446c Mon Sep 17 00:00:00 2001
From: principe <principe@cimne.upc.edu>
Date: Thu, 16 May 2024 13:22:35 +0200
Subject: [PATCH] Implemented DistributedInterpolable to evaluate
 DistributedCellFields at arbitrary points

---
 src/CellData.jl          | 59 +++++++++++++++++++++++++++++++++++++++-
 src/GridapDistributed.jl |  1 +
 test/CellDataTests.jl    | 13 +++++++++
 3 files changed, 72 insertions(+), 1 deletion(-)

diff --git a/src/CellData.jl b/src/CellData.jl
index 0814b398..033f7011 100644
--- a/src/CellData.jl
+++ b/src/CellData.jl
@@ -352,4 +352,61 @@ function Gridap.Arrays.evaluate!(cache, ::DistributedCellField, ::DistributedCel
   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
+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)
+
+struct DistributedInterpolable{A} <: Function
+  interps::A
+end
+local_views(a::DistributedInterpolable) = a.interps
+
+function Interpolable(f::DistributedCellField;kwargs...) 
+  interps = map(local_views(f)) do fun
+      Interpolable(fun,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)
+  end
+  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
+      try
+          evaluate!(cache,fi,x)
+      catch
+          -Inf
+      end
+  end
+  reduce(max,y)
+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
+      end
+      return w
+  end
+  reduce((v,w)->broadcast(max,v,w),y)
+end
diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl
index 0472ef61..8c20a975 100644
--- a/src/GridapDistributed.jl
+++ b/src/GridapDistributed.jl
@@ -29,6 +29,7 @@ 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.CellData: Interpolable
 
 export FullyAssembledRows
 export SubAssembledRows
diff --git a/test/CellDataTests.jl b/test/CellDataTests.jl
index 3bc15911..3e7127a0 100644
--- a/test/CellDataTests.jl
+++ b/test/CellDataTests.jl
@@ -58,6 +58,19 @@ function main(distribute,parts)
   u3 = CellField(2.0,Ω)
   u = _my_op∘(u1,u2,u3)
 
+  order = 1
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe)
+  uh = interpolate_everywhere(x->x[1]+x[2],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) ==0.2
+  @test uh(x2) ==1.0
+  @test uh(v) ==[0.2,1.0,1.8]
+ 
 end
 
 end # module