gridap · santiagobadia · Jun 2, 2021 · Jan 26, 2021 · Jan 26, 2021 · Jan 4, 2021
diff --git a/NEWS.md b/NEWS.md
@@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Experimental support for mixed dimensional PDEs. Since PR [#567](https://github.com/gridap/Gridap.jl/pull/567).
 - Added `get_cell_dof_basis(model,cell_reffes,::Conformity)` and `get_cell_shapefuns(model,cell_reffes,::Conformity)`. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579).
 - Implemented `get_cell_dof_basis` and `get_cell_shapefuns` for global RT FE spaces in a new file `DivConformingFESpaces.jl`. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579).
+- Added support to allow evaluation of FE functions at arbitrary points. Since PR [#523](https://github.com/gridap/Gridap.jl/pull/523).
+- Implemented `compute_cell_points_from_vector_of_points` to build `CellPoint` from a vector of points. Since PR [#523](https://github.com/gridap/Gridap.jl/pull/523).
 
 ### Changed
 - Major refactoring in the handling of blocks (e.g. in multi-field and skeleton terms). The new code follows a much more simple approach based in the new type `ArrayBlock`. Since PR [#583](https://github.com/gridap/Gridap.jl/pull/583).
@@ -28,7 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - All code associated with with `BlockArrayCoo`. Since PR [#583](https://github.com/gridap/Gridap.jl/pull/583).
 - Module `Gridap.Integration` has been deleted and its contents have been merged into `Gridap.ReferenceFEs` module.
 - Types `SparseMatrixCSR` and `SymSparseMatrixCSR` have been moved to the registered package [`SparseMatricesCSR`](https://github.com/gridap/SparseMatricesCSR.jl). To use them simply add `SparseMatricesCSR` into your environment and type `using SparseMatricesCSR`. Since  Since PR [#568](https://github.com/gridap/Gridap.jl/pull/568).
-- Removed `PushForwardMap` and all code depending upon it. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579). 
+- Removed `PushForwardMap` and all code depending upon it. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579).
 
 ## [0.15.5] - 2021-05-25
 
@@ -43,7 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [0.15.3] - 2021-03-16
 
-### Added 
+### Added
 - `get_cell_map` now returns array of `AffineMap` for linear grids of simplices. Needed to compute Laplacian operator, inverse maps etc. Since PR [#553](https://github.com/gridap/Gridap.jl/pull/553).
 
 ### Fixed
@@ -83,7 +85,7 @@ This version is a major (backwards-incompatible) refactoring of the project whic
  - Bug-fix for 32-bit Julia: Replace all occurences of Int64 by Int. Since PR [#445](https://github.com/gridap/Gridap.jl/pull/445).
  - Bug-fix for 32-bit Julia. Using inttype=Int keyword argument for JSON parsing. Since PR [#456](https://github.com/gridap/Gridap.jl/pull/456).
 
-## [0.14.1] - 2020-09-17 
+## [0.14.1] - 2020-09-17
 
 ### Added
  - Added VectorWithEntryInserted and VectorWithEntryRemoved. Since PR [#401](https://github.com/gridap/Gridap.jl/pull/401/).
@@ -102,7 +104,7 @@ This version is a major (backwards-incompatible) refactoring of the project whic
  - First and second argument switch in `update_state_variables!` in order to have function-first style. Since PR [#376](https://github.com/gridap/Gridap.jl/pull/376/).
  - Table struct has been generalized such that data and ptrs arrays can be of an arbitrary type extending AbstractArray. Since PR [#310](https://github.com/gridap/Gridap.jl/pull/310/)
  - `interpolate, interpolate!, interpolate_dirichlet...` switched argument order to function first style. For instance `interpolate(u, V)` instead of `interpolate(V, u)`
- 
+
 ### Added
  - Allowing the construction of an `HomogeneousTrialFESpace` from a `TrialFESpace`. Since PR [#384](https://github.com/gridap/Gridap.jl/pull/384).
  - Support for automatic differentiation of residuals and Jacobians in multi-field computations since PR [#383](https://github.com/gridap/Gridap.jl/pull/383/).
@@ -113,7 +115,7 @@ This version is a major (backwards-incompatible) refactoring of the project whic
 
 ## [0.13.4] - 2020-08-23
 
-### Added 
+### Added
   - New `FilteredCellArray` since PR [#372](https://github.com/gridap/Gridap.jl/pull/372/).
 
 ## [0.13.3] - 2020-08-12

diff --git a/Project.toml b/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"

diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl
@@ -10,16 +10,20 @@ 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
 using Gridap.Geometry
 
+using NearestNeighbors
+using StaticArrays
+
 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
@@ -57,6 +61,10 @@ export CellDof
 export get_normal_vector
 export get_cell_measure
 
+export make_inverse_table
+export point_to_cell!
+export compute_cell_points_from_vector_of_points
+
 export DomainContribution
 export num_domains
 export get_domains

diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl
@@ -1,4 +1,3 @@
-
 """
 A single point or an array of points on the cells of a Triangulation
 CellField objects can be evaluated efficiently at CellPoint instances.
@@ -208,20 +207,182 @@ get_data(f::GenericCellField) = f.cell_field
 get_triangulation(f::GenericCellField) = f.trian
 DomainStyle(::Type{GenericCellField{DS}}) where DS = DS()
 
-# Evaluation of CellFields
 
-(a::CellField)(x) = evaluate(a,x)
+"""
+   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)
+  cache1 = _point_to_cell_cache(trian)
+
+  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(trian::Triangulation)
+  topo = GridTopology(trian)
+  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
+end
+
+function point_to_cell!(cache, 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)
-  @notimplemented """\n
-  Evaluation of a CellField at a given Point is not implemented yet.
+  cache1,cache2 = cache
+  cell_f_cache, f_cache, cell_f, f₀ = cache2
+  @check f === f₀ "Wrong cache"
+
+  cell = point_to_cell!(cache1, 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)
 
-  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.
-  """
+  return j2is,i2lis
 end
 
+# Efficient version:
+function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point})
+  cache1,cache2 = cache
+  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,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
+
+function compute_cell_points_from_vector_of_points(xs::AbstractVector{<:Point}, trian::Triangulation, domain_style::PhysicalDomain)
+    cache1 = _point_to_cell_cache(trian)
+    x_to_cell(x) = point_to_cell!(cache1, x)
+    point_to_cell = map(x_to_cell, xs)
+    ncells = num_cells(trian)
+    cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells)
+    cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points)
+    cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain())
+end
+
+(a::CellField)(x) = evaluate(a,x)
+
 function evaluate!(cache,f::CellField,x::CellPoint)
   _f, _x = _to_common_domain(f,x)
   cell_field = get_data(_f)
@@ -334,11 +495,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

diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl
@@ -15,12 +15,13 @@ 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 FillArrays
 using SparseArrays
 using LinearAlgebra
-import BlockArrays
+using BlockArrays
 
 export num_fields
 export compute_field_offsets

diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl
@@ -1,4 +1,3 @@
-
 struct MultiFieldCellField{DS<:DomainStyle} <: CellField
   single_fields::Vector{<:CellField}
   trian::Triangulation