Repository navigation
fix(assemble): select boundary facets by incidence and add facet normals - #69
Merged
Merged
Conversation
FacetAssembler picked every facet whose nodes all carry the boundary mask, so an interior facet with all nodes on the boundary (the diagonal of a corner cell, a chord of a curved boundary) was integrated too, once from each cell. On a 4x4 triangulated square the perimeter came out as 5.41 instead of 4; on a 3x3x3 Kuhn cube the surface as 15.8 instead of 6. A facet is now selected only if a single cell references it and all of its nodes are flagged (topology.boundary_facet_mask, which also matches facets shared by cells of different types). Facet integrands can take `n`, the outward unit normal at the facet quadrature point (Transformation.facet_normal, J^-T N normalised), in FacetAssembler and FacetBilinearAssembler; a point_data key named `n` keeps precedence. Fixed on the way, all exercised by the new tests: - the reference "outward" normals (element/normal.py) pointed inward; - Transformation.facet_quadrature was built in float32 and cast, which limited float64 facet integrals to ~1e-7; - the prism/pyramid path of FacetAssembler.__call__ failed (five-way unpack of a four-tuple, unflattened points, cell JxW instead of FxW, J J instead of J J^T in detF); - from_elements no longer overwrites the caller's connectivity dict. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Dst8zhPLnxVpZMBhR4Pvfs Signed-off-by: Shizheng Wen <wenshizheng0710@gmail.com>
Shizheng-Wen
added a commit
that referenced
this pull request
Oct 20, 2026
…#71) A thick-walled cylinder under internal pressure (Hill 1950), in plane strain, with J2 plasticity and linear isotropic hardening (perfect plasticity by default). Every load step is solved by Newton's method with the consistent tangent of the radial return, so the iterations converge quadratically: three to five per plastic load step. The plastic state lives at the quadrature points. The stress and the tangent moduli reach an internal-force NodeAssembler and a tangent ElementAssembler as quadrature_data (#70), and the pressure is a FacetAssembler load along the outward normal n (#69). Displacement control follows the collapse plateau; pressure control is an option. The quarter ring of six-node triangles is built with NumPy, so the example needs no gmsh. Tests check the return map and its tangent against forward-mode autodiff, the assembled stiffness against finite differences, the Lame elastic response, first yield at the inner wall, the limit pressure (2 / sqrt 3) sigma_y ln(b / a) of perfect plasticity (within 1 %), and CUDA against CPU. Docs: gallery page and solid-mechanics index entry, with zh_CN translations. Claude-Session: https://claude.ai/code/session_01Dst8zhPLnxVpZMBhR4Pvfs Signed-off-by: Shizheng Wen <wenshizheng0710@gmail.com> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
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
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.
What
FacetAssembler(andFacetBilinearAssembler) selected every facet whose nodes all carry the boundary mask. An interior facet whose nodes all lie on the boundary was therefore integrated too, once from each of its two cells: the diagonal of a corner cell, or a chord between two nodes of a curved boundary. Measured on structured meshes:A facet is now selected only if a single cell references it and all of its nodes are flagged. The new
topology.boundary_facet_maskdoes the incidence check. It also matches facets shared by cells of different types and handles prisms and pyramids. It only looks at the cells that touch the mask, so it stays cheap on large meshes.Facet integrands can also take
n, the outward unit normal at the facet quadrature point. It is available inFacetAssemblerandFacetBilinearAssembler, through the newTransformation.facet_normal(J⁻ᵀN, normalised), and follows curved P2 facets point by point. Pressure loads such as-p nneed it. Apoint_datakey namednkeeps precedence, so existing forms behave as before.Fixed on the way
The new tests exercise each of these:
element/normal.pyreturned inward normals as "outward". Onlynanson_scaleused them, and it only takes their norm. The TODOs intests/element/test_normal.pynow check the direction.Transformation.facet_quadraturewas built in float32 and then cast, which limited float64 facet integrals to about 1e-7 relative error.FacetAssembler.__call__could not run: it unpacked a four-tuple into five names, passed unflattened points toeval_shape_val, used the cellJxWinstead ofFxW, and computed J·J instead of J·Jᵀ indetF.from_elementsreplaced the caller's connectivity dict entries with the boundary cells.Tests
tests/assemble/test_facet_selection.py(new) covers:∫ n ds = 0and∫ x·n ds = D|Ω|to round-off, also on curved P2 cells.update_points.ninFacetBilinearAssembler, andpoint_dataprecedence.Full suite on Clariden (aarch64) in a CPU conda env that has gmsh:
The new tests also pass on a GH200 node, CUDA test included. The docs build (EN + zh) adds no warnings (21 → 21).
nanson_scaleitself is unchanged. It still lacksdim=in its norm, still mixes dtypes, and its weights miss the reference-facet area factor. Nothing calls it.🤖 Generated with Claude Code
https://claude.ai/code/session_01Dst8zhPLnxVpZMBhR4Pvfs