Allow evaluating FE functions at arbitrary points

Project.toml

@@ -6,6 +6,7 @@ version = "0.15.1"
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
+Bernstein = "e9b0fb4c-9cb7-4f61-9c14-701a41c684d7"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -18,6 +19,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"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -27,6 +29,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [compat]
 AbstractTrees = "0.3.3"
 BSON = "0.2.5"
+Bernstein = "1.0.1"
 BlockArrays = "0.12.12, 0.13"
 Combinatorics = "1.0.0"
 DocStringExtensions = "0.8.1"
@@ -38,6 +41,7 @@ JLD2 = "0.1.11, 0.3"
 JSON = "0.21.0"
 LineSearches = "7.0.1"
 NLsolve = "4.3.0"
+NearestNeighbors = "0.4.8"
 QuadGK = "2.3.1, 2.4"
 StaticArrays = "0.12.1, 1.0"
 WriteVTK = "1.7, 1.8"

src/CellData/CellData.jl

@@ -8,6 +8,9 @@ module CellData
 using Test
 using DocStringExtensions
 using FillArrays
+using Bernstein
+using NearestNeighbors
+using StaticArrays
 
 using Gridap.Helpers
 using Gridap.Arrays
@@ -21,6 +24,7 @@ using Gridap.Integration
 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

@@ -113,7 +113,7 @@ function change_domain(a::CellField,::ReferenceDomain,::PhysicalDomain)
   trian = get_triangulation(a)
   cell_map = get_cell_map(trian)
   cell_invmap = lazy_map(inverse_map,cell_map)
-  cell_field_ref = get_data(cell_field)
+  cell_field_ref = get_data(a)
   cell_field_phys = lazy_map(Broadcasting(∘),cell_field_ref,cell_invmap)
   GenericCellField(cell_field_phys,trian,PhysicalDomain())
 end
@@ -204,14 +204,97 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS()
 
 (a::CellField)(x) = evaluate(a,x)
 
+struct CellFieldCache
+  kdtree::KDTree{SVector{D,T},Euclidean,T} where {D,T}
+end
+
+function return_cache(f::CellField,x::Point)
+  # Create k-d tree for mesh
+  trian = get_triangulation(f)
+  node_coordinates = get_node_coordinates(trian)
+  @assert !isempty(node_coordinates)
+  sample_coord = node_coordinates[1]
+  D = length(sample_coord)
+  T = eltype(sample_coord)
+  @assert D > 0
+  kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates))
+  kdtree::KDTree{SVector{D,T},Euclidean,T}
+  return CellFieldCache(kdtree)
+end
+
 function evaluate!(cache,f::CellField,x::Point)
-  @notimplemented """\n
-  Evaluation of a CellField at a given Point is not implemented yet.
+  cache::CellFieldCache
+
+  # Get k-d tree for mesh from cache
+  trian = get_triangulation(f)
+  node_coordinates = get_node_coordinates(trian)
+  @assert !isempty(node_coordinates)
+  sample_coord = node_coordinates[1]
+  D = length(sample_coord)
+  T = eltype(sample_coord)
+  @assert D > 0
+  kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T}
+
+  # Find neighbouring cells
+  id,dist = nn(kdtree, SVector(x.data))
+  cells = Int[]
+  # TODO: Use better algorithm
+  # Francesc Verdugo @fverdugo 14:07
+  # You have to start with a DiscreteModel, then
+  # topo=get_grid_topology(model)
+  # vertex_to_cells = get_faces(topo,0,num_cell_dims(model))
+  for (c,ids) in enumerate(trian.cell_node_ids)
+    id ∈ ids && push!(cells, c)
+  end
+  @assert !isempty(cells)
+
+  # Find polytope
+  reffes = trian.reffes
+  @assert length(reffes) == 1
+  reffe = reffes[1].reffe
+  polytope = reffe.polytope
+  extrusion = polytope.extrusion
+  @assert length(extrusion) == D
+
+  # Loop over neighbouring cells and calculate distances
+  dists = T[]
+
+  if all(e == TET_AXIS for e in extrusion)
+    # Calculate barycentric coordinates for a point x in the given cell
+    function bary(cell, x)
+      local ids = trian.cell_node_ids[cell]
+      @assert length(ids) == D+1 "simplex has correct number of vertices"
+      local ys = zero(SVector{D+1,SVector{D,T}})
+      for (n,id) in enumerate(ids)
+        local y = SVector(node_coordinates[id].data)
+        ys = setindex(ys, y, n)
+      end
+      local λ = cartesian2barycentric(ys, SVector(x.data))
+      return λ
+    end
 
-  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.
-  """
+    for c in cells
+      λ = bary(c, x)
+      # Negative distances are outside, positive distance are inside the simplex
+      dist = min(minimum(λ), 1 - maximum(λ))
+      push!(dists, dist)
+    end
+
+    dist, m = findmax(dists)
+    cell = cells[m]
+    λ = bary(cell, x)
+
+    # Ensure the point is inside one of the cells, up to round-off errors
+    @assert dist ≥ -1000 * eps(T) "Point is not inside a cell"
+
+    # TODO: Actually evaluate basis functions
+    values = f.cell_dof_values[cell]
+    fx = sum(λ[n] * values[n] for n in 1:D+1)
+  else
+    @assert false "Unsupported polytope"
+  end
+
+  return fx
 end
 
 function evaluate!(cache,f::CellField,x::CellPoint)
@@ -324,11 +407,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/Fields/AffineMaps.jl

@@ -69,3 +69,11 @@ function lazy_map(
   cell_∇∇bt = lazy_map(transpose,cell_∇∇b)
   cell_∇∇bt
 end
+
+function inverse_map(f::AffineMap)
+  Jt = f.gradient
+  y0 = f.origin
+  invJt = inv(Jt)
+  x0 = -y0⋅invJt
+  AffineMap(invJt,x0)
+end