Fix PointEvalHandler on true subdomains #1182

termi-official · 2025-04-01T19:22:48Z

Fixes #1181 .

codecov · 2025-04-01T19:28:11Z

Codecov Report

Attention: Patch coverage is 97.14286% with 1 line in your changes missing coverage. Please review.

Project coverage is 93.76%. Comparing base (87b0888) to head (f4ea8b1).

Files with missing lines	Patch %	Lines
src/PointEvalHandler.jl	96.55%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1182      +/-   ##
==========================================
+ Coverage   93.75%   93.76%   +0.01%     
==========================================
  Files          39       39              
  Lines        6241     6275      +34     
==========================================
+ Hits         5851     5884      +33     
- Misses        390      391       +1

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

kimauth

Thanks! I left a few questions, perhaps we can get a slightly nicer user interface 🙂

kimauth · 2025-04-03T14:10:38Z

src/Dofs/DofHandler.jl

+    end
+    return getncells(get_grid(dh)) == num_cells
+end
+


Wouldn't you have the same problem with a field that's defined on a subdomain, even if all cells are part of some SubDofHandler? In that case this check isn't sufficient, is it?

Good catch, you are right. Let me fix that.

kimauth · 2025-04-03T14:53:31Z

src/PointEvalHandler.jl

    return cell_dicts
 end

+function _get_node_cell_map(grid::AbstractGrid, subdomains::Vector{OrderedSet{Int}})


Do we need this to be its own function or could we just modify the original _get_node_cell_map to accept an optional keyword argument cellset that defaults to all cells? And why does subdomains need to be a vector here?

Do we need this to be its own function or could we just modify the original _get_node_cell_map to accept an optional keyword argument cellset that defaults to all cells?

I simply prefer to use multiple dispatch over kwargs/if-else constructs using isa, because it is easier to extend later when we add new features and feels cleaner.

why does subdomains need to be a vector here?

The field to which the point is associated might span multiple subdomains (e.g. we have the cellsets named "A", "B" and "C", then we might want to search only on "B" and "C"). Think e.g. about discretizing the torso of some kind of animal, then each organ is potentially at least one subdomain. Does that make sense?

Is the performance advantage sufficient to motivate the difference to other interfaces? E.g. that we don't allow users to pass a vector of facetsets to Dirichlet, but require their union. IMO the same is reasonable here, but we should allow an Ferrite.AbstractVecOrSet{<:Integer}, allow defaulting to 1:getncells(grid). Then it would be the same function in both cases. I don't think a kwarg is required here since this is internal.

Is the performance advantage sufficient to motivate the difference to other interfaces?

Sorry, but that does not make sense for me. First, I do not see what this has to do with performance advantages. Second, we already use Vector internally in PointEvalHandler (see

Ferrite.jl/src/PointEvalHandler.jl

Line 238 in 6a14d7a

evaluate_at_points!(out_vals, ph, dh, dof_vals, fname, func_interpolations)

) and I just sticked to that design. Third, it does not make sense to have users providing the exact split subdomains in one part of the API (DofHandler), just merge these subdomains into a single set again to call into the PointEvalHandler.

All other handlers provide add! and close! and PointEvalHandler does not on the user-facing side. I might need to change the design of the PointEvalHandler interface to match the other interfaces in an attempt to resolve another bug Kim found https://github.com/Ferrite-FEM/Ferrite.jl/pull/1182/files/f4ea8b1a2aa2dc61d1beed2c0e60f91e38f07222#r2027088188 in this branch.

We don't allow users to pass a vector of facetsets to Dirichlet, but require their union.

I am not sure if that is true. The following code works without problem, right?

ch = ... add!(ch, Dirichlet(:v, ∂Ω1, x -> sin(x[1])) add!(ch, Dirichlet(:v, ∂Ω2, x -> sin(x[1])) close!(ch)

Or is there something I am missing that requires us to use:

ch = ... add!(ch, Dirichlet(:v, union(∂Ω1, ∂Ω2) x -> sin(x[1])) close!(ch)

? If that is a problem, can we fix it without merging these sets?

I would even argue for simplicity here to just have some second constructor here PointEvalHandler(dh::DofHandler, points::AbstractVector{Vec{dim,T}}, field_name::Symbol; ... ) which takes care of extracting the right cellsets from the DofHandler.

I am not sure if that is true. The following code works without problem, right?

ch = ... add!(ch, Dirichlet(:v, ∂Ω1, x -> sin(x[1])) add!(ch, Dirichlet(:v, ∂Ω2, x -> sin(x[1])) close!(ch)

Sorry, maybe I misunderstood, but I thought you are proposing supporting the equivalent of

add!(ch, Dirichlet(:v, [∂Ω1, ∂Ω2], x -> sin(x[1]))

First, I do not see what this has to do with performance advantages.

The only one I could see was avoiding creating a new set, which should be rather cheap in most cases.

Second, we already use Vector internally in PointEvalHandler

I don't follow this - there is no cellset/subdomain/subdomains internally now, right?
I meant that instead of restricting the domains to be described by OrderedSet, we could also allow e.g. a range - but you are right that a Vector would be slow, so perhaps cellset::Union{AbstractSet{<:Integer}, UnitRange{<:Integer}} is better!

I might need to change the design of the PointEvalHandler interface to match the other interfaces

That sounds like a good idea to me, interested to see how that would look! I just didn't think the passing of a Vector of sets was aligned with such an interface, the ability to add! would definitively be more Ferrite-like!

I also like the short-hand version with the fieldname :)

I don't follow this - there is no cellset/subdomain/subdomains internally now, right?

Not directly, which is part of the problem. See

Ferrite.jl/src/PointEvalHandler.jl

Lines 237 to 238 in 6a14d7a

func_interpolations = get_func_interpolations(dh, fname)

evaluate_at_points!(out_vals, ph, dh, dof_vals, fname, func_interpolations)

I will try to bring this PR over the finish line early next week.

kimauth · 2025-04-03T14:54:41Z

src/PointEvalHandler.jl

+    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}


Would be nice if we could have a keword argument cellset here instead of subdomains, see my comment below.

I agree that a kwarg makes the user-experience better and we don't lock in on arg ordering.

Even if an arg defaulting to 1:getncells(grid) is alright from the convention that args say what and kwargs say how.

KnutAM

Nice to get this fixed 🚀

KnutAM · 2025-04-11T12:00:15Z

src/PointEvalHandler.jl

+        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


Isn't this enough, or was a performance optimization required?

Suggested change

_check_contained_subdomain(cellidx, subdomains) || continue

any(set -> cellidx in set, subdomains) || continue

But if following other comments, would of course also remove the any...

KnutAM · 2025-04-11T12:08:14Z

src/PointEvalHandler.jl

+    @inline function _check_contained_subdomain(i, subdomains)
+        for subdomain in subdomains
+            i ∈ subdomain && return true
+        end
+        return false
+    end


Suggested change

@inline function _check_contained_subdomain(i, subdomains)

for subdomain in subdomains

i ∈ subdomain && return true

end

return false

end

KnutAM · 2025-04-11T12:20:11Z

src/PointEvalHandler.jl

    return cell_dicts
 end

+function _get_node_cell_map(grid::AbstractGrid, subdomains::Vector{OrderedSet{Int}})


Is the performance advantage sufficient to motivate the difference to other interfaces? E.g. that we don't allow users to pass a vector of facetsets to Dirichlet, but require their union. IMO the same is reasonable here, but we should allow an Ferrite.AbstractVecOrSet{<:Integer}, allow defaulting to 1:getncells(grid). Then it would be the same function in both cases. I don't think a kwarg is required here since this is internal.

KnutAM · 2025-04-11T12:24:51Z

src/PointEvalHandler.jl

        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


Suggested change

@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

@warn "DofHandler is defined on a true subdomain only, but the PointEvalHandler searched cells on the full grid. The evaluation assignment might be broken!" maxlog = 1

KnutAM · 2025-04-11T12:27:11Z

src/PointEvalHandler.jl

+    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}


I agree that a kwarg makes the user-experience better and we don't lock in on arg ordering.

Even if an arg defaulting to 1:getncells(grid) is alright from the convention that args say what and kwargs say how.

KnutAM · 2025-04-11T12:27:55Z

test/test_pointevaluation.jl

@@ -1,4 +1,4 @@
-using Ferrite, Test
+using Ferrite, Test, Logging, OrderedCollections


Why is this even here, I don't think we do this consistently?

I am not sure why these are missing from the test. Running only parts of the tests is right now quite painful, as we include some of the deps in the runner and sometime some of the deps in the file containing the testsets. Would you agree that it would be better if the different files containing the test sections would be self-contained, so we can run them directly for testing instead of modifying the test codes every time we only want to test part of the API?

I totally agree that this can be very useful, I was just not aware that we have this as a convention to do!

Would of course then be nicer if it would error if one doesn't do that, otherwise it is quite hard to make sure it is done consistently. This also includes e.g. including test_utils.jl. Perhaps it could rather make sense to have a file initialize_tests.jl that is included at the beginning of runtests.jl, and running individual tests can then always be done by including that file before the specific file (or snippets of code)? Of course a separate discussion...

Yea we should probably follow something like this here https://github.com/SciML/OrdinaryDiffEq.jl/blob/1804a1bbc436c88580310316e69ff5a091f65eb1/test/runtests.jl , especially if we start expanding convergence rate tests and GPU test coverage.

termi-official added 3 commits April 1, 2025 21:19

Fix PointEvalHandler on true subdomains

e61b317

Dirt.

8a7cb0c

Make Runic happy

f4ea8b1

termi-official added the awaiting review PR is finished from the authors POV, waiting for feedback label Apr 1, 2025

termi-official requested a review from kimauth April 1, 2025 20:16

termi-official mentioned this pull request Apr 1, 2025

FIMH2021 verification JuliaHealth/Thunderbolt.jl#209

Draft

kimauth reviewed Apr 3, 2025

View reviewed changes

termi-official removed the awaiting review PR is finished from the authors POV, waiting for feedback label Apr 3, 2025

KnutAM reviewed Apr 11, 2025

View reviewed changes

termi-official marked this pull request as draft May 6, 2025 10:57

	func_interpolations = get_func_interpolations(dh, fname)
	evaluate_at_points!(out_vals, ph, dh, dof_vals, fname, func_interpolations)

	_check_contained_subdomain(cellidx, subdomains) \|\| continue
	any(set -> cellidx in set, subdomains) \|\| continue

	@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
	@warn "DofHandler is defined on a true subdomain only, but the PointEvalHandler searched cells on the full grid. The evaluation assignment might be broken!" maxlog = 1

		@@ -1,4 +1,4 @@
		using Ferrite, Test
		using Ferrite, Test, Logging, OrderedCollections

Fix PointEvalHandler on true subdomains #1182

Are you sure you want to change the base?

Fix PointEvalHandler on true subdomains #1182

Uh oh!

Conversation

termi-official commented Apr 1, 2025

Uh oh!

codecov bot commented Apr 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

kimauth left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KnutAM Apr 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KnutAM Apr 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KnutAM left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KnutAM Apr 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov bot commented Apr 1, 2025 •

edited

Loading

KnutAM Apr 11, 2025 •

edited

Loading

KnutAM Apr 11, 2025 •

edited

Loading

KnutAM Apr 11, 2025 •

edited

Loading