Fix PointEvalHandler on true subdomains

src/Dofs/DofHandler.jl

@@ -108,6 +108,14 @@ mutable struct DofHandler{dim, G <: AbstractGrid{dim}} <: AbstractDofHandler
     ndofs::Int
 end
 
+function is_defined_on_full_domain(dh::DofHandler)
+    num_cells = 0
+    for sdh in dh.subdofhandlers
+        num_cells += length(sdh.cellset)
+    end
+    return getncells(get_grid(dh)) == num_cells
+end
+
 """
     DofHandler(grid::Grid)
 

src/PointEvalHandler.jl

@@ -5,7 +5,7 @@
 end
 
 """
-    PointEvalHandler(grid::Grid, points::AbstractVector{Vec{dim,T}}; kwargs...) where {dim, T}
+    PointEvalHandler(grid::Grid, points::AbstractVector{Vec{dim,T}} [, set = subdomains::Vector{OrderedSet{Int}}] ; kwargs...) where {dim, T}
 
 The `PointEvalHandler` can be used for function evaluation in *arbitrary points* in the
 domain -- not just in quadrature points or nodes.
@@ -33,6 +33,7 @@
     grid::G
     cells::Vector{Union{Nothing, Int}}
     local_coords::Vector{Union{Nothing, Vec{1, T}, Vec{2, T}, Vec{3, T}}}
+    constructed_on_subdomain::Bool
 end
 
 function Base.show(io::IO, ::MIME"text/plain", ph::PointEvalHandler)
@@ -58,7 +59,13 @@
 function PointEvalHandler(grid::AbstractGrid{dim}, points::AbstractVector{Vec{dim, T}}; search_nneighbors = 3, warn::Bool = true, strategy = NewtonLineSearchPointFinder()) where {dim, T}
     node_cell_dicts = _get_node_cell_map(grid)
     cells, local_coords = _get_cellcoords(points, grid, node_cell_dicts, search_nneighbors, warn, strategy)
-    return PointEvalHandler(grid, cells, local_coords)
+    return PointEvalHandler(grid, cells, local_coords, false)
+end
+
+function PointEvalHandler(grid::AbstractGrid{dim}, points::AbstractVector{Vec{dim, T}}, subdomains::Vector{OrderedSet{Int}}; search_nneighbors = 3, warn::Bool = true, strategy = NewtonLineSearchPointFinder()) where {dim, T}
+    node_cell_dicts = _get_node_cell_map(grid, subdomains)
+    cells, local_coords = _get_cellcoords(points, grid, node_cell_dicts, search_nneighbors, warn, strategy)
+    return PointEvalHandler(grid, cells, local_coords, true)
 end
 
 function _get_cellcoords(points::AbstractVector{Vec{dim, T}}, grid::AbstractGrid, node_cell_dicts::Dict{C, Dict{Int, Vector{Int}}}, search_nneighbors, warn, strategy::NewtonLineSearchPointFinder) where {dim, T <: Real, C}
@@ -202,6 +209,31 @@
     return cell_dicts
 end
 
+function _get_node_cell_map(grid::AbstractGrid, subdomains::Vector{OrderedSet{Int}})
+    @inline function _check_contained_subdomain(i, subdomains)
+        for subdomain in subdomains
+            i ∈ subdomain && return true
+        end
+        return false
+    end
+    cells = getcells(grid)
+    C = eltype(cells) # possibly abstract
+    cell_dicts = Dict{Type{<:C}, Dict{Int, Vector{Int}}}()
+    ctypes = Set{Type{<:C}}(typeof(c) for c in cells)
+    for ctype in ctypes
+        cell_dict = cell_dicts[ctype] = Dict{Int, Vector{Int}}()
+        for (cellidx, cell) in enumerate(cells)
+            cell isa ctype || continue
+            _check_contained_subdomain(cellidx, subdomains) || continue
+            for node in cell.nodes
+                v = get!(Vector{Int}, cell_dict, node)
+                push!(v, cellidx)
+            end
+        end
+    end
+    return cell_dicts
+end
+
 """
     evaluate_at_points(ph::PointEvalHandler, dh::AbstractDofHandler, dof_values::Vector{T}, [fieldname::Symbol]) where T
     evaluate_at_points(ph::PointEvalHandler, proj::L2Projector, dof_values::Vector{T}) where T
@@ -226,6 +258,9 @@
         ph::PointEvalHandler{<:Any, T1}, dh::AbstractDofHandler, dof_vals::AbstractVector{T2},
         fname::Symbol = find_single_field(dh)
     ) where {T1, T2}
+    if !is_defined_on_full_domain(dh) && !ph.constructed_on_subdomain
+        @warn "DofHandler is defined on a true subdomain only, but the PointEvalHandler serached cells on the full grid. The evaluation assignment might be broken!" maxlog = 1
+    end
     npoints = length(ph.cells)
     # Figure out the value type by creating a dummy PointValues
     ip = getfieldinterpolation(dh, find_field(dh, fname))

test/test_pointevaluation.jl

@@ -1,4 +1,4 @@
-using Ferrite, Test
+using Ferrite, Test, Logging, OrderedCollections
 
 function test_pe_scalar_field()
     # isoparametric approximation
@@ -345,16 +345,15 @@ function test_pe_mixed_grid()
     # construct projector
     projector = L2Projector(ip_quad, mesh; set = getcellset(mesh, "quads"))
 
-    points = [Vec((x, 2x)) for x in range(0.0; stop = 1.0, length = 10)]
+    points = [Vec((x, 2x)) for x in range(0.0; stop = 0.5, length = 5)]
 
     # first alternative: L2Projection to dofs
     projector_values = project(projector, qp_vals_quads, qr)
     @test_logs min_level = Logging.Warn PointEvalHandler(mesh, points)
-    ph = PointEvalHandler(mesh, points)
+    ph = PointEvalHandler(mesh, points, [getcellset(mesh, "quads")])
     @test all(x -> x !== nothing, ph.cells)
     vals = evaluate_at_points(ph, projector, projector_values)
     @test vals[1:5] ≈ f.(points[1:5])
-    @test all(isnan, vals[6:end])
 
     # second alternative: assume a vector field :v
     dh = DofHandler(mesh)
@@ -427,6 +426,25 @@ function test_pe_first_point_missing()
     return
 end
 
+function test_pe_true_subdomain()
+    # Generate a subdofhandler on two elements
+    grid = generate_grid(Triangle, (4, 4), Vec((-2.0, -2.0)), Vec((2.0, 2.0)))
+    addcellset!(grid, "hole", x -> norm(x) ≤ 1.0)
+    subdomain = getcellset(grid, "hole")
+    dh = DofHandler(grid)
+    sdh = SubDofHandler(dh, subdomain)
+    add!(sdh, :v, Lagrange{RefTriangle, 1}())
+    close!(dh)
+
+    # Evaluate at all nodes of all triangles in "hole"
+    nodes = Vec{2, Float64}.([[0.0, -1.0], [0.0, 0.0], [-1.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
+    ph = PointEvalHandler(grid, nodes, [subdomain]; warn = true)
+    result = Ferrite.evaluate_at_points(ph, dh, zeros(ndofs(dh)))
+
+    @test !any(isnan.(result))
+    return
+end
+
 @testset "PointEvalHandler" begin
     @testset "scalar field" begin
         test_pe_scalar_field()
@@ -458,6 +476,10 @@ end
     @testset "failure cases" begin
         test_pe_first_point_missing()
     end
+
+    @testset "True Subdomain (issue #1181)" begin
+        test_pe_true_subdomain()
+    end
 end
 
 @testset "PointValues" begin