Error with the interpolation of 3D model

[dryudryu]

I've been struggling to evaluate a 3D solution at some points. Although I am sampling at a position within my mesh range, I got an error saying that I'm attempting to sample outside the active cell. Please refer to the following script and error message.
Here is a link for the msh file ("SimpleBox3D.msh"). https://gist.github.com/dryudryu/0545a5d1ea5270a84f9c4430da41866c#file-simplebox3d-msh
If this doesn't work, one can generate .msh by running this .jl
https://gist.github.com/dryudryu/ab7aaec1b7df7f39d0413e23055ff67e

Thanks!

using Gridap
using GridapGmsh
using Gridap.Fields
using Gridap.Geometry

model = GmshDiscreteModel("../models/SimpleBox3D.msh") # range [-2:2, -3:3, -0.5:0.5]

u(x) = x

order = 0
reffe = ReferenceFE(nedelec, Float64, order)

V = TestFESpace(model, reffe, vector_type = Vector{ComplexF64})

Ω = Triangulation(model)
uh3d = interpolate(u,V)

const cache = return_cache(uh3d,Point(0.,0.,0.))
f((x,y)) = evaluate!(cache,uh3d,Point(x,y,0))
writevtk(Ω, "sol2d", nsubcells=4, cellfields=["uh3d"=>f])

The following is error

AssertionError: Point (0.4521502328465496, -2.6749148894773853, 0.0) is not inside any active cell

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/Gridap/TyUsh/src/Helpers/Macros.jl:60 [inlined]
  [2] _point_to_cell!(cache::Tuple{Gridap.CellData.KDTreeSearch, NearestNeighbors.KDTree{StaticArraysCore.SVector{3, Float64}, Distances.Euclidean, Float64}, Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}, Vector{Int8}, Vector{Gridap.ReferenceFEs.ExtrusionPolytope{3}}, Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(affine_map), 1, Tuple{Base.OneTo{Int64}}}, AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}, Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}, x::VectorValue{3, Float64})
    @ Gridap.CellData ~/.julia/packages/Gridap/TyUsh/src/CellData/CellFields.jl:301
  [3] evaluate!(cache::Tuple{Tuple{Gridap.CellData.KDTreeSearch, NearestNeighbors.KDTree{StaticArraysCore.SVector{3, Float64}, Distances.Euclidean, Float64}, Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}, Vector{Int8}, Vector{Gridap.ReferenceFEs.ExtrusionPolytope{3}}, Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(affine_map), 1, Tuple{Base.OneTo{Int64}}}, AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}, Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{ComplexF64, Vector{ComplexF64}}, Tuple{Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{ComplexF64}}}, Tuple{Tuple{Nothing, Nothing, Tuple{Nothing, Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int64, Vector{Int64}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, TensorValue{3, 3, Float64, 9}}}}}, Gridap.Arrays.IndexItemPair{Int64, ConstantField{TensorValue{3, 3, Float64, 9}}}}}}, Gridap.Arrays.IndexItemPair{Int64, Gridap.Fields.BroadcastOpFieldArray{Gridap.ReferenceFEs.CoVariantPiolaMap, Gridap.Fields.OperationField{Gridap.ReferenceFEs.CoVariantPiolaMap, Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}}}}}, Gridap.Arrays.IndexItemPair{Int64, Gridap.Fields.LinearCombinationField{Vector{ComplexF64}, Gridap.Fields.BroadcastOpFieldArray{Gridap.ReferenceFEs.CoVariantPiolaMap, Gridap.Fields.OperationField{Gridap.ReferenceFEs.CoVariantPiolaMap, Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}}}}, Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int64, Vector{Int64}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, TensorValue{3, 3, Float64, 9}}}, Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int64, Vector{Int64}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, VectorValue{3, Float64}}}}}, Gridap.Arrays.IndexItemPair{Int64, AffineMap{3, 3, Float64, 9}}}}}, Gridap.Arrays.IndexItemPair{Int64, AffineMap{3, 3, Float64, 9}}}}}, Gridap.Arrays.IndexItemPair{Int64, Gridap.Fields.OperationField{Gridap.Fields.LinearCombinationField{Vector{ComplexF64}, Gridap.Fields.BroadcastOpFieldArray{Gridap.ReferenceFEs.CoVariantPiolaMap, Gridap.Fields.OperationField{Gridap.ReferenceFEs.CoVariantPiolaMap, Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}}, Tuple{AffineMap{3, 3, Float64, 9}}}}}, Tuple{Tuple{Tuple{Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Tuple{Tuple{Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedMatrix{VectorValue{3, Float64}, Matrix{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedMatrix{VectorValue{3, Float64}, Matrix{VectorValue{3, Float64}}}, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Gridap.Arrays.CachedMatrix{Float64, Matrix{Float64}}}, Gridap.Polynomials.MonomialBasis{3, VectorValue{3, Float64}}}, Vector{VectorValue{3, Float64}}}, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}}, Nothing}}, Nothing}, Tuple{Nothing}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{typeof(∘)}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.OperationField{Gridap.Fields.LinearCombinationField{Vector{ComplexF64}, Gridap.Fields.BroadcastOpFieldArray{Gridap.ReferenceFEs.CoVariantPiolaMap, Gridap.Fields.OperationField{Gridap.ReferenceFEs.CoVariantPiolaMap, Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}}, Tuple{AffineMap{3, 3, Float64, 9}}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(linear_combination), 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.LinearCombinationField{Vector{ComplexF64}, Gridap.Fields.BroadcastOpFieldArray{Gridap.ReferenceFEs.CoVariantPiolaMap, Gridap.Fields.OperationField{Gridap.ReferenceFEs.CoVariantPiolaMap, Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Gridap.Arrays.PosNegReindex{Vector{ComplexF64}, Vector{ComplexF64}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{ComplexF64}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Operation{Gridap.ReferenceFEs.CoVariantPiolaMap}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.BroadcastOpFieldArray{Gridap.ReferenceFEs.CoVariantPiolaMap, Gridap.Fields.OperationField{Gridap.ReferenceFEs.CoVariantPiolaMap, Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, ConstantField{TensorValue{3, 3, Float64, 9}}}}, 1, Tuple{Gridap.Arrays.CompressedArray{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}, 1, Vector{Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.NedelecPrebasisOnSimplex{3}}}, Vector{Int8}}, Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(constant_field), 1, Tuple{Base.OneTo{Int64}}}, ConstantField{TensorValue{3, 3, Float64, 9}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}}}}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(inverse_map), 1, Tuple{Base.OneTo{Int64}}}, AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(affine_map), 1, Tuple{Base.OneTo{Int64}}}, AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}}}}}, Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}}}, f::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, x::VectorValue{3, Float64})
    @ Gridap.CellData ~/.julia/packages/Gridap/TyUsh/src/CellData/CellFields.jl:309
  [4] f(::VectorValue{3, Float64})
    @ Main ./In[8]:18
  [5] evaluate!
    @ ~/.julia/packages/Gridap/TyUsh/src/Arrays/Maps.jl:93 [inlined]
  [6] evaluate!
    @ ~/.julia/packages/Gridap/TyUsh/src/Fields/FieldsInterfaces.jl:175 [inlined]
  [7] evaluate!(c::Tuple{Gridap.Arrays.CachedVector{VectorValue{3, ComplexF64}, Vector{VectorValue{3, ComplexF64}}}, Nothing}, f::GenericField{typeof(f)}, x::Vector{VectorValue{3, Float64}})
    @ Gridap.Fields ~/.julia/packages/Gridap/TyUsh/src/Fields/FieldsInterfaces.jl:155
  [8] _getindex_and_call!
    @ ~/.julia/packages/Gridap/TyUsh/src/Arrays/LazyArrays.jl:327 [inlined]
  [9] getindex!
    @ ~/.julia/packages/Gridap/TyUsh/src/Arrays/LazyArrays.jl:214 [inlined]
 [10] _prepare_node_to_coords_fill!(node_to_coords::Vector{VectorValue{3, ComplexF64}}, cell_to_offset::Vector{Int64}, cache::Tuple{Tuple{Nothing, Tuple{Gridap.Arrays.CachedVector{VectorValue{3, ComplexF64}, Vector{VectorValue{3, ComplexF64}}}, Nothing}, Tuple{Tuple{Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int64, Vector{Int64}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, TensorValue{3, 3, Float64, 9}}}, Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int64, Vector{Int64}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, VectorValue{3, Float64}}}}}, Gridap.Arrays.IndexItemPair{Int64, AffineMap{3, 3, Float64, 9}}}, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Nothing}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, ComplexF64}}}}, cell_to_points::Gridap.Arrays.LazyArray{FillArrays.Fill{GenericField{typeof(f)}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, ComplexF64}}, 1, Tuple{Gridap.Arrays.LazyArray{Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(affine_map), 1, Tuple{Base.OneTo{Int64}}}, AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}}})
    @ Gridap.Visualization ~/.julia/packages/Gridap/TyUsh/src/Visualization/VisualizationData.jl:65
 [11] _prepare_node_to_coords(cell_to_points::Gridap.Arrays.LazyArray{FillArrays.Fill{GenericField{typeof(f)}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, ComplexF64}}, 1, Tuple{Gridap.Arrays.LazyArray{Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(affine_map), 1, Tuple{Base.OneTo{Int64}}}, AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}}})
    @ Gridap.Visualization ~/.julia/packages/Gridap/TyUsh/src/Visualization/VisualizationData.jl:57
 [12] _prepare_pdata(trian::BodyFittedTriangulation{3, 3, UnstructuredDiscreteModel{3, 3, Float64, Oriented}, UnstructuredGrid{3, 3, Float64, Oriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}}, cellfields::Vector{Pair{String, typeof(f)}}, samplingpoints::Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}})
    @ Gridap.Visualization ~/.julia/packages/Gridap/TyUsh/src/Visualization/VisualizationData.jl:141
 [13] visualization_data(trian::BodyFittedTriangulation{3, 3, UnstructuredDiscreteModel{3, 3, Float64, Oriented}, UnstructuredGrid{3, 3, Float64, Oriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}}, filebase::String; order::Int64, nsubcells::Int64, celldata::Dict{Any, Any}, cellfields::Vector{Pair{String, typeof(f)}})
    @ Gridap.Visualization ~/.julia/packages/Gridap/TyUsh/src/Visualization/VisualizationData.jl:47
 [14] writevtk(::BodyFittedTriangulation{3, 3, UnstructuredDiscreteModel{3, 3, Float64, Oriented}, UnstructuredGrid{3, 3, Float64, Oriented, Nothing}, Gridap.Arrays.IdentityVector{Int64}}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:nsubcells, :cellfields), Tuple{Int64, Vector{Pair{String, typeof(f)}}}}})
    @ Gridap.Visualization ~/.julia/packages/Gridap/TyUsh/src/Visualization/Vtk.jl:4
 [15] top-level scope
    @ In[8]:19
 [16] eval
    @ ./boot.jl:373 [inlined]
 [17] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196

[stevengj]

I can reproduce. You don't even need the cache or the writevtk to get an error. Simply evaluating uh3d at an "unlucky" point fails:

julia> uh3d(Point(0.4521502328465496, -2.6749148894773853, 0.0))
ERROR: AssertionError: Point (0.4521502328465496, -2.6749148894773853, 0.0) is not inside any active cell

even though very nearby points succeed (changing z from 0.0 to ±0.01):

julia> uh3d(Point(0.4521502328465496, -2.6749148894773853, 0.01))
VectorValue{3, ComplexF64}(0.42891370112403415 + 0.0im, -2.63640975080864 + 0.0im, 0.016364254492839403 + 0.0im)

julia> uh3d(Point(0.4521502328465496, -2.6749148894773853, -0.01))
VectorValue{3, ComplexF64}(0.3982795100619132 + 0.0im, -2.6202275972552007 + 0.0im, -0.004346949958980106 + 0.0im)

On the other hand, it also fails for z = ±0.001, so it's not simply some rounding error that is happening exactly at z=0.0.

[fverdugo]

Hi @Balaje can you take a quick look to this?

[oriolcg]

Hi @dryudryu and @stevengj , try to use more nearest vertices in the search algorithm. In unstructured meshes it is common to find points whose nearest vertex is not a vertex of the element containing such point. You can modify that by passing an argument to the KDTreeSearch and using interpolable function:

uh3d = interpolate(u,V)
search_method = KDTreeSearch(num_nearest_vertices=5)
uh3d_ = interpolable(uh3d;search_method=search_method)
uh3d_(Point(x,y,z))

Please, let me know if that solves the problem.

[fverdugo]

@oriolcg thanks for pointing this out!

I can confirm that it works:

using Gridap
using GridapGmsh
using Gridap.Fields
using Gridap.Geometry
using Gridap.CellData
using Gridap.FESpaces

model = GmshDiscreteModel("SimpleBox3D.msh") # range [-2:2, -3:3, -0.5:0.5]
u(x) = x
order = 0
reffe = ReferenceFE(nedelec, Float64, order)
V = TestFESpace(model, reffe, vector_type = Vector{ComplexF64})
Ω = Triangulation(model)
uh3d = interpolate(u,V)
search_method = KDTreeSearch(num_nearest_vertices=5)
uh3d_ = Interpolable(uh3d;searchmethod=search_method)
uh3d_(Point(0.4521502328465496, -2.6749148894773853, 0.0))

[Balaje]

Hi @fverdugo ,

Sorry for the late reply. I tried to run the code and load the mesh file using GridapGmsh and I get the following error:

Info    : Reading 'SimpleBox3D.msh'...
Info    : 45 entities
Info    : 3438 nodes
Info    : 14977 elements
Info    : Done reading 'SimpleBox3D.msh'
ERROR: LoadError: BoundsError: attempt to access 2486-element Vector{Int32} at index [0]
Stacktrace:
  [1] getindex
    @ ./array.jl:801 [inlined]
  [2] _face_to_cells_fill!(face_to_cells_data::Vector{Int64}, face_to_cells_ptrs::Vector{Int32}, cell_to_faces_data::Vector{Int64}, cell_to_faces_ptrs::Vector{Int32})
    @ Gridap.Geometry ~/.julia/packages/Gridap/uOmS3/src/Geometry/GridTopologies.jl:902
  [3] _face_to_cells(cell_to_faces_data::Vector{Int64}, cell_to_faces_ptrs::Vector{Int32}, nfaces::Int64)
    @ Gridap.Geometry ~/.julia/packages/Gridap/uOmS3/src/Geometry/GridTopologies.jl:875
  [4] generate_cells_around
    @ ~/.julia/packages/Gridap/uOmS3/src/Geometry/GridTopologies.jl:563 [inlined]
  [5] UnstructuredGridTopology(vertex_coordinates::Vector{VectorValue{3, Float64}}, cell_vertices::Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}, cell_type::Vector{Int8}, polytopes::Vector{Gridap.ReferenceFEs.ExtrusionPolytope{3}}, orientation_style::Oriented)
    @ Gridap.Geometry ~/.julia/packages/Gridap/uOmS3/src/Geometry/UnstructuredGridTopologies.jl:36
  [6] _generate_grid_topology_from_grid(grid::UnstructuredGrid{3, 3, Float64, Oriented, Nothing}, cell_to_vertices::Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}, vertex_to_node::Vector{Int64})
    @ Gridap.Geometry ~/.julia/packages/Gridap/uOmS3/src/Geometry/UnstructuredGridTopologies.jl:148
  [7] UnstructuredGridTopology(grid::UnstructuredGrid{3, 3, Float64, Oriented, Nothing}, cell_to_vertices::Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}, vertex_to_node::Vector{Int64})
    @ Gridap.Geometry ~/.julia/packages/Gridap/uOmS3/src/Geometry/UnstructuredGridTopologies.jl:131
  [8] GmshDiscreteModel(mshfile::String; renumber::Bool)
    @ GridapGmsh ~/.julia/packages/GridapGmsh/Q8ZwW/src/GmshDiscreteModels.jl:23
  [9] GmshDiscreteModel(mshfile::String)
    @ GridapGmsh ~/.julia/packages/GridapGmsh/Q8ZwW/src/GmshDiscreteModels.jl:7
 [10] top-level scope
    @ ~/Repositories/test_ip.jl:8
 [11] include(fname::String)
    @ Base.MainInclude ./client.jl:444
 [12] top-level scope
    @ REPL[1]:1

I get this error when I use the SimpleBox3D.msh file directly. The error remains even if I run the code provided by @dryudryu to generate the mesh. However when I use the cube.msh file from the GridapGmsh tests, it seems to run okay.

[fverdugo]

Hi @Balaje! Make sure that you use the latest version of GridapGmsh (I think it is at least 0.5.1) and also the latest version of gmsh provided by BinaryBuilder. The model SimpleBox3D.msh provably has periodic boundary conditions, and I found and fixed some errors related to this recently.

[Balaje]

Hi @fverdugo ,

Yes! I can confirm that it works too.

using Gridap
using GridapGmsh
using Gridap.Fields
using Gridap.Geometry
using Gridap.CellData
using Gridap.FESpaces

#model = GmshDiscreteModel("./SimpleBox3D.msh") # range [-2:2, -3:3, -0.5:0.5]
model = GmshDiscreteModel("./SimpleBox3D_1.msh") # range [-2:2, -3:3, -0.5:0.5]
u(x) = x
order = 0
reffe = ReferenceFE(nedelec, Float64, order)
V = TestFESpace(model, reffe, vector_type = Vector{ComplexF64})
Ω = Triangulation(model)
uh3d = interpolate(u,V)
pt = Point(0.4521502328465496, -2.6749148894773853, 0.0)
search_method = KDTreeSearch(num_nearest_vertices=5)
uh3d_ = Interpolable(uh3d;searchmethod=search_method)
@show uh3d_(pt) # Works!
@show uh3d(pt) # Error...

Thanks.

[fverdugo]

Should we perhaps inclrease the default value of num_nearest_vertices ? Or precompute a value that is guaranteed to work?

[Balaje]

I can't think of an optimal solution from the top of my mind. I have a couple of thoughts though.

I think @oriolcg presently gets the value of the num_nearest_vertices and uses it in the _point_to_cell! function

Gridap.jl/src/CellData/CellFields.jl

Line 263 in 8f842f8

for (id,dist) in zip(knn(kdtree, SVector(Tuple(x)), searchmethod.num_nearest_vertices, true)...)

How about starting with a low value of num_nearest_vertices and then increasing the value after checking if the search failed? We could then keep the computational cost to a minimum. We could do it up to a certain level (say 5 or 6 nearest neighbours), and if it did not succeed, then conclude that the point lies outside the triangulation.

We could have a third option of overriding the above feature and letting the user input a higher value if needed.

[oriolcg]

I think this is how it is implemented now. We just need to increase the default value of num_nearest_vertices.

[Balaje]

@oriolcg As I understand, the problem is choosing the value of num_nearest_vertices accordingly so that the algorithm finds the correct cell. We could have num_nearest_vertices = 20 or something so that enough points are found, but wouldn't that make running knn(kdtree, point, 20, true) expensive to find the vertices in the mesh? It may not be too expensive, as I think the tree is already in the cache. But I don't know.

[stevengj]

I've implemented this sort of thing before, and the right thing to do is to make a KD tree that searches elements, not vertices.

That is, you recursively halve space by planes, as in a normal KD tree, but for each separating plane you then make a list of elements that lies on either side of the plane. And for any element that intersects the separating plane, you put it in both nodes of the tree.

For example, here is a version that I wrote for a BEM code: https://github.com/HomerReid/scuff-em/blob/9c6d0cb7695463af803dee8d04cdae939740cdcc/libs/libscuff/PointInObject.cc

• One trick is that I usually first preprocess the elements/objects by computing a Cartesian bounding box for each object (i.e. aligned with the xyz axes). The reason is that it's extremely quick to check whether a point is in a Cartesian box or whether a Cartesian plane intersects the box. You only have to do the more expensive point-in-polyhedron test if the the answer is true (if the point is not in the bounding box then it is definitely not in the polyhedron).

[stevengj]

I forgot that I've also implemented kd-tree searching of objects in Julia https://github.com/stevengj/GeometryPrimitives.jl/blob/master/src/util/kdtree.jl