Fix Neumann BC assembly for TETRA4 and vector fields#243
Merged
Conversation
TETRA4 is a 4-node linear tetrahedron (p=1), not quadratic (p=2). The incorrect default would have triggered a p=2 Tet interpolation, assembling 10-DOF-per-element data structures for a 4-node element.
The connectivity array is packed with NNPE_vol entries per element. The old code used NNPE_surf as the stride, causing element e > 1 to read global node indices from the wrong positions. Additionally, it read NNPE_surf consecutive entries rather than mapping through the face's local node indices, so the returned global nodes were wrong for any face that doesn't use the first NNPE_surf nodes of the element. Fix: stride by NNPE_vol and select entries via boundary_dofs(ref_fe, side), which gives the 1-based local node indices within the volume element that belong to the requested face. Affects all element types (HEX8, TETRA4, etc.) wherever Neumann BCs are assembled.
Replaces two placeholder stubs (no assertions, commented out) with tests that have a known analytical solution: -∇²u = 0 on [0,1]² u = 0 at x=0 (Dirichlet, :left) ∂u/∂n = 1 at x=1 (Neumann, :right) ∂u/∂n = 0 at y=0,1 (natural) Exact solution: u(x,y) = x. Linear elements (QUAD4, TRI3) represent this exactly, so the tolerance is set to 1e-10. Exercises the full NBC assembly pipeline: surface_connectivity (stride fix), MappedH1OrL2SurfaceInterpolants (Jacobian fix), and the surface quadrature coordinate fix in RFE — none of which were covered by any existing test.
The indexin call in BCBookKeeping for sidesets had its arguments reversed: indexin(block_id_map, sideset_elements) returned positions of block elements within the sideset list, giving local indices [1,2] when the correct block-local element indices were [3,4]. Fix: swap to indexin(sideset_elements, block_id_map) to get the position of each sideset element within the block's element ID map. Without this fix, Neumann BCs were assembled at wrong DOFs. Also relax NBC test tolerance from 1e-10 to 1e-6 to accommodate floating-point rounding from mesh coordinate arithmetic on various backends (AMDGPU float32 etc.).
Replace scalar-only `Nvec * f_val` with `reduce(vcat, ntuple(i -> Nvec[i] * f_val, length(Nvec)))`,
which correctly computes the outer product for both scalar (ND=1) and vector (ND>1) fields.
Previously the assembly only worked for scalar fields; for 3D solid mechanics (ND=3) the
`Nvec * SVector{3}` multiplication was undefined.
Contributor
|
@lxmota I need to make a release of RFE so these test changes can pass. Please don't merge this right away. |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #243 +/- ##
=======================================
Coverage 62.43% 62.43%
=======================================
Files 60 60
Lines 3381 3381
=======================================
Hits 2111 2111
Misses 1270 1270 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
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
BCBookKeepingsideset element index mapping:indexinarguments were reversed, causing NBC to assemble at the wrong elements (e.g., elements 1–2 instead of 3–4 for the right boundary).Nvec * f_valmultiplication is only defined for scalar fields (ND=1). Replace withreduce(vcat, ntuple(i -> Nvec[i] * f_val, length(Nvec))), which correctly computes the outer product for both scalar (ND=1) and vector (ND>1) fields (e.g., 3D solid mechanics).surface_connectivityto use volume element stride: The surface connectivity kernel was using the surface NNPE as a stride, causing wrong node lookups for elements beyond the first.g = −flux).