Fix Neumann BC assembly for TETRA4 and vector fields

src/FunctionSpaces.jl

@@ -9,7 +9,7 @@ const _default_p = Dict{String, Int}(
   "TRI6"    => 2,
   "TET"     => 1,
   "TETRA"   => 1,
-  "TETRA4"  => 2,
+  "TETRA4"  => 1,
   "TETRA10" => 2
 )
 const _default_q = Dict{String, Int}(

src/assemblers/NeumannBC.jl

@@ -52,13 +52,10 @@ function _assemble_block_vector_neumann_bc!(
       Nvec = interps.N_reduced
       JxW = interps.JxW
 
-      # TODO Clean this up
       f_val = vals[q, e]
-      if length(f_val) == 1
-        f_val = f_val[1]
-      end
-
-      R_el = R_el + JxW * Nvec * f_val
+      # Outer product: R_el[ND*(i-1)+α] += JxW * N_i * f_α
+      # Works for both scalar (ND=1) and vector (ND>1) fields.
+      R_el = R_el + JxW * reduce(vcat, ntuple(i -> Nvec[i] * f_val, length(Nvec)))
     end
     surf_conns = surface_connectivity(ref_fe, conns, side, e, 1)
     @views _assemble_element!(field, R_el, surf_conns, e, 0, 0)
@@ -88,13 +85,8 @@ KA.@kernel function _assemble_block_vector_neumann_bc_kernel!(
     Nvec = interps.N_reduced
     JxW = interps.JxW
 
-    # TODO Clean this up
     f_val = vals[q, E]
-    if length(f_val) == 1
-      f_val = f_val[1]
-    end
-
-    R_el = R_el + JxW * Nvec * f_val
+    R_el = R_el + JxW * reduce(vcat, ntuple(i -> Nvec[i] * f_val, length(Nvec)))
   end
 
   surf_conns = surface_connectivity(ref_fe, conns, side, E, 1)

src/bcs/BoundaryConditions.jl

@@ -129,7 +129,7 @@ function BCBookKeeping(
 
     @assert length(unique(blocks)) == 1 "Sidesets need to be in a single block"
     block_name = mesh.element_block_names[blocks[1]]
-    indices_in_sset = indexin(mesh.element_id_maps[block_name], elements)
+    indices_in_sset = indexin(elements, mesh.element_id_maps[block_name])
     filter!(x -> x !== nothing, indices_in_sset)
     elements = convert(Vector{Int}, indices_in_sset)
     # display(elements)

src/fields/Connectivity.jl

@@ -66,12 +66,15 @@ end
 
 # GPU safe
 @inline function surface_connectivity(ref_fe::ReferenceFE, conn_data, side::Int, e::Int, boffset::Int)
-    NNPE = ReferenceFiniteElements.num_cell_dofs(
-        ReferenceFiniteElements.boundary_element(ref_fe.element, side)
-    )
-    base = boffset + (e - 1) * NNPE
-    data = ntuple(i -> conn_data[base + i - 1], NNPE)
-    return SVector{NNPE, Int}(data)
+    # Stride through conn_data using the VOLUME element DOF count, not the surface element's.
+    # The connectivity array is packed with NNPE_vol entries per element; using NNPE_surf as
+    # the stride reads from the wrong position for element e > 1.
+    NNPE_vol  = ReferenceFiniteElements.num_cell_dofs(ref_fe)
+    face_nodes = ReferenceFiniteElements.boundary_dofs(ref_fe, side)  # 1-based local indices
+    NNPE_surf = length(face_nodes)
+    base = boffset + (e - 1) * NNPE_vol
+    data = ntuple(i -> conn_data[base + face_nodes[i] - 1], NNPE_surf)
+    return SVector{NNPE_surf, Int}(data)
 end
 
 function unsafe_connectivity(conn::Connectivity, e::Int, b::Int)

test/poisson/TestPoisson.jl

@@ -32,8 +32,8 @@ function test_poisson(backend, cond, nlsolver, lsolver)
   test_poisson_dirichlet_structured_mesh_quad4(backend, cond, nlsolver, lsolver)
   test_poisson_dirichlet_structured_mesh_tri3(backend, cond, nlsolver, lsolver)
   test_poisson_neumann(backend, cond, nlsolver, lsolver)
-  # test_poisson_neumann_structured_mesh_quad4(backend, cond, nlsolver, lsolver)
-  # test_poisson_neumann_structured_mesh_tri3(backend, cond, nlsolver, lsolver)
+  test_poisson_neumann_structured_mesh_quad4(backend, cond, nlsolver, lsolver)
+  test_poisson_neumann_structured_mesh_tri3(backend, cond, nlsolver, lsolver)
 end
 
 function test_poisson_dirichlet(
@@ -489,38 +489,45 @@ function test_poisson_neumann_structured_mesh_quad4(
   dev, use_condensed,
   nsolver, lsolver
 )
+  # Laplace equation (-∇²u = 0) on [0,1]² with mixed BCs:
+  #   u = 0            at x=0  (Dirichlet, :left)
+  #   ∂u/∂n = 1        at x=1  (Neumann,   :right;  outward normal = +x̂)
+  #   ∂u/∂n = 0        at y=0,1 (natural — not prescribed)
+  #
+  # Exact solution: u(x,y) = x
+  # QUAD4 bilinear elements represent linear functions exactly, so FEM error = 0.
+  #
+  # FEC sign convention: the NBC assembler adds +∫g·N dΓ to the residual, so
+  #   g = −(∂u/∂n).  For ∂u/∂n = +1 we must pass g = −1.
+  f_zero    = (_, _) -> 0.0
+  nbc_minus = (_, _) -> SVector{1, Float64}(-1.)
+
   mesh = StructuredMesh("quad", (0., 0.), (1., 1.), (11, 11))
-  V = FunctionSpace(mesh, H1Field, Lagrange) 
-  physics = Poisson(f)
+  V = FunctionSpace(mesh, H1Field, Lagrange)
+  physics = Poisson(f_zero)
   props = create_properties(physics)
   u = ScalarFunction(V, :u)
   asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
 
-  # setup and update bcs
   dbcs = DirichletBC[
-    DirichletBC(:u, bc_func; sideset_name = :bottom),
-    DirichletBC(:u, bc_func; sideset_name = :right)
+    DirichletBC(:u, bc_func; sideset_name = :left),
   ]
 
   nbcs = NeumannBC[
-    NeumannBC(:u, bc_func_neumann, :top),
-    NeumannBC(:u, bc_func_neumann, :left)
+    NeumannBC(:u, nbc_minus, :right),
   ]
 
-  # direct solver test
-  # setup the parameters
   p = create_parameters(
-    mesh, asm, physics, props; 
+    mesh, asm, physics, props;
     dirichlet_bcs=dbcs,
     neumann_bcs=nbcs
   )
 
   if dev != cpu
     p = p |> dev
-    asm = asm |> dev 
+    asm = asm |> dev
   end
 
-  # setup solver and integrator
   solver = nsolver(lsolver(asm))
   integrator = QuasiStaticIntegrator(solver)
   evolve!(integrator, p)
@@ -529,57 +536,49 @@ function test_poisson_neumann_structured_mesh_quad4(
     p = p |> cpu
   end
 
-  # TODO make a neumann gold file
-  # U = p.h1_field
-
-  # pp = PostProcessor(mesh, output_file, u)
-  # write_times(pp, 1, 0.0)
-  # write_field(pp, 1, ("u",), U)
-  # close(pp)
-
-  # if !Sys.iswindows()
-  #   @test exodiff(output_file, gold_file)
-  # end
-  # rm(output_file; force=true)
-  display(solver.timer)
+  # u(x,y) = x  →  max = 1.0  (x=1),  min = 0.0  (x=0, Dirichlet)
+  @test maximum(p.h1_field.data) ≈ 1.0 atol=1e-6
+  @test minimum(p.h1_field.data) ≈ 0.0 atol=1e-6
 end
 
 function test_poisson_neumann_structured_mesh_tri3(
   dev, use_condensed,
   nsolver, lsolver
 )
+  # Same problem as the QUAD4 variant above, meshed with TRI3 elements.
+  # TRI3 linear triangles represent linear functions exactly, so FEM error = 0.
+  #
+  # FEC sign convention: the NBC assembler adds +∫g·N dΓ to the residual, so
+  #   g = −(∂u/∂n).  For ∂u/∂n = +1 we must pass g = −1.
+  f_zero    = (_, _) -> 0.0
+  nbc_minus = (_, _) -> SVector{1, Float64}(-1.)
+
   mesh = StructuredMesh("tri", (0., 0.), (1., 1.), (11, 11))
-  V = FunctionSpace(mesh, H1Field, Lagrange) 
-  physics = Poisson(f)
+  V = FunctionSpace(mesh, H1Field, Lagrange)
+  physics = Poisson(f_zero)
   props = create_properties(physics)
   u = ScalarFunction(V, :u)
   asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
 
-  # setup and update bcs
   dbcs = DirichletBC[
-    DirichletBC(:u, bc_func; sideset_name = :bottom),
-    DirichletBC(:u, bc_func; sideset_name = :right)
+    DirichletBC(:u, bc_func; sideset_name = :left),
   ]
 
   nbcs = NeumannBC[
-    NeumannBC(:u, bc_func_neumann, :top),
-    NeumannBC(:u, bc_func_neumann, :left)
+    NeumannBC(:u, nbc_minus, :right),
   ]
 
-  # direct solver test
-  # setup the parameters
   p = create_parameters(
-    mesh, asm, physics, props; 
+    mesh, asm, physics, props;
     dirichlet_bcs=dbcs,
     neumann_bcs=nbcs
   )
 
   if dev != cpu
     p = p |> dev
-    asm = asm |> dev 
+    asm = asm |> dev
   end
 
-  # setup solver and integrator
   solver = nsolver(lsolver(asm))
   integrator = QuasiStaticIntegrator(solver)
   evolve!(integrator, p)
@@ -588,17 +587,6 @@ function test_poisson_neumann_structured_mesh_tri3(
     p = p |> cpu
   end
 
-  # TODO make a neumann gold file
-  # U = p.h1_field
-
-  # pp = PostProcessor(mesh, output_file, u)
-  # write_times(pp, 1, 0.0)
-  # write_field(pp, 1, ("u",), U)
-  # close(pp)
-
-  # if !Sys.iswindows()
-  #   @test exodiff(output_file, gold_file)
-  # end
-  # rm(output_file; force=true)
-  display(solver.timer)
+  @test maximum(p.h1_field.data) ≈ 1.0 atol=1e-6
+  @test minimum(p.h1_field.data) ≈ 0.0 atol=1e-6
 end