From 85900fe680b0d67abad3fddfcf38544f154a8556 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 11 May 2026 02:26:23 +0300 Subject: [PATCH] feat(mesh): add structured Quad4/Seg2/Cook meshes and face-aware extract_surface MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Structured mesh helpers now cover 2D quads in the xy plane, a 1D Seg2 chain embedded in R^3, and Cook’s membrane on a bilinear Quad4 map, so plane tests and benchmarks can build reference geometry without ad hoc connectivity. extract_surface gains an explicit volume-local face index and uses faces(T) when the face vertex count matches the surface topology, which avoids picking the wrong face by default for low-order hex meshes that list true face corners. - structured.jl: SPDX header; add create_structured_line_mesh(Seg2; x0,x1,nx,y,z) - structured.jl: add create_structured_box_mesh(Quad4; xmin,xmax,…,z) with :xmin/:xmax/:ymin/:ymax node sets - structured.jl: add create_cook_membrane_mesh(Quad4, nx, ny; scale) with classical Cook corner map and same node-set convention - mesh.jl: extract_surface(mesh, set; local_face=1) selects vol_faces[local_face]; build face_conn from face vertex indices when length matches n_face_nodes, else keep first-n legacy fallback for high-order volume/surface mismatch - README.md: document Quad4/Seg2/Cook entry points; clarify mesh I/O is not in core load path (see src/io/README.md) --- src/mesh/README.md | 16 ++-- src/mesh/mesh.jl | 42 +++++---- src/mesh/structured.jl | 200 ++++++++++++++++++++++++++++++++++++++++- 3 files changed, 234 insertions(+), 24 deletions(-) diff --git a/src/mesh/README.md b/src/mesh/README.md index 3d7ba45..416918b 100644 --- a/src/mesh/README.md +++ b/src/mesh/README.md @@ -12,10 +12,11 @@ optimisations used by the rest of the package. type, constructors with validation, the inverse connectivity (`node -> elements`) needed by node-based assembly, helpers for sets and surface extraction. -- `structured.jl` — `create_structured_box_mesh`, - `create_unit_cube_mesh`, `create_cantilever_mesh`, - `create_thin_plate_mesh`. Boundary node sets are populated - automatically (`:xmin`, `:xmax`, …). +- `structured.jl` — `create_structured_box_mesh` (`Hex8` or `Quad4`), + `create_structured_line_mesh` (`Seg2`), `create_unit_cube_mesh`, + `create_cook_membrane_mesh` (`Quad4`, Cook skew panel), + `create_cantilever_mesh`, `create_thin_plate_mesh`. Boundary node sets are + populated automatically (`:xmin`, `:xmax`, …). - `refine.jl` — `LongestEdgeBisection` and the `refine` entry point. Used for h-convergence studies. @@ -66,9 +67,10 @@ longest edge of each element, doubling the element count per level. ## I/O -Mesh import lives in `src/io/` (currently the self-contained Gmsh -reader). VTK / XDMF output is not implemented in the new path; the -legacy results writers under `src/legacy/` cover existing tests. +Mesh import from external tools (Gmsh, Netgen, Abaqus, …) is **not** part of +the core load path; use optional extensions or companion packages that build +`Mesh{…}` (see `src/io/README.md`). VTK / XDMF output is not implemented in the +new path; legacy results writers under `src/legacy/` cover older tests. ## Related code diff --git a/src/mesh/mesh.jl b/src/mesh/mesh.jl index c2cce4a..fed090e 100644 --- a/src/mesh/mesh.jl +++ b/src/mesh/mesh.jl @@ -311,7 +311,7 @@ surface_topology(::Type{Hexahedron{20}}) = Quadrilateral{8} surface_topology(::Type{Hexahedron{27}}) = Quadrilateral{9} """ - extract_surface(mesh::Mesh{N,T}, face_set::Symbol) -> Mesh{Nface,FaceT} + extract_surface(mesh::Mesh{N,T}, face_set::Symbol, local_face::Int=1) -> Mesh{Nface,FaceT} Extract a surface mesh from volume elements. The boundary-face topology is looked up via [`surface_topology`](@ref), which currently supports @@ -319,21 +319,24 @@ is looked up via [`surface_topology`](@ref), which currently supports # Arguments - `mesh::Mesh{N,T}`: Volume mesh. -- `face_set::Symbol`: Element set whose elements should contribute their - boundary face. The set must already exist in `mesh.element_sets`. +- `face_set::Symbol`: Element set whose elements should contribute a face + (one surface element per volume element in the set). +- `local_face::Int`: Which volume-local face to use, `1 … nfaces(T)` in the + order returned by `faces(T())` (e.g. `1` is the first `Face` for `Hex8`, + typically the ``z = z_{\\min}`` bilinear face in structured meshes). # Returns - `Mesh{Nface,FaceT}` whose nodes alias the volume mesh's node array. -# Limitations -The current implementation takes the first `nnodes(FaceT)` connectivity -entries of each volume element as a face. This is correct only when the -caller has already arranged volume elements so that the first -`nnodes(FaceT)` nodes form the boundary face (e.g. extruded prism layers). -Topology-aware face extraction using the per-volume face tables remains a -known TODO; see the corresponding session log. +When the face table lists exactly `nnodes(surface_topology(T))` volume-local +vertices (true for `Tet4` / `Hex8`), connectivity is built from those indices +(topology-aware). For higher-order volumes whose face description lists only +corners while the surface topology needs more nodes (e.g. `Hex20` with +`Quad8`), the implementation falls back to taking the first `nnodes(FaceT)` +entries of the volume connectivity (legacy behaviour; may not match a true +quadratic face). """ -function extract_surface(mesh::Mesh{N,T}, face_set::Symbol) where {N,T<:AbstractTopology{N}} +function extract_surface(mesh::Mesh{N,T}, face_set::Symbol, local_face::Int=1) where {N,T<:AbstractTopology{N}} if !hasmethod(surface_topology, Tuple{Type{T}}) error("extract_surface: no surface_topology trait defined for $T. " * "Supported volume topologies: Tetrahedron{4|10}, Hexahedron{8|20|27}.") @@ -344,13 +347,22 @@ function extract_surface(mesh::Mesh{N,T}, face_set::Symbol) where {N,T<:Abstract @assert haskey(mesh.element_sets, face_set) "Element set $face_set not found" face_elements = mesh.element_sets[face_set] - # Extract surface connectivity. NOTE: this assumes the first - # `n_face_nodes` entries of each volume element form a face. See the - # function docstring; topology-aware extraction is on the backlog. + vol_faces = faces(T()) + nf = length(vol_faces) + if !(1 ≤ local_face ≤ nf) + error("extract_surface: local_face must be in 1:$nf for topology $T, got $local_face") + end + loc = vol_faces[local_face].vertices + use_topology = length(loc) == n_face_nodes + surface_conn = NTuple{n_face_nodes,UInt32}[] for elem_id in face_elements elem_conn = mesh.connectivity[elem_id] - face_conn = ntuple(i -> elem_conn[i], n_face_nodes) + if use_topology + face_conn = ntuple(j -> elem_conn[Int(loc[j])]::UInt32, n_face_nodes) + else + face_conn = ntuple(i -> elem_conn[i], n_face_nodes) + end push!(surface_conn, face_conn) end diff --git a/src/mesh/structured.jl b/src/mesh/structured.jl index 7e9626a..2d7545e 100644 --- a/src/mesh/structured.jl +++ b/src/mesh/structured.jl @@ -1,5 +1,5 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +# SPDX-FileCopyrightText: 2015-2026 Jukka Aho +# SPDX-License-Identifier: MIT """ create_structured_box_mesh(::Type{Hex8}; @@ -320,3 +320,199 @@ function create_thin_plate_mesh( ymin=0.0, ymax=width, ny=ny, zmin=0.0, zmax=thickness, nz=nz) end + +""" + create_structured_line_mesh(::Type{Seg2}; + x0::Float64=0.0, x1::Float64=1.0, nx::Int=1, + y::Float64=0.0, z::Float64=0.0) -> Mesh{Seg2} + +Structured 1D mesh of `Seg2` elements embedded in 3D space along **x** +(from `x0` to `x1` with `nx` elements). Nodes use `Vec{3,Float64}` with +constant `y` and `z` (defaults 0), matching the rest of the mesh stack. + +# Node sets +- `:all`, `:xmin` (first node), `:xmax` (last node) +""" +function create_structured_line_mesh( + ::Type{Seg2}; + x0::Float64=0.0, + x1::Float64=1.0, + nx::Int=1, + y::Float64=0.0, + z::Float64=0.0, +) + @assert nx ≥ 1 "nx must be ≥ 1" + @assert x1 > x0 "x1 must be > x0" + + xs = range(x0, x1, length=nx + 1) + nodes = Vec{3,Float64}[Vec(x, y, z) for x in xs] + + connectivity = NTuple{2,UInt32}[] + for i in 1:nx + push!(connectivity, (UInt32(i), UInt32(i + 1))) + end + + element_sets = Dict{Symbol,Set{UInt32}}(:all => Set(UInt32(1):UInt32(nx))) + node_sets = Dict{Symbol,Set{UInt32}}( + :all => Set(UInt32(1):UInt32(length(nodes))), + :xmin => Set((UInt32(1),)), + :xmax => Set((UInt32(length(nodes)),)), + ) + + return Mesh{Seg2}(nodes, connectivity; element_sets=element_sets, node_sets=node_sets) +end + +""" + create_structured_box_mesh(::Type{Quad4}; + xmin=0.0, xmax=1.0, nx::Int=1, + ymin=0.0, ymax=1.0, ny::Int=1, + z::Float64=0.0) -> Mesh{Quad4} + +Structured tensor-product mesh of bilinear quads in the **xy** plane at +fixed `z` (default 0). Connectivity matches the `Hex8` bottom-face winding +(`Quad4` reference order). + +# Node sets +- `:all`, `:xmin`, `:xmax`, `:ymin`, `:ymax` (same naming as structured bricks; + no `z` faces for a single-layer 2D mesh) +""" +function create_structured_box_mesh( + ::Type{Quad4}; + xmin::Float64=0.0, + xmax::Float64=1.0, + nx::Int=1, + ymin::Float64=0.0, + ymax::Float64=1.0, + ny::Int=1, + z::Float64=0.0, +) + @assert nx ≥ 1 && ny ≥ 1 "nx and ny must be ≥ 1" + @assert xmax > xmin && ymax > ymin "box bounds must be increasing" + + xs = range(xmin, xmax, length=nx + 1) + ys = range(ymin, ymax, length=ny + 1) + + nodes = Vec{3,Float64}[] + for j in 1:(ny + 1) + for i in 1:(nx + 1) + push!(nodes, Vec(xs[i], ys[j], z)) + end + end + + node_index(i::Int, j::Int) = UInt32((j - 1) * (nx + 1) + i) + + connectivity = NTuple{4,UInt32}[] + for j in 1:ny, i in 1:nx + n1 = node_index(i, j) + n2 = node_index(i + 1, j) + n3 = node_index(i + 1, j + 1) + n4 = node_index(i, j + 1) + push!(connectivity, (n1, n2, n3, n4)) + end + + element_sets = Dict{Symbol,Set{UInt32}}(:all => Set(UInt32(1):UInt32(length(connectivity)))) + + node_sets = Dict{Symbol,Set{UInt32}}() + node_sets[:all] = Set(UInt32(1):UInt32(length(nodes))) + + xmin_nodes = Set{UInt32}() + xmax_nodes = Set{UInt32}() + for j in 1:(ny + 1) + push!(xmin_nodes, node_index(1, j)) + push!(xmax_nodes, node_index(nx + 1, j)) + end + node_sets[:xmin] = xmin_nodes + node_sets[:xmax] = xmax_nodes + + ymin_nodes = Set{UInt32}() + ymax_nodes = Set{UInt32}() + for i in 1:(nx + 1) + push!(ymin_nodes, node_index(i, 1)) + push!(ymax_nodes, node_index(i, ny + 1)) + end + node_sets[:ymin] = ymin_nodes + node_sets[:ymax] = ymax_nodes + + return Mesh{Quad4}(nodes, connectivity; element_sets=element_sets, node_sets=node_sets) +end + +""" + create_cook_membrane_mesh(::Type{Quad4}, nx::Int, ny::Int; scale::Float64=1e-3) -> Mesh{Quad4} + +Structured `nx × ny` bilinear `Quad4` mesh on **Cook's membrane** reference geometry +(classical skew panel benchmark). + +The physical domain is the convex quadrilateral with corners (millimetres before scaling): +(0, 0), (48, 44), (48, 60), (0, 44). + +Mapped from the parametric unit square `[0, 1]²` with the same bilinear map as a single +`Quad4` element: corners SW, SE, NE, NW at `(ξ, η) ∈ {(0,0), (1,0), (1,1), (0,1)}`. +Coordinates are multiplied by `scale` (default `1e-3`, i.e. millimetres to metres). + +# Node sets + +Same labels as `create_structured_box_mesh(Quad4; xmin=0, xmax=1, ymin=0, ymax=1, …)`: +`:xmin` is the clamped Cook **left** edge, `:xmax` the **right** edge (typical traction side), +plus `:ymin`, `:ymax`, `:all`. + +# References + +Cook, R. D., *Improved Two-Dimensional Finite Element*, Journal of Applied Mechanics +**40** (1973). The geometry and skew bending/shear mode are widely reproduced in FE +textbooks and software validation suites. +""" +function create_cook_membrane_mesh(::Type{Quad4}, nx::Int, ny::Int; scale::Float64=1e-3) + @assert nx ≥ 1 && ny ≥ 1 "nx and ny must be ≥ 1" + + function _cook_xy_mm(ξ::Float64, η::Float64) + x = 48.0 * ξ + y = ξ * (1.0 - η) * 44.0 + ξ * η * 60.0 + (1.0 - ξ) * η * 44.0 + return scale * x, scale * y + end + + nodes = Vec{3,Float64}[] + for j in 1:(ny + 1) + for i in 1:(nx + 1) + ξ = (i - 1) / nx + η = (j - 1) / ny + x, y = _cook_xy_mm(ξ, η) + push!(nodes, Vec(x, y, 0.0)) + end + end + + node_index(i::Int, j::Int) = UInt32((j - 1) * (nx + 1) + i) + + connectivity = NTuple{4,UInt32}[] + for j in 1:ny, i in 1:nx + n1 = node_index(i, j) + n2 = node_index(i + 1, j) + n3 = node_index(i + 1, j + 1) + n4 = node_index(i, j + 1) + push!(connectivity, (n1, n2, n3, n4)) + end + + element_sets = Dict{Symbol,Set{UInt32}}(:all => Set(UInt32(1):UInt32(length(connectivity)))) + + node_sets = Dict{Symbol,Set{UInt32}}() + node_sets[:all] = Set(UInt32(1):UInt32(length(nodes))) + + xmin_nodes = Set{UInt32}() + xmax_nodes = Set{UInt32}() + for j in 1:(ny + 1) + push!(xmin_nodes, node_index(1, j)) + push!(xmax_nodes, node_index(nx + 1, j)) + end + node_sets[:xmin] = xmin_nodes + node_sets[:xmax] = xmax_nodes + + ymin_nodes = Set{UInt32}() + ymax_nodes = Set{UInt32}() + for i in 1:(nx + 1) + push!(ymin_nodes, node_index(i, 1)) + push!(ymax_nodes, node_index(i, ny + 1)) + end + node_sets[:ymin] = ymin_nodes + node_sets[:ymax] = ymax_nodes + + return Mesh{Quad4}(nodes, connectivity; element_sets=element_sets, node_sets=node_sets) +end