Strengthen/bug-fix compute_facet_owners

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

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Alberto F. Martin <alberto.f.martin@anu.edu.au>"]
-version = "0.20.5"
+version = "0.20.6"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/FESpaces/Pullbacks.jl

@@ -188,21 +188,20 @@ 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)
 
   nfacets = num_faces(topo, Dc-1)
   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"
+    owners[facet] = select_nbor(facet_cells)
   end
-
   return owners
 end
 
-
 #################
 # DOFScalingMap #
 #################

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