diff --git a/NEWS.md b/NEWS.md index 4a9476cfd..618bb6ab9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- Make `compute_facet_owners(...)` more flexible, allowing the user to provide a function to select the owner from neighboring cells. Since PR[#1291](https://github.com/gridap/Gridap.jl/pull/1291). + ## [0.20.5] - 2026-04-28 ### Fixed diff --git a/src/FESpaces/Pullbacks.jl b/src/FESpaces/Pullbacks.jl index 7c565bf7c..357129804 100644 --- a/src/FESpaces/Pullbacks.jl +++ b/src/FESpaces/Pullbacks.jl @@ -188,7 +188,7 @@ function evaluate!(cache,k::NormalSignMap,reffe,facet_own_dofs,cell) return Diagonal(dof_sign) end -function compute_facet_owners(model::DiscreteModel{Dc}) where {Dc} +function compute_facet_owners(model::DiscreteModel{Dc}, select_nbor=maximum) where {Dc} topo = get_grid_topology(model) facet_to_cell = get_faces(topo, Dc-1, Dc) @@ -196,13 +196,17 @@ function compute_facet_owners(model::DiscreteModel{Dc}) where {Dc} owners = Vector{Int32}(undef, nfacets) for facet in 1:nfacets facet_cells = view(facet_to_cell, facet) - owners[facet] = first(facet_cells) + @check !isempty(facet_cells) "Facet $facet has no adjacent cells" + selected_owner = select_nbor(facet_cells) + @check selected_owner isa Integer "select_nbor must return an integer owner for facet $facet, got $(typeof(selected_owner))" + owner = Int(selected_owner) + @check owner != 0 "select_nbor returned invalid owner 0 for facet $facet; expected one of $(collect(facet_cells))" + @check owner in facet_cells "select_nbor returned invalid owner $owner for facet $facet; expected one of $(collect(facet_cells))" + owners[facet] = Int32(owner) end - return owners end - ################# # DOFScalingMap # ################# diff --git a/test/FESpacesTests/DivConformingFESpacesTests.jl b/test/FESpacesTests/DivConformingFESpacesTests.jl index 2d5c72e99..d6706f6f0 100644 --- a/test/FESpacesTests/DivConformingFESpacesTests.jl +++ b/test/FESpacesTests/DivConformingFESpacesTests.jl @@ -294,4 +294,41 @@ test_div_v_q_equiv(U,V,P,Q,Ω) end +# This test checks that the facet owner cell criterion +# underlying the normal sign map is consistent no matter +# how cells sharing a facet are listed in the model. +@testset "NormalSignMap" begin + + # Create domain + domain = (0,1,0,1) + cells = (2,1) + model = CartesianDiscreteModel(domain,cells) + order = 0 + reffe = ReferenceFE(raviart_thomas,Float64,order) + V1 = TestFESpace(model,reffe) + uh1 = FEFunction(V1,ones(num_free_dofs(V1))) + + topo = get_grid_topology(model) + facet_cells = get_faces(topo,1,2) + interior_facets = 0 + for facet in 1:(length(facet_cells.ptrs)-1) + first_cell = facet_cells.ptrs[facet] + last_cell = facet_cells.ptrs[facet+1] - 1 + if last_cell - first_cell + 1 == 2 + facet_cells.data[first_cell], facet_cells.data[last_cell] = + facet_cells.data[last_cell], facet_cells.data[first_cell] + interior_facets += 1 + end + end + @test interior_facets == 1 + V2 = TestFESpace(model,reffe) + uh2 = FEFunction(V2,get_free_dof_values(uh1)) + + data_uh1 = get_data(uh1) + data_uh2 = get_data(uh2) + + @test evaluate(data_uh1[1],[Point(0.5,0.5)])[1] ≈ evaluate(data_uh2[1],[Point(0.5,0.5)])[1] + @test evaluate(data_uh1[2],[Point(0.5,0.5)])[1] ≈ evaluate(data_uh2[2],[Point(0.5,0.5)])[1] +end + end # module