Fix MappedH1OrL2SurfaceInterpolants Jacobian for surface quadrature#98
Merged
Conversation
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.
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.
…rmalization 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).
cmhamel
approved these changes
Mar 16, 2026
Contributor
Author
|
@JuliaRegistrator register |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
det_J = norm(X_diff)(physical edge length) for all 1D boundaries, which is only correct for shifted elements parameterized on [0,1]. For unshifted edges on [−1,1], the reference length is 2, sodet_J = physical_length / 2. This caused a 2× force error for QUAD4 Neumann BCs.det_Jby the appropriate reference area.SVector{2}hardcoded dimension: The physical quadrature point coordinateX_qwas constructed asSVector{2, ...}, failing for 3D elements. Replace withX_face * N_reducedwhich infers the correct dimension.