Allow evaluating FE functions at arbitrary points

Project.toml

@@ -18,6 +18,7 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -40,6 +41,7 @@ JLD2 = "0.1.11, 0.3, 0.4"
 JSON = "0.21.0"
 LineSearches = "7.0.1"
 NLsolve = "4.3.0"
+NearestNeighbors = "0.4.8"
 QuadGK = "2.3.1, 2.4"
 SparseMatricesCSR = "0.6"
 StaticArrays = "0.12.1, 1.0"

src/CellData/CellData.jl

@@ -10,8 +10,8 @@ using DocStringExtensions
 using FillArrays
 
 using Gridap.Helpers
-using Gridap.Arrays
 using Gridap.Algebra
+using Gridap.Arrays
 using Gridap.TensorValues
 using Gridap.Fields
 using Gridap.ReferenceFEs
@@ -20,6 +20,7 @@ using Gridap.Geometry
 import Gridap.Arrays: lazy_append
 import Gridap.Arrays: get_array
 import Gridap.Arrays: evaluate!
+import Gridap.Arrays: return_cache
 import Gridap.Fields: gradient
 import Gridap.Fields: ∇∇
 import Gridap.Fields: integrate

src/CellData/CellFields.jl

@@ -212,16 +212,6 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS()
 
 (a::CellField)(x) = evaluate(a,x)
 
-function evaluate!(cache,f::CellField,x::Point)
-  @notimplemented """\n
-  Evaluation of a CellField at a given Point is not implemented yet.
-
-  This is a feature that we want to have at some point in Gridap.
-  If you are ready to help with this implementation, please contact the
-  Gridap administrators.
-  """
-end
-
 function evaluate!(cache,f::CellField,x::CellPoint)
   _f, _x = _to_common_domain(f,x)
   cell_field = get_data(_f)
@@ -334,11 +324,11 @@ struct OperationCellField{DS} <: CellField
          r = Fields.BroadcastingFieldOpMap(op.op)(axi...)
       catch
         @unreachable """\n
-        It is not possible to perform operation $(op.op) on the given cell fields.
+        It is not possible to perform the operation "$(op.op)" on the given cell fields.
 
-        See the caught error for more information. (If you are using Visual
+        See the caught error for more information. (If you are using the Visual
           Studio Code REPL you might not see the caught error, please use the
-          command-line REPL).
+          command-line REPL instead).
         """
       end
     end

src/MultiField/MultiField.jl

@@ -15,12 +15,21 @@ using Gridap.TensorValues
 using Gridap.CellData
 using Gridap.Fields
 
+using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values
 using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis
+using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, ExtrusionPolytope,
+                           get_faces, get_node_coordinates, get_polytope,
+                           num_cell_dims
 
 using FillArrays
 using SparseArrays
 using LinearAlgebra
-import BlockArrays
+using BlockArrays
+using NearestNeighbors
+using StaticArrays
+
+import Gridap.Arrays: evaluate!
+import Gridap.Arrays: return_cache
 
 export num_fields
 export compute_field_offsets

src/MultiField/MultiFieldCellFields.jl

@@ -39,3 +39,167 @@ num_fields(a::MultiFieldCellField) = length(a.single_fields)
 Base.getindex(a::MultiFieldCellField,i::Integer) = a.single_fields[i]
 Base.iterate(a::MultiFieldCellField)  = iterate(a.single_fields)
 Base.iterate(a::MultiFieldCellField,state)  = iterate(a.single_fields,state)
+
+# Evaluation of CellFields
+
+# This code lives here (instead of in the CellData module) because we
+# need access to MultiField functionality.
+
+"""
+   dist = distance(polytope::ExtrusionPolytope,
+                   inv_cmap::Field,
+                   x::Point)
+
+Calculate distance from point `x` to the polytope. The polytope is
+given by its type and by the inverse cell map, i.e. by the map from
+the physical to the reference space.
+
+Positive distances are outside the polytope, negative distances are
+inside the polytope.
+
+The distance is measured in an unspecified norm, currently the L∞
+norm.
+"""
+function distance(polytope::ExtrusionPolytope, inv_cmap::Field, x::Point)
+  extrusion = polytope.extrusion
+  isempty(extrusion) && return zero(eltype(x))
+  p = inv_cmap(x)
+  if all(e == HEX_AXIS for e in extrusion)
+    # Boundaries are at `a=0` and `a=1` in each direction
+    return maximum(max(0 - a, a - 1) for a in p)
+  elseif all(e == TET_AXIS for e in extrusion)
+    # Calculate barycentric coordinates
+    λ = Point(p..., 1 - sum(p))
+    return maximum(-λ)
+  else
+    @notimplemented "Only hypercubes and simplices are implemented so far"
+  end
+end
+
+function return_cache(f::CellField,x::Point)
+  trian = get_triangulation(f)
+  topo = GridTopology(trian)    # Note: this is expensive
+  vertex_coordinates = Geometry.get_vertex_coordinates(topo)
+  kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates))
+  D = num_cell_dims(trian)
+  vertex_to_cells = get_faces(topo,0,D)
+  cell_to_ctype = get_cell_type(trian)
+  ctype_to_reffe = get_reffes(trian)
+  ctype_to_polytope = map(get_polytope, ctype_to_reffe)
+  cell_map = get_cell_map(trian)
+  cache1 = kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map
+
+  cell_f = get_array(f)
+  cell_f_cache = array_cache(cell_f)
+  cf = testitem(cell_f)
+  f_cache = return_cache(cf,x)
+  cache2 = cell_f_cache, f_cache, cell_f, f
+
+  return cache1,cache2
+end
+
+function point_to_cell!(cache,f::CellField,x::Point)
+  kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache
+
+  # Find nearest vertex
+  id,dist = nn(kdtree, SVector(Tuple(x)))
+
+  # Find all neighbouring cells
+  cells = vertex_to_cells[id]
+  @assert !isempty(cells)
+
+  # Calculate the distance from the point to all the cells. Without
+  # round-off, and with non-overlapping cells, the distance would be
+  # negative for exactly one cell and positive for all other ones. Due
+  # to round-off, the distance can be slightly negative or slightly
+  # positive for points near cell boundaries, in particular near
+  # vertices. In this case, choose the cell with the smallest
+  # distance, and check that the distance (if positive) is at most at
+  # round-off level.
+  T = eltype(dist)
+  function cell_distance(cell::Integer)
+    ctype = cell_to_ctype[cell]
+    polytope = ctype_to_polytope[ctype]
+    cmap = cell_map[cell]
+    inv_cmap = inverse_map(cmap)
+    return distance(polytope, inv_cmap, x)
+  end
+  # findmin, without allocating an array
+  cell = zero(eltype(cells))
+  dist = T(Inf)
+  for jcell in cells
+    jdist = cell_distance(jcell)
+    if jdist < dist
+      cell = jcell
+      dist = jdist
+    end
+  end
+  # Ensure the point is inside one of the cells, up to round-off errors
+  @check dist ≤ 1000eps(T) "Point is not inside any cell"
+
+  return cell
+end
+
+function evaluate!(cache,f::CellField,x::Point)
+  cache1,cache2 = cache
+  cell_f_cache, f_cache, cell_f, f₀ = cache2
+  @check f === f₀ "Wrong cache"
+
+  cell = point_to_cell!(cache1,f,x)
+  cf = getindex!(cell_f_cache, cell_f, cell)
+  fx = evaluate!(f_cache, cf, x)
+  return fx
+end
+
+return_cache(f::CellField,xs::AbstractVector{<:Point}) = return_cache(f,testitem(xs))
+
+# # Simple version:
+# function evaluate!(cache,f::CellField,xs::AbstractVector{<:Point})
+#   return map(x->evaluate!(cache,f,x), xs)
+# end
+
+function make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int)
+  ni = length(i2j)
+  @assert nj≥0
+
+  p = sortperm(i2j)
+  # i2j[p] is a sorted array of js
+
+  data = p
+  ptrs = Array{Int}(undef, nj+1)
+  i2lis = Array{Int}(undef, ni)
+  jprev = zero(eltype(i2j))
+  liprev = 0
+  for (n,i) in enumerate(p)
+    j = i2j[i]
+    @assert jprev≤j≤nj
+    li = (j==jprev ? liprev : 0) + 1
+    ptrs[jprev+1:j] .= n
+    i2lis[i] = li
+    jprev = j
+    liprev = li
+  end
+  ptrs[jprev+1:nj+1] .= ni+1
+  j2is = Table(data,ptrs)  
+
+  return j2is,i2lis
+end
+
+# Efficient version:
+function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point})
+  cache1,cache2 = return_cache(f,testitem(point_to_x))
+  kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache1
+  cell_f_cache, f_cache, cell_f, f₀ = cache2
+  @check f === f₀ "Wrong cache"
+
+  ncells = length(cell_map)
+  x_to_cell(x) = point_to_cell!(cache1,f,x)
+  point_to_cell = map(x_to_cell,point_to_x)
+  cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells)
+  cell_to_xs = lazy_map(Broadcasting(Reindex(point_to_x)),cell_to_points)
+  cell_to_f = get_array(f)
+  cell_to_fxs = lazy_map(evaluate,cell_to_f,cell_to_xs)
+  point_to_fxs = lazy_map(Reindex(cell_to_fxs),point_to_cell)
+  point_to_fx = lazy_map(getindex,point_to_fxs,point_to_lpoint)
+  collect(point_to_fx)          # Collect into a plain array
+end

test/CellDataTests/CellFieldsTests.jl

@@ -4,10 +4,12 @@ using Test
 using BlockArrays
 using FillArrays
 using Gridap.Arrays
-using Gridap.Fields
 using Gridap.TensorValues
+using Gridap.Fields
+using Gridap.ReferenceFEs
 using Gridap.Geometry
 using Gridap.CellData
+using Gridap.FESpaces
 
 domain = (0,1,0,1)
 cells = (3,3)
@@ -152,9 +154,6 @@ test_array(hx_N,collect(hx_N))
 
 
 
-
-
-
 #np = 3
 #ndofs = 4
 #

test/FieldsTests/runtests.jl

@@ -20,4 +20,6 @@ using Test
 
 @testset "InverseFields" begin include("InverseFieldsTests.jl") end
 
+@testset "InverseFields" begin include("InverseFieldsTests.jl") end
+
 end

test/MultiFieldTests/MultiFieldCellFieldsTests.jl

@@ -9,6 +9,8 @@ using Gridap.Fields
 using Gridap.CellData
 using Gridap.MultiField
 using Gridap.TensorValues
+using Random
+using StaticArrays
 using Test
 
 domain = (0,1,0,1)
@@ -116,6 +118,50 @@ cellmat1_Γ = integrate(((n⋅dv.⁺)-dq.⁻)*((n⋅du.⁺)+dp.⁻),quad_Γ)
 cellmat2_Γ = integrate((n⋅dv.⁺)*(n⋅du.⁺)+(n⋅dv.⁺)*dp.⁻-dq.⁻*(n⋅du.⁺)-dq.⁻*dp.⁻,quad_Γ)
 test_array(cellmat1_Γ,cellmat2_Γ,≈)
 
+# Test function evaluation
+
+# Set reproducible random number seed
+Random.seed!(0)
+@testset "evaluating functions" for D in 1:5
+    xmin = 0
+    xmax = 1
+    domain = repeat([xmin, xmax], D)
+    ncells = 3
+    partition = repeat([ncells], D)
+    model = CartesianDiscreteModel(domain, partition)
+    # TODO: test both with and without this
+    model = simplexify(model)
+
+    order = 2
+    reffe = ReferenceFE(lagrangian, Float64, order)
+    V = FESpace(model, reffe)
+
+    coeff0 = rand(Float64)
+    coeffs = rand(SVector{D,Float64})
+    f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0
+    # TODO: use this mechanism instead to project
+    # Francesc Verdugo @fverdugo 13:11
+    # a(u,v) = ∫( u*v )dΩ
+    # l(v) = a(f,v)
+    # Solve a fe problem with this weak form
+    # See also tutorial 10, "Isotropic damage model", section "L2
+    # projection", function "project"
+    fh = interpolate_everywhere(f, V)
+    fhcache = return_cache(fh, VectorValue(zeros(D)...))
+
+    xs = [VectorValue(rand(D)...) for i in 1:10]
+    for x in xs
+        x = VectorValue(rand(D)...)
+        fx = f(x)
+        fhx = evaluate!(fhcache, fh, x)
+        @test fhx ≈ fx
+    end
+    fhxs = fh(xs)
+    @test fhxs ≈ f.(xs)
+end
+
+
+
 #a = cellmat_Γ
 #using BenchmarkTools
 #cache = array_cache(a)