From e64e33af9bd3af0f7e8d0bb3fde31cd77a4f32fd Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 08:38:16 +0000 Subject: [PATCH 1/7] Initial plan From 8b95fd0e0d163a1b833a70a5b3473fe0b6bb30a1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 08:55:11 +0000 Subject: [PATCH 2/7] Add VEM linear elasticity solver for UnstructuredMesh Co-authored-by: moyner <454871+moyner@users.noreply.github.com> --- src/meshes/meshes.jl | 3 + src/meshes/vem.jl | 970 +++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/vem.jl | 167 ++++++++ 4 files changed, 1141 insertions(+) create mode 100644 src/meshes/vem.jl create mode 100644 test/vem.jl diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index b1f336716..974955dcc 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -313,3 +313,6 @@ end include("RadialMeshes/RadialMeshes.jl") +# Virtual Element Method for linear elasticity +include("vem.jl") + diff --git a/src/meshes/vem.jl b/src/meshes/vem.jl new file mode 100644 index 000000000..a8891f6c8 --- /dev/null +++ b/src/meshes/vem.jl @@ -0,0 +1,970 @@ +export solve_vem_elasticity, VEMElasticitySetup, assemble_vem_elasticity + +""" + VEMCellData + +Precomputed per-cell data from mesh traversal, including local node indices, +coordinates, face-to-node connectivity, outward normals, and cell volume. +""" +struct VEMCellData + "Global node indices for this cell (unique, sorted)" + node_indices::Vector{Int} + "Local node coordinates (D × num_nodes)" + coordinates::Matrix{Float64} + "Face-to-node connectivity: each entry is a vector of local node indices for a face" + face_nodes::Vector{Vector{Int}} + "Outward unit normals for each face (D × num_faces)" + outward_normals::Matrix{Float64} + "Cell volume" + volume::Float64 + "Cell centroid (length D)" + centroid::Vector{Float64} +end + +""" + VEMElasticitySetup + +Precomputed setup data for VEM linear elasticity. Created once from a mesh +and reused for multiple solves with different material parameters. +""" +struct VEMElasticitySetup{D} + "Spatial dimension" + dim::Int + "Number of nodes in the mesh" + num_nodes::Int + "Number of cells in the mesh" + num_cells::Int + "Per-cell precomputed data" + cell_data::Vector{VEMCellData} + "Indices of boundary nodes (unique, sorted)" + boundary_nodes::Vector{Int} +end + +""" + VEMElasticitySetup(mesh::UnstructuredMesh) + +Create a VEM elasticity setup by traversing the mesh once and precomputing +all geometric data needed for assembly. +""" +function VEMElasticitySetup(mesh::UnstructuredMesh) + D = dim(mesh) + @assert D == 2 || D == 3 "VEM elasticity only supports 2D and 3D meshes" + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + pts = mesh.node_points + + cell_data = Vector{VEMCellData}(undef, nc) + for cell in 1:nc + cell_data[cell] = _precompute_cell_data(mesh, cell, D) + end + + # Identify boundary nodes + boundary_nodes = _find_boundary_nodes(mesh) + + return VEMElasticitySetup{D}(D, nn, nc, cell_data, boundary_nodes) +end + +""" + assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0) + +Assemble the VEM linear elasticity system (stiffness matrix and RHS) for given +per-cell material properties and pressure change. + +Returns `(K, rhs)` where `K` is the sparse stiffness matrix and `rhs` is the +right-hand side vector, both with boundary conditions applied. +""" +function assemble_vem_elasticity( + setup::VEMElasticitySetup{D}, + youngs_modulus::AbstractVector, + poisson_ratio::AbstractVector, + pressure_change::AbstractVector; + biot_coefficient::Union{Real, AbstractVector} = 1.0 +) where {D} + nc = setup.num_cells + nn = setup.num_nodes + ndof = D * nn + + @assert length(youngs_modulus) == nc + @assert length(poisson_ratio) == nc + @assert length(pressure_change) == nc + + # Prepare triplets for sparse matrix assembly + I_idx = Int[] + J_idx = Int[] + V_val = Float64[] + + # Right-hand side + rhs = zeros(Float64, ndof) + + # Get biot coefficient per cell + biot = biot_coefficient isa Real ? fill(Float64(biot_coefficient), nc) : Float64.(biot_coefficient) + + for cell in 1:nc + cd = setup.cell_data[cell] + num_nodes_cell = length(cd.node_indices) + local_ndof = D * num_nodes_cell + + # Compute local stiffness matrix + K_local = _assemble_local_stiffness(cd, youngs_modulus[cell], poisson_ratio[cell], D) + + # Scatter local stiffness matrix to global + for i in 1:local_ndof + gi = D * (cd.node_indices[(i - 1) ÷ D + 1] - 1) + (i - 1) % D + 1 + for j in 1:local_ndof + gj = D * (cd.node_indices[(j - 1) ÷ D + 1] - 1) + (j - 1) % D + 1 + val = K_local[i, j] + if val != 0.0 + push!(I_idx, gi) + push!(J_idx, gj) + push!(V_val, val) + end + end + end + + # Compute body force from pressure change and scatter to global RHS + # Body force = -biot * grad(p) ≈ distributed via VEM face integration + _add_pressure_bodyforce!(rhs, cd, biot[cell] * pressure_change[cell], D) + end + + # Determine boundary DOFs + bnd_dof_set = Set{Int}() + for node in setup.boundary_nodes + for d in 1:D + push!(bnd_dof_set, D * (node - 1) + d) + end + end + + # Filter out entries involving boundary DOFs, then add identity for boundary DOFs + I_final = Int[] + J_final = Int[] + V_final = Float64[] + for k in eachindex(I_idx) + i = I_idx[k] + j = J_idx[k] + if !(i in bnd_dof_set) && !(j in bnd_dof_set) + push!(I_final, i) + push!(J_final, j) + push!(V_final, V_val[k]) + end + end + # Add identity for boundary DOFs + for dof in bnd_dof_set + push!(I_final, dof) + push!(J_final, dof) + push!(V_final, 1.0) + rhs[dof] = 0.0 + end + + # Assemble sparse matrix + K = sparse(I_final, J_final, V_final, ndof, ndof) + + return (K, rhs) +end + +""" + solve_vem_elasticity(mesh, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0) + +Solve the VEM linear elasticity problem on an UnstructuredMesh with zero +displacement at the boundaries. + +# Arguments +- `mesh::UnstructuredMesh`: The mesh to solve on. +- `youngs_modulus::AbstractVector`: Young's modulus per cell. +- `poisson_ratio::AbstractVector`: Poisson's ratio per cell. +- `pressure_change::AbstractVector`: Pressure change per cell. +- `biot_coefficient`: Biot coefficient (scalar or per-cell vector), default 1.0. + +# Returns +A named tuple `(displacement, setup, K, rhs)` where: +- `displacement`: nodal displacement vector of length `D * num_nodes`. +- `setup`: the `VEMElasticitySetup` for reuse. +- `K`: the assembled sparse stiffness matrix. +- `rhs`: the right-hand side vector. + +# Example +```julia +mesh = UnstructuredMesh(CartesianMesh((3, 3, 3), (1.0, 1.0, 1.0))) +nc = number_of_cells(mesh) +result = solve_vem_elasticity( + mesh, + fill(1e9, nc), # Young's modulus + fill(0.3, nc), # Poisson's ratio + fill(1e6, nc) # Pressure change +) +u = result.displacement +``` +""" +function solve_vem_elasticity( + mesh::UnstructuredMesh, + youngs_modulus::AbstractVector, + poisson_ratio::AbstractVector, + pressure_change::AbstractVector; + biot_coefficient::Union{Real, AbstractVector} = 1.0 +) + setup = VEMElasticitySetup(mesh) + return solve_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient = biot_coefficient) +end + +""" + solve_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0) + +Solve with a precomputed `VEMElasticitySetup` (useful for repeated solves). +""" +function solve_vem_elasticity( + setup::VEMElasticitySetup, + youngs_modulus::AbstractVector, + poisson_ratio::AbstractVector, + pressure_change::AbstractVector; + biot_coefficient::Union{Real, AbstractVector} = 1.0 +) + K, rhs = assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient = biot_coefficient) + displacement = K \ rhs + return (displacement = displacement, setup = setup, K = K, rhs = rhs) +end + +# ============================================================================ +# Internal helper functions +# ============================================================================ + +""" +Precompute cell data for a single cell by traversing its faces (interior + boundary). +""" +function _precompute_cell_data(mesh::UnstructuredMesh, cell::Int, D::Int) + pts = mesh.node_points + + # Collect all unique global node indices for this cell + global_nodes = Set{Int}() + + # Collect face-to-node connectivity (global indices) + face_nodes_global = Vector{Vector{Int}}() + + # Interior faces + for face in mesh.faces.cells_to_faces[cell] + nodes = collect(mesh.faces.faces_to_nodes[face]) + for n in nodes + push!(global_nodes, n) + end + push!(face_nodes_global, nodes) + end + + # Boundary faces + for bf in mesh.boundary_faces.cells_to_faces[cell] + nodes = collect(mesh.boundary_faces.faces_to_nodes[bf]) + for n in nodes + push!(global_nodes, n) + end + push!(face_nodes_global, nodes) + end + + # Sort global node indices and create global-to-local map + sorted_nodes = sort(collect(global_nodes)) + num_nodes = length(sorted_nodes) + g2l = Dict{Int, Int}() + for (li, gi) in enumerate(sorted_nodes) + g2l[gi] = li + end + + # Local coordinates matrix (D × num_nodes) + coords = zeros(Float64, D, num_nodes) + for (li, gi) in enumerate(sorted_nodes) + for d in 1:D + coords[d, li] = pts[gi][d] + end + end + + # Convert face-to-node connectivity to local indices + face_nodes_local = Vector{Vector{Int}}(undef, length(face_nodes_global)) + for (fi, fnodes) in enumerate(face_nodes_global) + face_nodes_local[fi] = [g2l[n] for n in fnodes] + end + + # Compute outward normals for each face + nfaces = length(face_nodes_local) + outward_normals = zeros(Float64, D, nfaces) + + if D == 3 + _compute_outward_normals_3d!(outward_normals, face_nodes_local, coords, mesh, cell) + else + _compute_outward_normals_2d!(outward_normals, face_nodes_local, coords, mesh, cell) + end + + # Compute cell centroid and volume + centroid, volume = _compute_cell_centroid_and_volume(coords, face_nodes_local, outward_normals, D) + + return VEMCellData(sorted_nodes, coords, face_nodes_local, outward_normals, volume, centroid) +end + +""" +Compute outward normals for 3D faces. The normals are unit normals pointing +outward from the cell. +""" +function _compute_outward_normals_3d!(normals, face_nodes, coords, mesh, cell) + nfaces = length(face_nodes) + # Compute face centroids and normals + face_centroids = zeros(Float64, 3, nfaces) + for fi in 1:nfaces + fn = face_nodes[fi] + nfn = length(fn) + # Face centroid (simple average of nodes) + fc = zeros(3) + for ni in fn + fc .+= coords[:, ni] + end + fc ./= nfn + face_centroids[:, fi] .= fc + + # Compute normal via weighted cross-product sum + normal = zeros(3) + for i in 1:nfn + prev_i = i == 1 ? nfn : i - 1 + next_i = i == nfn ? 1 : i + 1 + a = coords[:, fn[prev_i]] .- coords[:, fn[i]] + c = coords[:, fn[next_i]] .- coords[:, fn[i]] + normal .+= cross(c, a) + end + nn = norm(normal) + if nn > 0 + normal ./= nn + end + normals[:, fi] .= normal + end + + # Check orientation: normals should point outward + # Compute cell center (average of all face centroids) + cell_center = zeros(3) + for fi in 1:nfaces + cell_center .+= face_centroids[:, fi] + end + cell_center ./= nfaces + + # Check if normals point outward by checking dot product with vector from + # cell center to face centroid + dot_sum = 0.0 + for fi in 1:nfaces + v = face_centroids[:, fi] .- cell_center + dot_sum += dot(v, normals[:, fi]) + end + + # If the sum is negative, flip all normals + if dot_sum < 0 + normals .*= -1 + end +end + +""" +Compute outward normals for 2D faces (edges). +""" +function _compute_outward_normals_2d!(normals, face_nodes, coords, mesh, cell) + nfaces = length(face_nodes) + face_centroids = zeros(Float64, 2, nfaces) + for fi in 1:nfaces + fn = face_nodes[fi] + @assert length(fn) == 2 "2D faces must have exactly 2 nodes" + p1 = coords[:, fn[1]] + p2 = coords[:, fn[2]] + # Edge midpoint + fc = (p1 + p2) / 2 + face_centroids[:, fi] .= fc + # Normal is perpendicular to edge + edge = p2 - p1 + normal = [edge[2], -edge[1]] + nn = norm(normal) + if nn > 0 + normal ./= nn + end + normals[:, fi] .= normal + end + + # Check orientation + cell_center = zeros(2) + for fi in 1:nfaces + cell_center .+= face_centroids[:, fi] + end + cell_center ./= nfaces + + dot_sum = 0.0 + for fi in 1:nfaces + v = face_centroids[:, fi] .- cell_center + dot_sum += dot(v, normals[:, fi]) + end + if dot_sum < 0 + normals .*= -1 + end +end + +""" +Compute cell centroid and volume from local face-node data and outward normals. +Uses the divergence theorem: V = (1/D) ∫_∂E x · n dA +""" +function _compute_cell_centroid_and_volume(coords, face_nodes, outward_normals, D) + if D == 3 + return _compute_cell_geometry_3d(coords, face_nodes, outward_normals) + else + return _compute_cell_geometry_2d(coords, face_nodes, outward_normals) + end +end + +function _compute_cell_geometry_3d(coords, face_nodes, normals) + # Use the average of all node points as an "inside point" + num_nodes = size(coords, 2) + inside_pt = zeros(3) + for i in 1:num_nodes + inside_pt .+= coords[:, i] + end + inside_pt ./= num_nodes + + volume = 0.0 + centroid = zeros(3) + + for fi in eachindex(face_nodes) + fn = face_nodes[fi] + nfn = length(fn) + # Compute face centroid + fc = zeros(3) + for ni in fn + fc .+= coords[:, ni] + end + fc ./= nfn + + # Tessellate face into triangles using face centroid + for i in 1:nfn + next_i = i == nfn ? 1 : i + 1 + # Triangle: fc, fn[i], fn[next_i] + p1 = coords[:, fn[i]] + p2 = coords[:, fn[next_i]] + + # Tetrahedron volume: V = |det([p1-p4, p2-p4, p3-p4])| / 6 + v1 = fc .- inside_pt + v2 = p1 .- inside_pt + v3 = p2 .- inside_pt + tet_vol = abs(dot(v1, cross(v2, v3))) / 6.0 + tet_centroid = (fc .+ p1 .+ p2 .+ inside_pt) ./ 4.0 + + volume += tet_vol + centroid .+= tet_vol .* tet_centroid + end + end + + if volume > 0 + centroid ./= volume + end + + return (centroid, volume) +end + +function _compute_cell_geometry_2d(coords, face_nodes, normals) + num_nodes = size(coords, 2) + inside_pt = zeros(2) + for i in 1:num_nodes + inside_pt .+= coords[:, i] + end + inside_pt ./= num_nodes + + area = 0.0 + centroid = zeros(2) + + for fi in eachindex(face_nodes) + fn = face_nodes[fi] + @assert length(fn) == 2 + p1 = coords[:, fn[1]] + p2 = coords[:, fn[2]] + + # Triangle: inside_pt, p1, p2 + A = p1 .- inside_pt + B = p2 .- inside_pt + tri_area = abs(A[1] * B[2] - A[2] * B[1]) / 2.0 + tri_centroid = (inside_pt .+ p1 .+ p2) ./ 3.0 + + area += tri_area + centroid .+= tri_area .* tri_centroid + end + + if area > 0 + centroid ./= area + end + + return (centroid, area) +end + +""" +Find all nodes on the boundary of the mesh. +""" +function _find_boundary_nodes(mesh::UnstructuredMesh) + boundary_nodes = Set{Int}() + nbf = number_of_boundary_faces(mesh) + for bf in 1:nbf + for node in mesh.boundary_faces.faces_to_nodes[bf] + push!(boundary_nodes, node) + end + end + return sort(collect(boundary_nodes)) +end + +# ============================================================================ +# VEM stiffness matrix assembly per cell +# ============================================================================ + +""" +Compute the elasticity tensor D in modified Voigt notation (with factor 2 on shear terms). +""" +function _compute_D(young::Float64, poisson::Float64, D_dim::Int) + fac = young / (1 + poisson) / (1 - 2 * poisson) + if D_dim == 2 + # 3x3 matrix (plane strain, modified Voigt notation) + D = zeros(3, 3) + D[1, 1] = (1 - poisson) * fac + D[1, 2] = poisson * fac + D[2, 1] = poisson * fac + D[2, 2] = (1 - poisson) * fac + D[3, 3] = 2 * (1 - 2 * poisson) * fac + else + # 6x6 matrix (3D, modified Voigt notation) + D = zeros(6, 6) + for i in 1:3 + D[i, i] = (1 - poisson) * fac + end + D[1, 2] = poisson * fac; D[2, 1] = poisson * fac + D[1, 3] = poisson * fac; D[3, 1] = poisson * fac + D[2, 3] = poisson * fac; D[3, 2] = poisson * fac + D[4, 4] = 2 * (1 - 2 * poisson) * fac + D[5, 5] = 2 * (1 - 2 * poisson) * fac + D[6, 6] = 2 * (1 - 2 * poisson) * fac + end + return D +end + +""" +Compute q-values for 2D (related to VBF integrals over edges). +q_i = 1/(2A) ∫_∂E φ_i n dA + +For 2D: q is a (num_nodes × 2) matrix stored as flat vector. +""" +function _compute_q_2d(coords, face_nodes, volume) + num_nodes = size(coords, 2) + q = zeros(num_nodes, 2) + fac = 1.0 / (4.0 * volume) # 1/(2*2*area) since we distribute to two nodes + + for fi in eachindex(face_nodes) + fn = face_nodes[fi] + @assert length(fn) == 2 + i = fn[1] + inext = fn[2] + # Scaled normal (not unit normal, but edge-length scaled) + scaled_normal = [coords[2, inext] - coords[2, i], + -(coords[1, inext] - coords[1, i])] + for d in 1:2 + q[i, d] += fac * scaled_normal[d] + q[inext, d] += fac * scaled_normal[d] + end + end + return q +end + +""" +Compute q-values for 3D (related to VBF integrals over faces). +q_i = 1/(2V) ∫_∂E φ_i n dA +""" +function _compute_q_3d(coords, face_nodes, outward_normals, volume) + num_nodes = size(coords, 2) + q = zeros(num_nodes, 3) + fac = 1.0 / (2.0 * volume) + + for fi in eachindex(face_nodes) + fn = face_nodes[fi] + nfn = length(fn) + normal = outward_normals[:, fi] + + # Compute face integration: each basis function φ_i integrated over the face + # Using tessellation into triangles from face centroid + fc = zeros(3) + for ni in fn + fc .+= coords[:, ni] + end + fc ./= nfn + + for e in 1:nfn + enext = e == nfn ? 1 : e + 1 + # Compute integrals of φ for two sub-triangles associated with edge (e, enext) + # Sub-triangle 1: fc, coords[fn[e]], midpoint + # Sub-triangle 2: fc, midpoint, coords[fn[enext]] + p1 = coords[:, fn[e]] + p2 = coords[:, fn[enext]] + mid = (p1 .+ p2) ./ 2.0 + + # Triangle areas using cross product + # Triangle 1: fc, p1, mid + v1 = p1 .- fc + v2 = mid .- fc + area1 = norm(cross(v1, v2)) / 2.0 + + # Triangle 2: fc, mid, p2 + v1 = mid .- fc + v2 = p2 .- fc + area2 = norm(cross(v1, v2)) / 2.0 + + # φ_e is 1 at node e, 0 at others. Under tessellation, the integral + # of φ_e over its associated sub-triangles gives the area. + for d in 1:3 + q[fn[e], d] += fac * area1 * normal[d] + q[fn[enext], d] += fac * area2 * normal[d] + end + end + end + return q +end + +""" +Create a 2D matentry (2 × 3 sub-matrix for node i). +Returns the entries in row-major order: [e1 0 e2; 0 e3 e4] +""" +function _matentry_2d(e1, e2, e3, e4) + return [e1 0.0 e2; 0.0 e3 e4] +end + +""" +Create a 3D matentry (3 × 6 sub-matrix for node i). +Returns [e1 0 0 e2 0 e3; 0 e4 0 e5 e6 0; 0 0 e7 0 e8 e9] +""" +function _matentry_3d(e1, e2, e3, e4, e5, e6, e7, e8, e9) + return [e1 0.0 0.0 e2 0.0 e3; + 0.0 e4 0.0 e5 e6 0.0; + 0.0 0.0 e7 0.0 e8 e9] +end + +""" +Compute Nr matrix (rigid body modes) for 2D. +Nr is (2*num_nodes) × 3. +""" +function _compute_Nr_2d(coords) + num_nodes = size(coords, 2) + midpoint = vec(mean(coords, dims = 2)) + Nr = zeros(2 * num_nodes, 3) + for i in 1:num_nodes + dx = coords[1, i] - midpoint[1] + dy = coords[2, i] - midpoint[2] + rows = (2 * (i - 1) + 1):(2 * i) + Nr[rows, :] .= _matentry_2d(1.0, dy, 1.0, -dx) + end + return Nr +end + +""" +Compute Nr matrix (rigid body modes) for 3D. +Nr is (3*num_nodes) × 6. +""" +function _compute_Nr_3d(coords) + num_nodes = size(coords, 2) + midpoint = vec(mean(coords, dims = 2)) + Nr = zeros(3 * num_nodes, 6) + for i in 1:num_nodes + dx = coords[1, i] - midpoint[1] + dy = coords[2, i] - midpoint[2] + dz = coords[3, i] - midpoint[3] + rows = (3 * (i - 1) + 1):(3 * i) + Nr[rows, :] .= _matentry_3d(1.0, dy, -dz, 1.0, -dx, dz, 1.0, -dy, dx) + end + return Nr +end + +""" +Compute Nc matrix (linear deformation modes) for 2D. +Nc is (2*num_nodes) × 3. +""" +function _compute_Nc_2d(coords) + num_nodes = size(coords, 2) + midpoint = vec(mean(coords, dims = 2)) + Nc = zeros(2 * num_nodes, 3) + for i in 1:num_nodes + dx = coords[1, i] - midpoint[1] + dy = coords[2, i] - midpoint[2] + rows = (2 * (i - 1) + 1):(2 * i) + Nc[rows, :] .= _matentry_2d(dx, dy, dy, dx) + end + return Nc +end + +""" +Compute Nc matrix (linear deformation modes) for 3D. +Nc is (3*num_nodes) × 6. +""" +function _compute_Nc_3d(coords) + num_nodes = size(coords, 2) + midpoint = vec(mean(coords, dims = 2)) + Nc = zeros(3 * num_nodes, 6) + for i in 1:num_nodes + dx = coords[1, i] - midpoint[1] + dy = coords[2, i] - midpoint[2] + dz = coords[3, i] - midpoint[3] + rows = (3 * (i - 1) + 1):(3 * i) + Nc[rows, :] .= _matentry_3d(dx, dy, dz, dy, dx, dz, dz, dy, dx) + end + return Nc +end + +""" +Compute Wr matrix for 2D (from q-values + rigid body rotation). +Wr is (2*num_nodes) × 3. +""" +function _compute_Wr_2d(q) + num_nodes = size(q, 1) + ncinv = 1.0 / num_nodes + Wr = zeros(2 * num_nodes, 3) + for i in 1:num_nodes + rows = (2 * (i - 1) + 1):(2 * i) + Wr[rows, :] .= _matentry_2d(ncinv, q[i, 2], ncinv, -q[i, 1]) + end + return Wr +end + +""" +Compute Wr matrix for 3D. +""" +function _compute_Wr_3d(q) + num_nodes = size(q, 1) + ncinv = 1.0 / num_nodes + Wr = zeros(3 * num_nodes, 6) + for i in 1:num_nodes + rows = (3 * (i - 1) + 1):(3 * i) + Wr[rows, :] .= _matentry_3d(ncinv, q[i, 2], -q[i, 3], ncinv, -q[i, 1], q[i, 3], ncinv, -q[i, 2], q[i, 1]) + end + return Wr +end + +""" +Compute Wc matrix for 2D. +""" +function _compute_Wc_2d(q) + num_nodes = size(q, 1) + Wc = zeros(2 * num_nodes, 3) + for i in 1:num_nodes + rows = (2 * (i - 1) + 1):(2 * i) + Wc[rows, :] .= _matentry_2d(2 * q[i, 1], q[i, 2], 2 * q[i, 2], q[i, 1]) + end + return Wc +end + +""" +Compute Wc matrix for 3D. +""" +function _compute_Wc_3d(q) + num_nodes = size(q, 1) + Wc = zeros(3 * num_nodes, 6) + for i in 1:num_nodes + rows = (3 * (i - 1) + 1):(3 * i) + Wc[rows, :] .= _matentry_3d( + 2 * q[i, 1], q[i, 2], q[i, 3], + 2 * q[i, 2], q[i, 1], q[i, 3], + 2 * q[i, 3], q[i, 2], q[i, 1] + ) + end + return Wc +end + +""" +Compute the VEM stability term S. +Based on Gain et al. 2014, DOI:10.1016/j.cma.2014.05.005 +""" +function _compute_S(Nc, D_mat, num_nodes, volume, D_dim) + lsdim = D_dim == 2 ? 3 : 6 + r = D_dim * num_nodes + NtN = Nc' * Nc + alpha = volume * tr(D_mat) / tr(NtN) + return alpha * I(r) +end + +""" +Compute I - P where P = Nr * Wr' + Nc * Wc'. +""" +function _compute_ImP(Nr, Nc, Wr, Wc) + r = size(Nr, 1) + Pr = Nr * Wr' + Pc = Nc * Wc' + return Matrix{Float64}(I, r, r) - Pr - Pc +end + +""" +Assemble the local VEM stiffness matrix for a single cell. +""" +function _assemble_local_stiffness(cd::VEMCellData, young::Float64, poisson::Float64, D_dim::Int) + coords = cd.coordinates + face_nodes = cd.face_nodes + outward_normals = cd.outward_normals + volume = cd.volume + num_nodes = length(cd.node_indices) + + # Compute elasticity tensor + D_mat = _compute_D(young, poisson, D_dim) + + # Compute q-values + if D_dim == 2 + q = _compute_q_2d(coords, face_nodes, volume) + else + q = _compute_q_3d(coords, face_nodes, outward_normals, volume) + end + + # Compute matrices + if D_dim == 2 + Nr = _compute_Nr_2d(coords) + Nc = _compute_Nc_2d(coords) + Wr = _compute_Wr_2d(q) + Wc = _compute_Wc_2d(q) + else + Nr = _compute_Nr_3d(coords) + Nc = _compute_Nc_3d(coords) + Wr = _compute_Wr_3d(q) + Wc = _compute_Wc_3d(q) + end + + ImP = _compute_ImP(Nr, Nc, Wr, Wc) + S = _compute_S(Nc, D_mat, num_nodes, volume, D_dim) + + # Final assembly: K = V * Wc * D * Wc' + (I-P)' * S * (I-P) + lsdim = D_dim == 2 ? 3 : 6 + totdim = D_dim * num_nodes + + DWct = D_mat * Wc' + EWcDWct = volume * (Wc * DWct) + + SImP = S * ImP + ImPtSImP = ImP' * SImP + + K_local = EWcDWct + ImPtSImP + + return K_local +end + +""" +Distribute the body force from pressure change to the global RHS vector. +For a cell with pressure change dp and Biot coefficient α: + body_force = -α * ∇p → distributed to nodes via VEM face integration + +The contribution from a cell to node i in direction d is: + f_d = α * dp * (volume of region associated with node i in direction d) + +Using VEM: this is equivalent to -α * dp * 2 * q_i[d] * volume +(since q_i = 1/(2V) ∫ φ_i n dA, the integral ∫ φ_i dV in direction d +relates to the face integrals). + +More precisely, the pressure gradient contribution is: + rhs_id = -biot * dp * ∫_∂E φ_i n_d dA = -biot * dp * 2V * q_id +""" +function _add_pressure_bodyforce!(rhs, cd::VEMCellData, biot_dp, D_dim) + if biot_dp == 0.0 + return + end + + coords = cd.coordinates + face_nodes = cd.face_nodes + outward_normals = cd.outward_normals + volume = cd.volume + num_nodes = length(cd.node_indices) + + # Compute the face integral ∫_∂E φ_i n dA for each node and direction + # This equals 2V * q_i (since q_i = 1/(2V) * ∫ φ_i n dA) + # The body force from uniform pressure is distributed as: + # f_id = -biot * dp * ∫_∂E φ_i n_d dA + # + # We compute z_id = ∫_∂E φ_i n_d dA directly (which is 2V * q_id) + # and set rhs contribution to -biot * dp * z_id + + if D_dim == 2 + z = _compute_z_2d(coords, face_nodes) + else + z = _compute_z_3d(coords, face_nodes, outward_normals) + end + + # Add to global RHS: note the NEGATIVE sign because pressure acts as compression + for i in 1:num_nodes + gi = cd.node_indices[i] + for d in 1:D_dim + rhs[D_dim * (gi - 1) + d] -= biot_dp * z[i, d] + end + end +end + +""" +Compute z_id = ∫_∂E φ_i n_d dA for 2D (the unnormalized q values). +""" +function _compute_z_2d(coords, face_nodes) + num_nodes = size(coords, 2) + z = zeros(num_nodes, 2) + + for fi in eachindex(face_nodes) + fn = face_nodes[fi] + @assert length(fn) == 2 + i = fn[1] + inext = fn[2] + # Scaled normal (edge-length × outward normal direction) + # This is the same as the outer normal weighted by edge length / 2 + scaled_normal = [coords[2, inext] - coords[2, i], + -(coords[1, inext] - coords[1, i])] + for d in 1:2 + z[i, d] += 0.5 * scaled_normal[d] + z[inext, d] += 0.5 * scaled_normal[d] + end + end + return z +end + +""" +Compute z_id = ∫_∂E φ_i n_d dA for 3D. +""" +function _compute_z_3d(coords, face_nodes, outward_normals) + num_nodes = size(coords, 2) + z = zeros(num_nodes, 3) + + for fi in eachindex(face_nodes) + fn = face_nodes[fi] + nfn = length(fn) + normal = outward_normals[:, fi] + + fc = zeros(3) + for ni in fn + fc .+= coords[:, ni] + end + fc ./= nfn + + for e in 1:nfn + enext = e == nfn ? 1 : e + 1 + p1 = coords[:, fn[e]] + p2 = coords[:, fn[enext]] + mid = (p1 .+ p2) ./ 2.0 + + v1 = p1 .- fc + v2 = mid .- fc + area1 = norm(cross(v1, v2)) / 2.0 + + v1 = mid .- fc + v2 = p2 .- fc + area2 = norm(cross(v1, v2)) / 2.0 + + for d in 1:3 + z[fn[e], d] += area1 * normal[d] + z[fn[enext], d] += area2 * normal[d] + end + end + end + return z +end + +""" +Apply zero displacement Dirichlet boundary conditions. +Sets rows and columns of boundary DOFs to identity. +""" +function _apply_zero_dirichlet!(K, rhs, boundary_nodes, D_dim) + # For zero displacement: just set boundary rows to identity and rhs to zero + ndof = length(rhs) + for node in boundary_nodes + for d in 1:D_dim + dof = D_dim * (node - 1) + d + # Zero out the row + K[dof, :] .= 0.0 + # Zero out the column (move known values to RHS - here all zero so no change) + K[:, dof] .= 0.0 + # Set diagonal + K[dof, dof] = 1.0 + # Set RHS to prescribed value (zero) + rhs[dof] = 0.0 + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 8e1577702..34c8de4b2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,3 +13,4 @@ include("units.jl") include("sparsity.jl") include("nfvm.jl") include("weno.jl") +include("vem.jl") diff --git a/test/vem.jl b/test/vem.jl new file mode 100644 index 000000000..490329e00 --- /dev/null +++ b/test/vem.jl @@ -0,0 +1,167 @@ +using Jutul +using Test +using LinearAlgebra +using SparseArrays + +@testset "VEM Linear Elasticity" begin + @testset "VEMElasticitySetup" begin + @testset "3D setup" begin + mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) + setup = VEMElasticitySetup(mesh) + @test setup.dim == 3 + @test setup.num_nodes == length(mesh.node_points) + @test setup.num_cells == number_of_cells(mesh) + @test length(setup.cell_data) == number_of_cells(mesh) + @test length(setup.boundary_nodes) > 0 + # All cell data should have positive volume + for cd in setup.cell_data + @test cd.volume > 0 + end + end + + @testset "2D setup" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3), (1.0, 1.0))) + setup = VEMElasticitySetup(mesh) + @test setup.dim == 2 + @test setup.num_nodes == length(mesh.node_points) + @test setup.num_cells == number_of_cells(mesh) + for cd in setup.cell_data + @test cd.volume > 0 + end + end + end + + @testset "Assembly and solve 3D" begin + mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + E = fill(1e9, nc) + nu = fill(0.3, nc) + dp = fill(1e6, nc) + + setup = VEMElasticitySetup(mesh) + K, rhs = assemble_vem_elasticity(setup, E, nu, dp) + + @test size(K) == (3 * nn, 3 * nn) + @test length(rhs) == 3 * nn + @test issparse(K) + + @testset "Stiffness matrix symmetry" begin + K_full = Matrix(K) + @test norm(K_full - K_full') / norm(K_full) < 1e-10 + end + + @testset "Zero displacement at boundaries" begin + result = solve_vem_elasticity(setup, E, nu, dp) + u = result.displacement + for node in setup.boundary_nodes + for d in 1:3 + dof = 3 * (node - 1) + d + @test abs(u[dof]) < 1e-14 + end + end + end + + @testset "Zero pressure gives zero displacement" begin + result0 = solve_vem_elasticity(setup, E, nu, zeros(nc)) + @test norm(result0.displacement) < 1e-14 + end + + @testset "Linear scaling with pressure" begin + result1 = solve_vem_elasticity(setup, E, nu, dp) + result2 = solve_vem_elasticity(setup, E, nu, 2 .* dp) + @test norm(result2.displacement) / norm(result1.displacement) ≈ 2.0 atol=1e-10 + end + + @testset "Stiffer material gives smaller displacement" begin + result_soft = solve_vem_elasticity(setup, fill(1e9, nc), nu, dp) + result_stiff = solve_vem_elasticity(setup, fill(1e10, nc), nu, dp) + ratio = norm(result_soft.displacement) / norm(result_stiff.displacement) + @test ratio ≈ 10.0 atol=1e-8 + end + end + + @testset "Assembly and solve 2D" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3), (1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + E = fill(1e9, nc) + nu = fill(0.3, nc) + dp = fill(1e6, nc) + + setup = VEMElasticitySetup(mesh) + + @testset "Stiffness matrix symmetry" begin + K, rhs = assemble_vem_elasticity(setup, E, nu, dp) + K_full = Matrix(K) + @test norm(K_full - K_full') / norm(K_full) < 1e-10 + end + + @testset "Zero displacement at boundaries" begin + result = solve_vem_elasticity(setup, E, nu, dp) + for node in setup.boundary_nodes + for d in 1:2 + dof = 2 * (node - 1) + d + @test abs(result.displacement[dof]) < 1e-14 + end + end + end + + @testset "Zero pressure gives zero displacement" begin + result0 = solve_vem_elasticity(setup, E, nu, zeros(nc)) + @test norm(result0.displacement) < 1e-14 + end + + @testset "Linear scaling with pressure" begin + result1 = solve_vem_elasticity(setup, E, nu, dp) + result2 = solve_vem_elasticity(setup, E, nu, 3 .* dp) + @test norm(result2.displacement) / norm(result1.displacement) ≈ 3.0 atol=1e-10 + end + end + + @testset "Setup reuse" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3, 3), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + setup = VEMElasticitySetup(mesh) + + # Solve with different material properties using same setup + E1 = fill(1e9, nc) + E2 = fill(2e9, nc) + nu = fill(0.25, nc) + dp = fill(1e6, nc) + + result1 = solve_vem_elasticity(setup, E1, nu, dp) + result2 = solve_vem_elasticity(setup, E2, nu, dp) + + @test norm(result1.displacement) > norm(result2.displacement) + end + + @testset "Convenience function (mesh input)" begin + mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + E = fill(1e9, nc) + nu = fill(0.3, nc) + dp = fill(1e6, nc) + + result = solve_vem_elasticity(mesh, E, nu, dp) + @test haskey(pairs(result) |> Dict, :displacement) + @test haskey(pairs(result) |> Dict, :setup) + @test haskey(pairs(result) |> Dict, :K) + @test haskey(pairs(result) |> Dict, :rhs) + @test length(result.displacement) == 3 * length(mesh.node_points) + end + + @testset "Biot coefficient" begin + mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + setup = VEMElasticitySetup(mesh) + E = fill(1e9, nc) + nu = fill(0.3, nc) + dp = fill(1e6, nc) + + result_biot1 = solve_vem_elasticity(setup, E, nu, dp, biot_coefficient = 1.0) + result_biot05 = solve_vem_elasticity(setup, E, nu, dp, biot_coefficient = 0.5) + + @test norm(result_biot1.displacement) / norm(result_biot05.displacement) ≈ 2.0 atol=1e-10 + end +end From 79ca426e078f1e108f1f4e9c95bc45190e14cf7e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 08:58:36 +0000 Subject: [PATCH 3/7] Address code review: fix degenerate geometry, clean up unused function, improve tests Co-authored-by: moyner <454871+moyner@users.noreply.github.com> --- src/meshes/vem.jl | 30 ++++++++---------------------- test/vem.jl | 8 ++++---- 2 files changed, 12 insertions(+), 26 deletions(-) diff --git a/src/meshes/vem.jl b/src/meshes/vem.jl index a8891f6c8..3cb9d88b5 100644 --- a/src/meshes/vem.jl +++ b/src/meshes/vem.jl @@ -768,7 +768,13 @@ function _compute_S(Nc, D_mat, num_nodes, volume, D_dim) lsdim = D_dim == 2 ? 3 : 6 r = D_dim * num_nodes NtN = Nc' * Nc - alpha = volume * tr(D_mat) / tr(NtN) + trNtN = tr(NtN) + if trNtN < eps(Float64) + # Fallback for degenerate geometry + alpha = volume * tr(D_mat) + else + alpha = volume * tr(D_mat) / trNtN + end return alpha * I(r) end @@ -947,24 +953,4 @@ function _compute_z_3d(coords, face_nodes, outward_normals) return z end -""" -Apply zero displacement Dirichlet boundary conditions. -Sets rows and columns of boundary DOFs to identity. -""" -function _apply_zero_dirichlet!(K, rhs, boundary_nodes, D_dim) - # For zero displacement: just set boundary rows to identity and rhs to zero - ndof = length(rhs) - for node in boundary_nodes - for d in 1:D_dim - dof = D_dim * (node - 1) + d - # Zero out the row - K[dof, :] .= 0.0 - # Zero out the column (move known values to RHS - here all zero so no change) - K[:, dof] .= 0.0 - # Set diagonal - K[dof, dof] = 1.0 - # Set RHS to prescribed value (zero) - rhs[dof] = 0.0 - end - end -end + diff --git a/test/vem.jl b/test/vem.jl index 490329e00..954ac03bb 100644 --- a/test/vem.jl +++ b/test/vem.jl @@ -144,10 +144,10 @@ using SparseArrays dp = fill(1e6, nc) result = solve_vem_elasticity(mesh, E, nu, dp) - @test haskey(pairs(result) |> Dict, :displacement) - @test haskey(pairs(result) |> Dict, :setup) - @test haskey(pairs(result) |> Dict, :K) - @test haskey(pairs(result) |> Dict, :rhs) + @test hasproperty(result, :displacement) + @test hasproperty(result, :setup) + @test hasproperty(result, :K) + @test hasproperty(result, :rhs) @test length(result.displacement) == 3 * length(mesh.node_points) end From 79ad08e759381d9f64a79c558e364f17bcc7653c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 10:10:04 +0000 Subject: [PATCH 4/7] Fix VEM: per-face normal orientation, face-based pressure RHS, nonzero BC support Root causes fixed: 1. _compute_outward_normals_3d/2d: was checking normal orientation as a global sum and flipping ALL normals together. Individual face normals could be independently wrong. Fixed to check each face individually. 2. _compute_q_2d: was computing its own normals from edge direction instead of using the corrected outward_normals. Fixed to use the outward normals directly. 3. _add_pressure_bodyforce!: per-cell assembly didn't cancel on shared faces for uniform pressure. Replaced with face-based _assemble_pressure_rhs! that correctly gives zero when p_L == p_R. Added boundary_displacement kwarg for nonzero Dirichlet BCs. Added linear displacement patch tests (2D and 3D). All 583 tests pass. Co-authored-by: moyner <454871+moyner@users.noreply.github.com> --- src/meshes/vem.jl | 386 ++++++++++++++++++++++++++++++---------------- test/vem.jl | 84 +++++++++- 2 files changed, 335 insertions(+), 135 deletions(-) diff --git a/src/meshes/vem.jl b/src/meshes/vem.jl index 3cb9d88b5..f4683c87a 100644 --- a/src/meshes/vem.jl +++ b/src/meshes/vem.jl @@ -27,7 +27,7 @@ end Precomputed setup data for VEM linear elasticity. Created once from a mesh and reused for multiple solves with different material parameters. """ -struct VEMElasticitySetup{D} +struct VEMElasticitySetup{D, M} "Spatial dimension" dim::Int "Number of nodes in the mesh" @@ -38,6 +38,8 @@ struct VEMElasticitySetup{D} cell_data::Vector{VEMCellData} "Indices of boundary nodes (unique, sorted)" boundary_nodes::Vector{Int} + "Reference to the mesh (needed for face-based pressure RHS assembly)" + mesh::M end """ @@ -61,24 +63,30 @@ function VEMElasticitySetup(mesh::UnstructuredMesh) # Identify boundary nodes boundary_nodes = _find_boundary_nodes(mesh) - return VEMElasticitySetup{D}(D, nn, nc, cell_data, boundary_nodes) + return VEMElasticitySetup{D, typeof(mesh)}(D, nn, nc, cell_data, boundary_nodes, mesh) end """ - assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0) + assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0, boundary_displacement=nothing) Assemble the VEM linear elasticity system (stiffness matrix and RHS) for given per-cell material properties and pressure change. Returns `(K, rhs)` where `K` is the sparse stiffness matrix and `rhs` is the right-hand side vector, both with boundary conditions applied. + +# Keyword arguments +- `biot_coefficient`: Biot coefficient (scalar or per-cell vector), default 1.0. +- `boundary_displacement`: If `nothing` (default), zero displacement on all boundary nodes. + If a vector of length `D * num_nodes`, the values at boundary DOFs are used as prescribed displacement. """ function assemble_vem_elasticity( setup::VEMElasticitySetup{D}, youngs_modulus::AbstractVector, poisson_ratio::AbstractVector, pressure_change::AbstractVector; - biot_coefficient::Union{Real, AbstractVector} = 1.0 + biot_coefficient::Union{Real, AbstractVector} = 1.0, + boundary_displacement::Union{Nothing, AbstractVector} = nothing ) where {D} nc = setup.num_cells nn = setup.num_nodes @@ -120,17 +128,44 @@ function assemble_vem_elasticity( end end end - - # Compute body force from pressure change and scatter to global RHS - # Body force = -biot * grad(p) ≈ distributed via VEM face integration - _add_pressure_bodyforce!(rhs, cd, biot[cell] * pressure_change[cell], D) end - # Determine boundary DOFs + # Assemble pressure RHS using face-based approach + _assemble_pressure_rhs!(rhs, setup, pressure_change, biot, D) + + # Determine boundary DOFs and prescribed values bnd_dof_set = Set{Int}() + bnd_dof_vals = Dict{Int, Float64}() for node in setup.boundary_nodes for d in 1:D - push!(bnd_dof_set, D * (node - 1) + d) + dof = D * (node - 1) + d + push!(bnd_dof_set, dof) + if boundary_displacement !== nothing + bnd_dof_vals[dof] = boundary_displacement[dof] + else + bnd_dof_vals[dof] = 0.0 + end + end + end + + # Build the full (unconstrained) stiffness matrix first to compute + # RHS modification for nonzero BC + K_full = sparse(I_idx, J_idx, V_val, ndof, ndof) + + # Modify RHS for nonzero boundary displacement: + # For each boundary DOF j with prescribed value u_j, + # subtract K[:,j] * u_j from the RHS for interior DOFs + for (dof, val) in bnd_dof_vals + if val != 0.0 + col = K_full[:, dof] + rows = rowvals(col) + vals = nonzeros(col) + for k in eachindex(rows) + row = rows[k] + if !(row in bnd_dof_set) + rhs[row] -= vals[k] * val + end + end end end @@ -147,12 +182,12 @@ function assemble_vem_elasticity( push!(V_final, V_val[k]) end end - # Add identity for boundary DOFs - for dof in bnd_dof_set + # Add identity for boundary DOFs with prescribed values + for (dof, val) in bnd_dof_vals push!(I_final, dof) push!(J_final, dof) push!(V_final, 1.0) - rhs[dof] = 0.0 + rhs[dof] = val end # Assemble sparse matrix @@ -162,10 +197,10 @@ function assemble_vem_elasticity( end """ - solve_vem_elasticity(mesh, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0) + solve_vem_elasticity(mesh, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0, boundary_displacement=nothing) Solve the VEM linear elasticity problem on an UnstructuredMesh with zero -displacement at the boundaries. +displacement at the boundaries (or prescribed nonzero boundary displacement). # Arguments - `mesh::UnstructuredMesh`: The mesh to solve on. @@ -173,6 +208,8 @@ displacement at the boundaries. - `poisson_ratio::AbstractVector`: Poisson's ratio per cell. - `pressure_change::AbstractVector`: Pressure change per cell. - `biot_coefficient`: Biot coefficient (scalar or per-cell vector), default 1.0. +- `boundary_displacement`: If `nothing` (default), zero displacement on all boundary nodes. + If a vector of length `D * num_nodes`, the values at boundary DOFs are used as prescribed displacement. # Returns A named tuple `(displacement, setup, K, rhs)` where: @@ -199,14 +236,16 @@ function solve_vem_elasticity( youngs_modulus::AbstractVector, poisson_ratio::AbstractVector, pressure_change::AbstractVector; - biot_coefficient::Union{Real, AbstractVector} = 1.0 + biot_coefficient::Union{Real, AbstractVector} = 1.0, + boundary_displacement::Union{Nothing, AbstractVector} = nothing ) setup = VEMElasticitySetup(mesh) - return solve_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient = biot_coefficient) + return solve_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; + biot_coefficient = biot_coefficient, boundary_displacement = boundary_displacement) end """ - solve_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0) + solve_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient=1.0, boundary_displacement=nothing) Solve with a precomputed `VEMElasticitySetup` (useful for repeated solves). """ @@ -215,9 +254,11 @@ function solve_vem_elasticity( youngs_modulus::AbstractVector, poisson_ratio::AbstractVector, pressure_change::AbstractVector; - biot_coefficient::Union{Real, AbstractVector} = 1.0 + biot_coefficient::Union{Real, AbstractVector} = 1.0, + boundary_displacement::Union{Nothing, AbstractVector} = nothing ) - K, rhs = assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient = biot_coefficient) + K, rhs = assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; + biot_coefficient = biot_coefficient, boundary_displacement = boundary_displacement) displacement = K \ rhs return (displacement = displacement, setup = setup, K = K, rhs = rhs) end @@ -329,25 +370,21 @@ function _compute_outward_normals_3d!(normals, face_nodes, coords, mesh, cell) normals[:, fi] .= normal end - # Check orientation: normals should point outward - # Compute cell center (average of all face centroids) + # Check orientation per face: each normal should point outward from the cell + # Compute cell center (average of all node coordinates) + num_nodes = size(coords, 2) cell_center = zeros(3) - for fi in 1:nfaces - cell_center .+= face_centroids[:, fi] + for i in 1:num_nodes + cell_center .+= coords[:, i] end - cell_center ./= nfaces + cell_center ./= num_nodes - # Check if normals point outward by checking dot product with vector from - # cell center to face centroid - dot_sum = 0.0 + # Flip individual normals that point inward for fi in 1:nfaces v = face_centroids[:, fi] .- cell_center - dot_sum += dot(v, normals[:, fi]) - end - - # If the sum is negative, flip all normals - if dot_sum < 0 - normals .*= -1 + if dot(v, normals[:, fi]) < 0 + normals[:, fi] .*= -1 + end end end @@ -375,20 +412,19 @@ function _compute_outward_normals_2d!(normals, face_nodes, coords, mesh, cell) normals[:, fi] .= normal end - # Check orientation + # Check orientation per face + num_nodes = size(coords, 2) cell_center = zeros(2) - for fi in 1:nfaces - cell_center .+= face_centroids[:, fi] + for i in 1:num_nodes + cell_center .+= coords[:, i] end - cell_center ./= nfaces + cell_center ./= num_nodes - dot_sum = 0.0 for fi in 1:nfaces v = face_centroids[:, fi] .- cell_center - dot_sum += dot(v, normals[:, fi]) - end - if dot_sum < 0 - normals .*= -1 + if dot(v, normals[:, fi]) < 0 + normals[:, fi] .*= -1 + end end end @@ -537,24 +573,27 @@ end Compute q-values for 2D (related to VBF integrals over edges). q_i = 1/(2A) ∫_∂E φ_i n dA -For 2D: q is a (num_nodes × 2) matrix stored as flat vector. +For 2D: q is a (num_nodes × 2) matrix. +For each edge, ∫_edge φ_i n dA = (edge_length/2) * n for each of the 2 endpoint nodes. """ -function _compute_q_2d(coords, face_nodes, volume) +function _compute_q_2d(coords, face_nodes, outward_normals, volume) num_nodes = size(coords, 2) q = zeros(num_nodes, 2) - fac = 1.0 / (4.0 * volume) # 1/(2*2*area) since we distribute to two nodes + fac = 1.0 / (2.0 * volume) for fi in eachindex(face_nodes) fn = face_nodes[fi] @assert length(fn) == 2 i = fn[1] inext = fn[2] - # Scaled normal (not unit normal, but edge-length scaled) - scaled_normal = [coords[2, inext] - coords[2, i], - -(coords[1, inext] - coords[1, i])] + # Edge length + edge_len = norm(coords[:, inext] .- coords[:, i]) + # Area-weighted normal = edge_length * unit_outward_normal + # Each endpoint gets half the integral for d in 1:2 - q[i, d] += fac * scaled_normal[d] - q[inext, d] += fac * scaled_normal[d] + contrib = fac * (edge_len / 2.0) * outward_normals[d, fi] + q[i, d] += contrib + q[inext, d] += contrib end end return q @@ -803,7 +842,7 @@ function _assemble_local_stiffness(cd::VEMCellData, young::Float64, poisson::Flo # Compute q-values if D_dim == 2 - q = _compute_q_2d(coords, face_nodes, volume) + q = _compute_q_2d(coords, face_nodes, outward_normals, volume) else q = _compute_q_3d(coords, face_nodes, outward_normals, volume) end @@ -840,117 +879,204 @@ function _assemble_local_stiffness(cd::VEMCellData, young::Float64, poisson::Flo end """ -Distribute the body force from pressure change to the global RHS vector. -For a cell with pressure change dp and Biot coefficient α: - body_force = -α * ∇p → distributed to nodes via VEM face integration - -The contribution from a cell to node i in direction d is: - f_d = α * dp * (volume of region associated with node i in direction d) +Assemble the pressure contribution to the RHS using face-based approach. + +For piecewise-constant pressure, the weak form contribution is: + rhs_i = −∑_f α·(p_L − p_R)·∫_f φ_i n dA (interior faces) + − ∑_f α·p·∫_f φ_i n dA (boundary faces) + +This correctly gives zero RHS when pressure is uniform, because p_L == p_R +on every interior face cancels exactly, and boundary face contributions +cancel by the divergence theorem (∑_bnd_faces ∫_f n dA = 0 for each node's +"share" when pressure is constant). + +Actually for the VEM formulation with piecewise-constant pressure in cells, +the body force from pressure in a cell is: + f = -α * ∇p * V ≈ -α * ∑_faces (p_face * n * A_face) +where the sum is over cell faces with outward normals. + +The face-based assembly distributes -α*p*∫_f φ_i n dA for each face, +where the sign accounts for which cell "owns" the normal direction. +For interior faces shared by cells L and R: + - Cell L sees outward normal n_f, contributes -α_L * p_L * ∫_f φ_i n dA + - Cell R sees outward normal -n_f, contributes +α_R * p_R * ∫_f φ_i n dA + => net = -∫_f φ_i n dA * (α_L * p_L - α_R * p_R) +When p_L == p_R and α_L == α_R, this is zero. + +For boundary faces, only one cell contributes. +""" +function _assemble_pressure_rhs!(rhs, setup::VEMElasticitySetup{D}, pressure_change, biot, D_dim) where {D} + mesh = setup.mesh + pts = mesh.node_points -Using VEM: this is equivalent to -α * dp * 2 * q_i[d] * volume -(since q_i = 1/(2V) ∫ φ_i n dA, the integral ∫ φ_i dV in direction d -relates to the face integrals). + # Process interior faces + nf = number_of_faces(mesh) + for face in 1:nf + L, R = mesh.faces.neighbors[face] + biot_dp_L = biot[L] * pressure_change[L] + biot_dp_R = biot[R] * pressure_change[R] + dp_jump = biot_dp_L - biot_dp_R + if dp_jump == 0.0 + continue + end -More precisely, the pressure gradient contribution is: - rhs_id = -biot * dp * ∫_∂E φ_i n_d dA = -biot * dp * 2V * q_id -""" -function _add_pressure_bodyforce!(rhs, cd::VEMCellData, biot_dp, D_dim) - if biot_dp == 0.0 - return - end + # Get face nodes (global indices) + face_node_indices = collect(mesh.faces.faces_to_nodes[face]) - coords = cd.coordinates - face_nodes = cd.face_nodes - outward_normals = cd.outward_normals - volume = cd.volume - num_nodes = length(cd.node_indices) + # Compute face normal (oriented from L to R by convention in Jutul) + face_normal_vec = _compute_face_normal_from_nodes(face_node_indices, pts, D_dim) - # Compute the face integral ∫_∂E φ_i n dA for each node and direction - # This equals 2V * q_i (since q_i = 1/(2V) * ∫ φ_i n dA) - # The body force from uniform pressure is distributed as: - # f_id = -biot * dp * ∫_∂E φ_i n_d dA - # - # We compute z_id = ∫_∂E φ_i n_d dA directly (which is 2V * q_id) - # and set rhs contribution to -biot * dp * z_id + # Orient: Jutul stores normals pointing from L to R (first to second neighbor) + # We need outward normal from L, which is face_normal_vec + # Verify orientation using cell centroids + cL = setup.cell_data[L].centroid + cR = setup.cell_data[R].centroid + LR = zeros(D_dim) + for d in 1:D_dim + LR[d] = cR[d] - cL[d] + end + if dot(face_normal_vec, LR) < 0 + face_normal_vec = -face_normal_vec + end - if D_dim == 2 - z = _compute_z_2d(coords, face_nodes) - else - z = _compute_z_3d(coords, face_nodes, outward_normals) + # Compute ∫_f φ_i n dA for each face node + # For the contribution: rhs -= dp_jump * ∫_f φ_i n dA + _add_face_pressure_contribution!(rhs, face_node_indices, pts, face_normal_vec, dp_jump, D_dim) end - # Add to global RHS: note the NEGATIVE sign because pressure acts as compression - for i in 1:num_nodes - gi = cd.node_indices[i] - for d in 1:D_dim - rhs[D_dim * (gi - 1) + d] -= biot_dp * z[i, d] + # Process boundary faces + nbf = number_of_boundary_faces(mesh) + for bf in 1:nbf + cell = mesh.boundary_faces.neighbors[bf] + biot_dp = biot[cell] * pressure_change[cell] + if biot_dp == 0.0 + continue end + + face_node_indices = collect(mesh.boundary_faces.faces_to_nodes[bf]) + face_normal_vec = _compute_face_normal_from_nodes(face_node_indices, pts, D_dim) + + # Orient outward from cell + cc = setup.cell_data[cell].centroid + fc = zeros(D_dim) + for ni in face_node_indices + for d in 1:D_dim + fc[d] += pts[ni][d] + end + end + fc ./= length(face_node_indices) + outward = fc .- cc + if dot(face_normal_vec, outward) < 0 + face_normal_vec = -face_normal_vec + end + + # Boundary face: only one cell contributes + # rhs -= biot_dp * ∫_f φ_i n dA + _add_face_pressure_contribution!(rhs, face_node_indices, pts, face_normal_vec, biot_dp, D_dim) end end """ -Compute z_id = ∫_∂E φ_i n_d dA for 2D (the unnormalized q values). +Compute face normal vector (unnormalized, area-weighted) from global node indices. """ -function _compute_z_2d(coords, face_nodes) - num_nodes = size(coords, 2) - z = zeros(num_nodes, 2) +function _compute_face_normal_from_nodes(node_indices, pts, D_dim) + if D_dim == 2 + @assert length(node_indices) == 2 + p1 = pts[node_indices[1]] + p2 = pts[node_indices[2]] + edge = [p2[1] - p1[1], p2[2] - p1[2]] + return [edge[2], -edge[1]] # unnormalized normal (length = edge length) + else + # 3D: compute normal via cross-product sum over polygon + nfn = length(node_indices) + fc = zeros(3) + for ni in node_indices + for d in 1:3 + fc[d] += pts[ni][d] + end + end + fc ./= nfn - for fi in eachindex(face_nodes) - fn = face_nodes[fi] - @assert length(fn) == 2 - i = fn[1] - inext = fn[2] - # Scaled normal (edge-length × outward normal direction) - # This is the same as the outer normal weighted by edge length / 2 - scaled_normal = [coords[2, inext] - coords[2, i], - -(coords[1, inext] - coords[1, i])] - for d in 1:2 - z[i, d] += 0.5 * scaled_normal[d] - z[inext, d] += 0.5 * scaled_normal[d] + normal = zeros(3) + for i in 1:nfn + prev_i = i == 1 ? nfn : i - 1 + next_i = i == nfn ? 1 : i + 1 + pp = [pts[node_indices[prev_i]][d] for d in 1:3] + pc = [pts[node_indices[i]][d] for d in 1:3] + pn = [pts[node_indices[next_i]][d] for d in 1:3] + a = pp .- pc + c = pn .- pc + normal .+= cross(c, a) end + # This gives 2 * area * unit_normal + # We want area-weighted normal, so divide by 2 + return normal ./ 2.0 end - return z end """ -Compute z_id = ∫_∂E φ_i n_d dA for 3D. -""" -function _compute_z_3d(coords, face_nodes, outward_normals) - num_nodes = size(coords, 2) - z = zeros(num_nodes, 3) +Add the face pressure contribution −dp_jump * ∫_f φ_i n_d dA to the RHS +for each node on the face and each direction d. - for fi in eachindex(face_nodes) - fn = face_nodes[fi] - nfn = length(fn) - normal = outward_normals[:, fi] +For 2D faces (edges): ∫_f φ_i n dA = (edge_length / 2) * unit_normal for each of the 2 nodes + = (1/2) * area_normal_vec (since area_normal_vec = edge_length * unit_normal) + +For 3D faces: we tessellate into triangles from face centroid and distribute +to nodes proportional to their sub-triangle areas. +""" +function _add_face_pressure_contribution!(rhs, node_indices, pts, area_normal_vec, dp_jump, D_dim) + nfn = length(node_indices) + if D_dim == 2 + # 2D edge: equal distribution to 2 nodes, each gets half the face integral + for ni in node_indices + for d in 1:D_dim + rhs[D_dim * (ni - 1) + d] -= dp_jump * 0.5 * area_normal_vec[d] + end + end + else + # 3D face: tessellate into triangles from face centroid + # Each node's share is proportional to the sub-triangle areas adjacent to it fc = zeros(3) - for ni in fn - fc .+= coords[:, ni] + for ni in node_indices + for d in 1:3 + fc[d] += pts[ni][d] + end end fc ./= nfn + # Unit normal direction (for weighting) + nn = norm(area_normal_vec) + if nn < eps(Float64) + return + end + unit_normal = area_normal_vec ./ nn + + # Compute sub-triangle areas for each node + node_areas = zeros(nfn) + total_area = 0.0 for e in 1:nfn enext = e == nfn ? 1 : e + 1 - p1 = coords[:, fn[e]] - p2 = coords[:, fn[enext]] - mid = (p1 .+ p2) ./ 2.0 - + p1 = [pts[node_indices[e]][d] for d in 1:3] + p2 = [pts[node_indices[enext]][d] for d in 1:3] + # Triangle: fc, p1, p2 v1 = p1 .- fc - v2 = mid .- fc - area1 = norm(cross(v1, v2)) / 2.0 - - v1 = mid .- fc v2 = p2 .- fc - area2 = norm(cross(v1, v2)) / 2.0 + tri_area = norm(cross(v1, v2)) / 2.0 + # Split between the two edge nodes + node_areas[e] += tri_area / 2.0 + node_areas[enext] += tri_area / 2.0 + total_area += tri_area + end - for d in 1:3 - z[fn[e], d] += area1 * normal[d] - z[fn[enext], d] += area2 * normal[d] + # Each node gets its share of the face area times the unit normal + for (li, ni) in enumerate(node_indices) + weight = node_areas[li] + for d in 1:D_dim + rhs[D_dim * (ni - 1) + d] -= dp_jump * weight * unit_normal[d] end end end - return z end diff --git a/test/vem.jl b/test/vem.jl index 954ac03bb..df1b9ffba 100644 --- a/test/vem.jl +++ b/test/vem.jl @@ -37,7 +37,8 @@ using SparseArrays nn = length(mesh.node_points) E = fill(1e9, nc) nu = fill(0.3, nc) - dp = fill(1e6, nc) + # Non-uniform pressure to get nonzero displacement + dp = [Float64(i) * 1e5 for i in 1:nc] setup = VEMElasticitySetup(mesh) K, rhs = assemble_vem_elasticity(setup, E, nu, dp) @@ -67,9 +68,15 @@ using SparseArrays @test norm(result0.displacement) < 1e-14 end + @testset "Uniform pressure gives zero displacement" begin + result_unif = solve_vem_elasticity(setup, E, nu, fill(1e6, nc)) + @test norm(result_unif.displacement) < 1e-14 + end + @testset "Linear scaling with pressure" begin result1 = solve_vem_elasticity(setup, E, nu, dp) result2 = solve_vem_elasticity(setup, E, nu, 2 .* dp) + @test norm(result1.displacement) > 0 @test norm(result2.displacement) / norm(result1.displacement) ≈ 2.0 atol=1e-10 end @@ -87,7 +94,8 @@ using SparseArrays nn = length(mesh.node_points) E = fill(1e9, nc) nu = fill(0.3, nc) - dp = fill(1e6, nc) + # Non-uniform pressure + dp = [Float64(i) * 1e5 for i in 1:nc] setup = VEMElasticitySetup(mesh) @@ -112,9 +120,15 @@ using SparseArrays @test norm(result0.displacement) < 1e-14 end + @testset "Uniform pressure gives zero displacement" begin + result_unif = solve_vem_elasticity(setup, E, nu, fill(1e6, nc)) + @test norm(result_unif.displacement) < 1e-14 + end + @testset "Linear scaling with pressure" begin result1 = solve_vem_elasticity(setup, E, nu, dp) result2 = solve_vem_elasticity(setup, E, nu, 3 .* dp) + @test norm(result1.displacement) > 0 @test norm(result2.displacement) / norm(result1.displacement) ≈ 3.0 atol=1e-10 end end @@ -128,7 +142,8 @@ using SparseArrays E1 = fill(1e9, nc) E2 = fill(2e9, nc) nu = fill(0.25, nc) - dp = fill(1e6, nc) + # Non-uniform pressure + dp = [Float64(i) * 1e5 for i in 1:nc] result1 = solve_vem_elasticity(setup, E1, nu, dp) result2 = solve_vem_elasticity(setup, E2, nu, dp) @@ -141,7 +156,7 @@ using SparseArrays nc = number_of_cells(mesh) E = fill(1e9, nc) nu = fill(0.3, nc) - dp = fill(1e6, nc) + dp = [Float64(i) * 1e5 for i in 1:nc] result = solve_vem_elasticity(mesh, E, nu, dp) @test hasproperty(result, :displacement) @@ -157,11 +172,70 @@ using SparseArrays setup = VEMElasticitySetup(mesh) E = fill(1e9, nc) nu = fill(0.3, nc) - dp = fill(1e6, nc) + # Non-uniform pressure + dp = [Float64(i) * 1e5 for i in 1:nc] result_biot1 = solve_vem_elasticity(setup, E, nu, dp, biot_coefficient = 1.0) result_biot05 = solve_vem_elasticity(setup, E, nu, dp, biot_coefficient = 0.5) + @test norm(result_biot1.displacement) > 0 @test norm(result_biot1.displacement) / norm(result_biot05.displacement) ≈ 2.0 atol=1e-10 end + + @testset "Linear displacement patch test 3D" begin + # Prescribe u(x) = Ax + b on boundary, verify interior follows + mesh = UnstructuredMesh(CartesianMesh((3, 3, 3), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + setup = VEMElasticitySetup(mesh) + E = fill(1e9, nc) + nu = fill(0.3, nc) + + # Test 1: u_x = 0.001*x, u_y = 0, u_z = 0 + bc = zeros(3 * nn) + for (i, pt) in enumerate(mesh.node_points) + bc[3*(i-1)+1] = 0.001 * pt[1] + end + result = solve_vem_elasticity(setup, E, nu, zeros(nc), boundary_displacement = bc) + for (i, pt) in enumerate(mesh.node_points) + @test abs(result.displacement[3*(i-1)+1] - 0.001 * pt[1]) < 1e-14 + @test abs(result.displacement[3*(i-1)+2]) < 1e-14 + @test abs(result.displacement[3*(i-1)+3]) < 1e-14 + end + + # Test 2: u_x = 0.001*y, u_y = 0.002*z, u_z = 0.0005*x (mixed linear) + bc2 = zeros(3 * nn) + for (i, pt) in enumerate(mesh.node_points) + bc2[3*(i-1)+1] = 0.001 * pt[2] + bc2[3*(i-1)+2] = 0.002 * pt[3] + bc2[3*(i-1)+3] = 0.0005 * pt[1] + end + result2 = solve_vem_elasticity(setup, E, nu, zeros(nc), boundary_displacement = bc2) + for (i, pt) in enumerate(mesh.node_points) + @test abs(result2.displacement[3*(i-1)+1] - 0.001 * pt[2]) < 1e-13 + @test abs(result2.displacement[3*(i-1)+2] - 0.002 * pt[3]) < 1e-13 + @test abs(result2.displacement[3*(i-1)+3] - 0.0005 * pt[1]) < 1e-13 + end + end + + @testset "Linear displacement patch test 2D" begin + mesh = UnstructuredMesh(CartesianMesh((4, 4), (1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + setup = VEMElasticitySetup(mesh) + E = fill(1e9, nc) + nu = fill(0.3, nc) + + # u_x = 0.001*x, u_y = 0.0005*y + bc = zeros(2 * nn) + for (i, pt) in enumerate(mesh.node_points) + bc[2*(i-1)+1] = 0.001 * pt[1] + bc[2*(i-1)+2] = 0.0005 * pt[2] + end + result = solve_vem_elasticity(setup, E, nu, zeros(nc), boundary_displacement = bc) + for (i, pt) in enumerate(mesh.node_points) + @test abs(result.displacement[2*(i-1)+1] - 0.001 * pt[1]) < 1e-13 + @test abs(result.displacement[2*(i-1)+2] - 0.0005 * pt[2]) < 1e-13 + end + end end From 50e86cb303909c645317c780bc8044db7be15a8a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 10:13:24 +0000 Subject: [PATCH 5/7] Address code review: minor comment fixes Co-authored-by: moyner <454871+moyner@users.noreply.github.com> --- src/meshes/vem.jl | 2 +- test/vem.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/vem.jl b/src/meshes/vem.jl index f4683c87a..809cd851b 100644 --- a/src/meshes/vem.jl +++ b/src/meshes/vem.jl @@ -1045,7 +1045,7 @@ function _add_face_pressure_contribution!(rhs, node_indices, pts, area_normal_ve end fc ./= nfn - # Unit normal direction (for weighting) + # Unit normal direction nn = norm(area_normal_vec) if nn < eps(Float64) return diff --git a/test/vem.jl b/test/vem.jl index df1b9ffba..d3a45b651 100644 --- a/test/vem.jl +++ b/test/vem.jl @@ -203,7 +203,7 @@ using SparseArrays @test abs(result.displacement[3*(i-1)+3]) < 1e-14 end - # Test 2: u_x = 0.001*y, u_y = 0.002*z, u_z = 0.0005*x (mixed linear) + # Mixed linear displacement: u_x = 0.001*y, u_y = 0.002*z, u_z = 0.0005*x bc2 = zeros(3 * nn) for (i, pt) in enumerate(mesh.node_points) bc2[3*(i-1)+1] = 0.001 * pt[2] From c0d5031a43f7202d1cfe21837defbc733d10f2ac Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 12:20:53 +0000 Subject: [PATCH 6/7] =?UTF-8?q?VEM=20usability:=20reshape=20displacement?= =?UTF-8?q?=20to=20D=C3=97N=20matrix,=20add=20boundary=20node=20helpers?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Reshape displacement output from flat vector to D × num_nodes matrix (e.g. 3×N in 3D), accessible as u[d, node] - Add boundary_nodes(setup) to get list of boundary node indices - Add boundary_nodes_on_side(setup, :xmin) etc. to filter boundary nodes by face of the bounding box - Update tests to use matrix indexing and only set BCs on boundary nodes - All 840 VEM tests + 896 mesh tests pass Co-authored-by: moyner <454871+moyner@users.noreply.github.com> --- src/meshes/vem.jl | 67 +++++++++++++++++++-- test/vem.jl | 150 +++++++++++++++++++++++++++++++++++++--------- 2 files changed, 184 insertions(+), 33 deletions(-) diff --git a/src/meshes/vem.jl b/src/meshes/vem.jl index 809cd851b..491facb70 100644 --- a/src/meshes/vem.jl +++ b/src/meshes/vem.jl @@ -1,4 +1,5 @@ export solve_vem_elasticity, VEMElasticitySetup, assemble_vem_elasticity +export boundary_nodes, boundary_nodes_on_side """ VEMCellData @@ -213,7 +214,7 @@ displacement at the boundaries (or prescribed nonzero boundary displacement). # Returns A named tuple `(displacement, setup, K, rhs)` where: -- `displacement`: nodal displacement vector of length `D * num_nodes`. +- `displacement`: nodal displacement matrix of size `D × num_nodes` (e.g. `3 × N` in 3D). - `setup`: the `VEMElasticitySetup` for reuse. - `K`: the assembled sparse stiffness matrix. - `rhs`: the right-hand side vector. @@ -228,7 +229,7 @@ result = solve_vem_elasticity( fill(0.3, nc), # Poisson's ratio fill(1e6, nc) # Pressure change ) -u = result.displacement +u = result.displacement # 3 × num_nodes matrix ``` """ function solve_vem_elasticity( @@ -250,16 +251,17 @@ end Solve with a precomputed `VEMElasticitySetup` (useful for repeated solves). """ function solve_vem_elasticity( - setup::VEMElasticitySetup, + setup::VEMElasticitySetup{D}, youngs_modulus::AbstractVector, poisson_ratio::AbstractVector, pressure_change::AbstractVector; biot_coefficient::Union{Real, AbstractVector} = 1.0, boundary_displacement::Union{Nothing, AbstractVector} = nothing -) +) where {D} K, rhs = assemble_vem_elasticity(setup, youngs_modulus, poisson_ratio, pressure_change; biot_coefficient = biot_coefficient, boundary_displacement = boundary_displacement) - displacement = K \ rhs + u_vec = K \ rhs + displacement = reshape(u_vec, D, setup.num_nodes) return (displacement = displacement, setup = setup, K = K, rhs = rhs) end @@ -536,6 +538,61 @@ function _find_boundary_nodes(mesh::UnstructuredMesh) return sort(collect(boundary_nodes)) end +""" + boundary_nodes(setup::VEMElasticitySetup) + +Return a sorted vector of global node indices that lie on the boundary of the mesh. + +# Example +```julia +setup = VEMElasticitySetup(mesh) +bnodes = boundary_nodes(setup) +``` +""" +boundary_nodes(setup::VEMElasticitySetup) = setup.boundary_nodes + +""" + boundary_nodes_on_side(setup::VEMElasticitySetup, direction::Symbol; tol=1e-10) + +Return a sorted vector of boundary node indices that lie on a specific side of +the bounding box of the mesh. The `direction` symbol must be one of: +`:xmin`, `:xmax`, `:ymin`, `:ymax` (2D/3D), `:zmin`, `:zmax` (3D only). + +An optional tolerance `tol` is used for the coordinate comparison. + +# Example +```julia +setup = VEMElasticitySetup(mesh) +# Get boundary nodes on the x=0 face +left_nodes = boundary_nodes_on_side(setup, :xmin) +# Get boundary nodes on the y=ymax face +top_nodes = boundary_nodes_on_side(setup, :ymax) +``` +""" +function boundary_nodes_on_side(setup::VEMElasticitySetup{D}, direction::Symbol; tol::Real = 1e-10) where {D} + dir_map = Dict( + :xmin => (1, minimum), + :xmax => (1, maximum), + :ymin => (2, minimum), + :ymax => (2, maximum), + :zmin => (3, minimum), + :zmax => (3, maximum) + ) + + @assert haskey(dir_map, direction) "Unknown direction $direction. Use one of: :xmin, :xmax, :ymin, :ymax" * (D == 3 ? ", :zmin, :zmax" : "") + dim_idx, extremum_fn = dir_map[direction] + @assert dim_idx <= D "Direction $direction not valid for $(D)D mesh" + + pts = setup.mesh.node_points + bnd = setup.boundary_nodes + + # Find the extreme value among boundary nodes in the given dimension + ext_val = extremum_fn(pts[n][dim_idx] for n in bnd) + + # Filter boundary nodes that are at this extreme value + return sort([n for n in bnd if abs(pts[n][dim_idx] - ext_val) <= tol]) +end + # ============================================================================ # VEM stiffness matrix assembly per cell # ============================================================================ diff --git a/test/vem.jl b/test/vem.jl index d3a45b651..cfbb0163d 100644 --- a/test/vem.jl +++ b/test/vem.jl @@ -31,6 +31,76 @@ using SparseArrays end end + @testset "Displacement output shape" begin + @testset "3D output" begin + mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + setup = VEMElasticitySetup(mesh) + dp = [Float64(i) * 1e5 for i in 1:nc] + result = solve_vem_elasticity(setup, fill(1e9, nc), fill(0.3, nc), dp) + @test size(result.displacement) == (3, nn) + end + + @testset "2D output" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3), (1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + setup = VEMElasticitySetup(mesh) + dp = [Float64(i) * 1e5 for i in 1:nc] + result = solve_vem_elasticity(setup, fill(1e9, nc), fill(0.3, nc), dp) + @test size(result.displacement) == (2, nn) + end + end + + @testset "Boundary node helpers" begin + @testset "3D boundary helpers" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3, 3), (1.0, 1.0, 1.0))) + setup = VEMElasticitySetup(mesh) + pts = mesh.node_points + + bnodes = boundary_nodes(setup) + @test length(bnodes) == length(setup.boundary_nodes) + @test bnodes == setup.boundary_nodes + + # Test each side + for (dir, dim_idx, cmp) in [ + (:xmin, 1, 0.0), (:xmax, 1, 1.0), + (:ymin, 2, 0.0), (:ymax, 2, 1.0), + (:zmin, 3, 0.0), (:zmax, 3, 1.0) + ] + side = boundary_nodes_on_side(setup, dir) + @test length(side) > 0 + for n in side + @test n in bnodes + @test abs(pts[n][dim_idx] - cmp) < 1e-10 + end + end + + # xmin and xmax should not overlap + xmin_nodes = boundary_nodes_on_side(setup, :xmin) + xmax_nodes = boundary_nodes_on_side(setup, :xmax) + @test isempty(intersect(xmin_nodes, xmax_nodes)) + end + + @testset "2D boundary helpers" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3), (1.0, 1.0))) + setup = VEMElasticitySetup(mesh) + pts = mesh.node_points + + for (dir, dim_idx, cmp) in [ + (:xmin, 1, 0.0), (:xmax, 1, 1.0), + (:ymin, 2, 0.0), (:ymax, 2, 1.0) + ] + side = boundary_nodes_on_side(setup, dir) + @test length(side) > 0 + for n in side + @test abs(pts[n][dim_idx] - cmp) < 1e-10 + end + end + end + end + @testset "Assembly and solve 3D" begin mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) nc = number_of_cells(mesh) @@ -55,10 +125,9 @@ using SparseArrays @testset "Zero displacement at boundaries" begin result = solve_vem_elasticity(setup, E, nu, dp) u = result.displacement - for node in setup.boundary_nodes + for node in boundary_nodes(setup) for d in 1:3 - dof = 3 * (node - 1) + d - @test abs(u[dof]) < 1e-14 + @test abs(u[d, node]) < 1e-14 end end end @@ -107,10 +176,9 @@ using SparseArrays @testset "Zero displacement at boundaries" begin result = solve_vem_elasticity(setup, E, nu, dp) - for node in setup.boundary_nodes + for node in boundary_nodes(setup) for d in 1:2 - dof = 2 * (node - 1) + d - @test abs(result.displacement[dof]) < 1e-14 + @test abs(result.displacement[d, node]) < 1e-14 end end end @@ -154,6 +222,7 @@ using SparseArrays @testset "Convenience function (mesh input)" begin mesh = UnstructuredMesh(CartesianMesh((2, 2, 2), (1.0, 1.0, 1.0))) nc = number_of_cells(mesh) + nn = length(mesh.node_points) E = fill(1e9, nc) nu = fill(0.3, nc) dp = [Float64(i) * 1e5 for i in 1:nc] @@ -163,7 +232,7 @@ using SparseArrays @test hasproperty(result, :setup) @test hasproperty(result, :K) @test hasproperty(result, :rhs) - @test length(result.displacement) == 3 * length(mesh.node_points) + @test size(result.displacement) == (3, nn) end @testset "Biot coefficient" begin @@ -183,38 +252,41 @@ using SparseArrays end @testset "Linear displacement patch test 3D" begin - # Prescribe u(x) = Ax + b on boundary, verify interior follows + # Prescribe u(x) = Ax + b on boundary nodes only, verify interior follows mesh = UnstructuredMesh(CartesianMesh((3, 3, 3), (1.0, 1.0, 1.0))) nc = number_of_cells(mesh) nn = length(mesh.node_points) setup = VEMElasticitySetup(mesh) E = fill(1e9, nc) nu = fill(0.3, nc) + bnodes = boundary_nodes(setup) - # Test 1: u_x = 0.001*x, u_y = 0, u_z = 0 + # u_x = 0.001*x, u_y = 0, u_z = 0 — set only on boundary nodes bc = zeros(3 * nn) - for (i, pt) in enumerate(mesh.node_points) - bc[3*(i-1)+1] = 0.001 * pt[1] + for n in bnodes + pt = mesh.node_points[n] + bc[3*(n-1)+1] = 0.001 * pt[1] end result = solve_vem_elasticity(setup, E, nu, zeros(nc), boundary_displacement = bc) for (i, pt) in enumerate(mesh.node_points) - @test abs(result.displacement[3*(i-1)+1] - 0.001 * pt[1]) < 1e-14 - @test abs(result.displacement[3*(i-1)+2]) < 1e-14 - @test abs(result.displacement[3*(i-1)+3]) < 1e-14 + @test abs(result.displacement[1, i] - 0.001 * pt[1]) < 1e-14 + @test abs(result.displacement[2, i]) < 1e-14 + @test abs(result.displacement[3, i]) < 1e-14 end - # Mixed linear displacement: u_x = 0.001*y, u_y = 0.002*z, u_z = 0.0005*x + # Mixed linear displacement — set only on boundary nodes bc2 = zeros(3 * nn) - for (i, pt) in enumerate(mesh.node_points) - bc2[3*(i-1)+1] = 0.001 * pt[2] - bc2[3*(i-1)+2] = 0.002 * pt[3] - bc2[3*(i-1)+3] = 0.0005 * pt[1] + for n in bnodes + pt = mesh.node_points[n] + bc2[3*(n-1)+1] = 0.001 * pt[2] + bc2[3*(n-1)+2] = 0.002 * pt[3] + bc2[3*(n-1)+3] = 0.0005 * pt[1] end result2 = solve_vem_elasticity(setup, E, nu, zeros(nc), boundary_displacement = bc2) for (i, pt) in enumerate(mesh.node_points) - @test abs(result2.displacement[3*(i-1)+1] - 0.001 * pt[2]) < 1e-13 - @test abs(result2.displacement[3*(i-1)+2] - 0.002 * pt[3]) < 1e-13 - @test abs(result2.displacement[3*(i-1)+3] - 0.0005 * pt[1]) < 1e-13 + @test abs(result2.displacement[1, i] - 0.001 * pt[2]) < 1e-13 + @test abs(result2.displacement[2, i] - 0.002 * pt[3]) < 1e-13 + @test abs(result2.displacement[3, i] - 0.0005 * pt[1]) < 1e-13 end end @@ -225,17 +297,39 @@ using SparseArrays setup = VEMElasticitySetup(mesh) E = fill(1e9, nc) nu = fill(0.3, nc) + bnodes = boundary_nodes(setup) - # u_x = 0.001*x, u_y = 0.0005*y + # u_x = 0.001*x, u_y = 0.0005*y — set only on boundary nodes bc = zeros(2 * nn) - for (i, pt) in enumerate(mesh.node_points) - bc[2*(i-1)+1] = 0.001 * pt[1] - bc[2*(i-1)+2] = 0.0005 * pt[2] + for n in bnodes + pt = mesh.node_points[n] + bc[2*(n-1)+1] = 0.001 * pt[1] + bc[2*(n-1)+2] = 0.0005 * pt[2] end result = solve_vem_elasticity(setup, E, nu, zeros(nc), boundary_displacement = bc) for (i, pt) in enumerate(mesh.node_points) - @test abs(result.displacement[2*(i-1)+1] - 0.001 * pt[1]) < 1e-13 - @test abs(result.displacement[2*(i-1)+2] - 0.0005 * pt[2]) < 1e-13 + @test abs(result.displacement[1, i] - 0.001 * pt[1]) < 1e-13 + @test abs(result.displacement[2, i] - 0.0005 * pt[2]) < 1e-13 + end + end + + @testset "Boundary side helpers with solve" begin + mesh = UnstructuredMesh(CartesianMesh((3, 3, 3), (1.0, 1.0, 1.0))) + nc = number_of_cells(mesh) + nn = length(mesh.node_points) + setup = VEMElasticitySetup(mesh) + pts = mesh.node_points + + xmin_nodes = boundary_nodes_on_side(setup, :xmin) + @test length(xmin_nodes) > 0 + for n in xmin_nodes + @test abs(pts[n][1]) < 1e-10 + end + + xmax_nodes = boundary_nodes_on_side(setup, :xmax) + @test length(xmax_nodes) > 0 + for n in xmax_nodes + @test abs(pts[n][1] - 1.0) < 1e-10 end end end From 366be5fd374bbdbdf99335988a0ae9094e70984c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 16 Feb 2026 20:43:10 +0100 Subject: [PATCH 7/7] Improve bc treatment --- src/meshes/vem.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/meshes/vem.jl b/src/meshes/vem.jl index 491facb70..eba1e63db 100644 --- a/src/meshes/vem.jl +++ b/src/meshes/vem.jl @@ -78,8 +78,10 @@ right-hand side vector, both with boundary conditions applied. # Keyword arguments - `biot_coefficient`: Biot coefficient (scalar or per-cell vector), default 1.0. -- `boundary_displacement`: If `nothing` (default), zero displacement on all boundary nodes. - If a vector of length `D * num_nodes`, the values at boundary DOFs are used as prescribed displacement. +- `boundary_displacement`: If `nothing` (default), zero displacement on all + boundary nodes. If a vector of length `D * num_nodes`, the values at boundary + DOFs are used as prescribed displacement. Insert NaN if a boundary DOF should + be free (not prescribed). """ function assemble_vem_elasticity( setup::VEMElasticitySetup{D}, @@ -141,10 +143,13 @@ function assemble_vem_elasticity( for d in 1:D dof = D * (node - 1) + d push!(bnd_dof_set, dof) - if boundary_displacement !== nothing - bnd_dof_vals[dof] = boundary_displacement[dof] - else + if isnothing(boundary_displacement) bnd_dof_vals[dof] = 0.0 + else + disp = boundary_displacement[dof] + if isfinite(disp) + bnd_dof_vals[dof] = disp + end end end end