From 4951699a152b42c6226504fea1849f7c818967c3 Mon Sep 17 00:00:00 2001 From: Alejandro Mota Date: Sun, 15 Mar 2026 22:35:33 -0700 Subject: [PATCH 1/3] Fix face 4 surface quadrature coordinates for Tet elements MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The last block in surface_quadrature_points_and_weights was writing face 4's (z=0) quadrature coordinates to ξ_return[:,:,1] instead of ξ_return[:,:,4], corrupting face 1 and leaving face 4 as all-zeros. This caused incorrect Neumann BC assembly on Tet faces 1 and 4. --- src/elements/Tet.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/elements/Tet.jl b/src/elements/Tet.jl index 1c2b05b..4303f11 100644 --- a/src/elements/Tet.jl +++ b/src/elements/Tet.jl @@ -222,9 +222,9 @@ function surface_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLoba ξ_return[2, :, 3] .= ξs[1, :] ξ_return[3, :, 3] .= ξs[2, :] # - ξ_return[1, :, 1] .= ξs[1, :] - ξ_return[2, :, 1] .= ξs[2, :] - ξ_return[3, :, 1] .= 0. + ξ_return[1, :, 4] .= ξs[1, :] + ξ_return[2, :, 4] .= ξs[2, :] + ξ_return[3, :, 4] .= 0. for n in 1:4 w_return[:, n] .= ws From ccf063e71567412078da9257e7ac975d22c7ae68 Mon Sep 17 00:00:00 2001 From: Alejandro Mota Date: Mon, 16 Mar 2026 10:11:10 -0700 Subject: [PATCH 2/3] Fix MappedH1OrL2SurfaceInterpolants Jacobian for 3D elements MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The old implementation computed det_J as the length of the edge between nodes 1 and 2 of the volume element (a 1D Jacobian), and hardcoded X_q as SVector{2}, both wrong for 3D volume elements whose faces are 2D (Tri for Tet, Quad for Hex). For 2D boundary elements (Tri/Quad), compute the area Jacobian via the cross product of the two face tangent vectors: t1 = X_face[:,2] - X_face[:,1] t2 = X_face[:,3] - X_face[:,1] det_J = |t1 × t2|, n = (t1 × t2) / det_J For 1D boundary elements (Edge), retain the existing edge-length Jacobian with the hard-coded reference normal. Also fixes: X_q now inferred from X_face * N_reduced (correct spatial dimension) instead of SVector{2,...} (hardcoded 2D); X_face is built from the face's local node indices via boundary_dofs, not from the first two columns of the volume element coordinate matrix. --- src/ReferenceFE.jl | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/ReferenceFE.jl b/src/ReferenceFE.jl index 8998048..c02eab8 100644 --- a/src/ReferenceFE.jl +++ b/src/ReferenceFE.jl @@ -477,21 +477,29 @@ function MappedH1OrL2SurfaceInterpolants(e::ReferenceFE, X, q::Integer, f::Integ NNPS = num_cell_dofs(boundary_element(e, f)) edge_nodes = SVector{NNPS, Int}(boundary_dofs(e, f)) N_reduced = SVector{NNPS, eltype(N)}(@views N[edge_nodes]) - n = boundary_normal(e, f) - # jacobian - X_diff = X[:, 2] - X[:, 1] - det_J = norm(X_diff) - # interpolate coordinates - edge_nodes = boundary_dofs(e, f) - X_q = SVector{2, eltype(X)}(@views X[:, edge_nodes] * N_reduced) - - # JxW + # Physical coordinates of the face nodes (columns of X indexed by local face-node ids) + X_face = X[:, edge_nodes] + + # Jacobian and outward unit normal, dispatched on boundary element dimension: + # 1D boundary (Edge): line Jacobian = edge length + # 2D boundary (Tri/Quad): area Jacobian = |t1 × t2|, normal computed from cross product + if dimension(boundary_element(e, f)) == 1 + X_diff = X_face[:, 2] - X_face[:, 1] + det_J = norm(X_diff) + n = boundary_normal(e, f) + else + t1 = X_face[:, 2] - X_face[:, 1] + t2 = X_face[:, 3] - X_face[:, 1] + n_raw = cross(t1, t2) + det_J = norm(n_raw) + n = n_raw / det_J + end + + # Physical coordinates at the quadrature point + X_q = X_face * N_reduced + JxW = det_J * w - - # TODO below incorrect. Not giving correct gradient - # or normal - # @show N return MappedH1OrL2SurfaceInterpolants(X_q, N, N_reduced, JxW, n) end From 1c60c2de91545503a4adbd4c1abcc2c9dde85076 Mon Sep 17 00:00:00 2001 From: Alejandro Mota Date: Mon, 16 Mar 2026 11:41:17 -0700 Subject: [PATCH 3/3] Fix MappedH1OrL2SurfaceInterpolants Jacobian for reference element normalization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The physical Jacobian must account for the reference element's measure: - 1D boundary (unshifted Edge on [-1,1]): det_J = physical_length / 2 - 1D boundary (shifted Edge on [0,1]): det_J = physical_length (unchanged) - 2D boundary (shifted Tri on [0,1]²): det_J = |t1×t2| (unchanged) - 2D boundary (unshifted Quad on [-1,1]²): det_J = |t1×t2| / 4 Without this fix, QUAD4 NBC forces were 2x too large (reference interval length 2 vs. physical length h), while TRI3 was already correct because the shifted edge uses weights that sum to 1 (not 2). --- src/ReferenceFE.jl | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/ReferenceFE.jl b/src/ReferenceFE.jl index c02eab8..207438f 100644 --- a/src/ReferenceFE.jl +++ b/src/ReferenceFE.jl @@ -481,19 +481,30 @@ function MappedH1OrL2SurfaceInterpolants(e::ReferenceFE, X, q::Integer, f::Integ # Physical coordinates of the face nodes (columns of X indexed by local face-node ids) X_face = X[:, edge_nodes] - # Jacobian and outward unit normal, dispatched on boundary element dimension: - # 1D boundary (Edge): line Jacobian = edge length - # 2D boundary (Tri/Quad): area Jacobian = |t1 × t2|, normal computed from cross product - if dimension(boundary_element(e, f)) == 1 + # Jacobian and outward unit normal, dispatched on boundary element dimension. + # The quadrature weights w come from the reference element (e.g. Gauss-Legendre + # on [-1,1] for unshifted edges, or on [0,1] for shifted edges). The physical + # Jacobian must map from reference measure to physical measure, so + # det_J = physical_length_or_area / reference_length_or_area. + # For shifted edges ([0,1], ref_length=1) det_J = physical_length. + # For unshifted edges ([-1,1], ref_length=2) det_J = physical_length / 2. + # Analogously for 2D faces (shifted Tri on [0,1]², area=1/2; unshifted Quad + # on [-1,1]², area=4). + be = boundary_element(e.element, f) + if dimension(be) == 1 X_diff = X_face[:, 2] - X_face[:, 1] - det_J = norm(X_diff) + ref_len = be.shifted ? 1.0 : 2.0 + det_J = norm(X_diff) / ref_len n = boundary_normal(e, f) else t1 = X_face[:, 2] - X_face[:, 1] t2 = X_face[:, 3] - X_face[:, 1] n_raw = cross(t1, t2) - det_J = norm(n_raw) - n = n_raw / det_J + # reference area: shifted Tri (right triangle on [0,1]²) has area 1/2; + # unshifted Quad on [-1,1]² has area 4. + ref_area = isa(be, AbstractTri) ? 0.5 : 4.0 + det_J = norm(n_raw) / ref_area + n = n_raw / norm(n_raw) end # Physical coordinates at the quadrature point