allow a not_in_domain_value in interpolate!

Author: chmerdon

src/cellfinder.jl

@@ -95,7 +95,7 @@ Upon return, xref contains the barycentric coordinates of the point in the seque
 !!! warning
     Currently implemented for simplex grids only.
 """
-function gFindLocal!(xref, CF::CellFinder{Tv, Ti}, x; icellstart = 1, trybrute = true, eps = 1.0e-14) where {Tv, Ti}
+function gFindLocal!(xref, CF::CellFinder{Tv, Ti}, x; icellstart = 1, stay_in_cell = false, trybrute = true, eps = 1.0e-14) where {Tv, Ti}
     # works for convex domains and simplices only !
     xCellFaces::Adjacency{Ti} = CF.xCellFaces
     xFaceCells::Adjacency{Ti} = CF.xFaceCells
@@ -131,6 +131,11 @@ function gFindLocal!(xref, CF::CellFinder{Tv, Ti}, x; icellstart = 1, trybrute =
             cx[j] = x[j] - L2Gb[j]
         end
         mapderiv!(invA, L2G, xref)
+
+        if stay_in_cell
+            return icell
+        end
+
         fill!(xreftest, 0)
         for j in 1:length(cx), k in 1:length(cx)
             xreftest[k] += invA[j, k] * cx[j]
@@ -258,7 +263,7 @@ end
 
 Mutating form of [`interpolate`](@ref)
 """
-function interpolate!(u_to::AbstractArray, grid_to, u_from::AbstractArray, grid_from; eps = 1.0e-14, trybrute = true)
+function interpolate!(u_to::AbstractArray, grid_to, u_from::AbstractArray, grid_from; eps = 1.0e-14, not_in_domain_value = NaN, trybrute = true)
     shuffle = [[2, 1], [3, 1, 2], [4, 1, 2, 3]]
 
     update!(u_to::AbstractVector, inode_to, λ, u_from::AbstractVector, inode_from) = u_to[inode_to] += λ * u_from[inode_from]
@@ -286,10 +291,13 @@ function interpolate!(u_to::AbstractArray, grid_to, u_from::AbstractArray, grid_
     icellstart = 1
     for inode_to in 1:nnodes_to
         @views icell_from = gFindLocal!(λ, cf, coord[:, inode_to]; icellstart, eps, trybrute)
-        @assert icell_from > 0
-        for i in 1:(dim + 1)
-            inode_from = cn_from[i, icell_from]
-            update!(u_to, inode_to, λ_shuffle[i], u_from, inode_from)
+        if icell_from <= 0
+            u_to[inode_to] = not_in_domain_value
+        else
+            for i in 1:(dim + 1)
+                inode_from = cn_from[i, icell_from]
+                update!(u_to, inode_to, λ_shuffle[i], u_from, inode_from)
+            end
         end
         icell_start = icell_from
     end