feat(test): First test rewritten for immutable elements (test_elasticity_1d)

Rewrote test_elasticity_1d.jl to follow immutable element pattern.
This is the first fully working test with the new architecture!

Changes:
1. test/test_elasticity_1d.jl:
   - Convert Dict node data to element-local tuple format
   - Wrap data in DVTI field objects (Discrete, Variable, Time-Invariant)
   - Create element with fields at construction: Element(Seg2, conn; fields=(...))
   - Fix Jacobian shape expectation (3×1 not 1×3 for 1D in 3D)

2. src/JuliaFEM.jl:
   - Add minimal jacobian() function for AbstractBasis (non-parametric)
   - Handles embedding (1D element in 3D space) correctly
   - Returns Matrix instead of Tensor for flexibility

3. src/elements/elements.jl:
   - Fix Jacobian computation to handle both Tuple and IntegrationPoint
   - Fix detJ calculation logic for embedded elements (check m not size(JT,2))
   - Correctly handle 1D elements: detJ = ||∂X/∂ξ||

Result: test_elasticity_1d.jl passes! ✓

This validates the immutable architecture:
- Element created with fields at construction
- No mutation needed during test
- Field system integration working (DVTI fields)
- Jacobian computation working for embedded elements
This commit is contained in:
Jukka Aho
2025-11-09 18:42:56 +02:00
parent 41e09b2c92
commit aab8b7d6ce
3 changed files with 55 additions and 12 deletions
+24
View File
@@ -222,6 +222,30 @@ include("basis/nurbs.jl")
# TODO: Rewrite for new AbstractBasis (non-parametric)
# include("basis/math.jl") # Uses AbstractBasis{dim} throughout (jacobian, grad, interpolate, etc.)
# TODO: Rewrite math functions for new AbstractBasis
# TEMPORARY: Define minimal jacobian function for testing
function jacobian(B::AbstractBasis, X::Vector{<:Vec}, xi::Vec)
dB = eval_dbasis!(B, xi)
@assert length(X) == length(dB)
# Compute J = dX/dξ: rows are physical dims, columns are parametric dims
# J[i,j] = ∂X_i/∂ξ_j = sum_k X_k[i] * dN_k/dξ_j
dim_physical = length(first(X))
dim_parametric = length(xi)
# Build Jacobian matrix manually for embedding case (e.g., 1D element in 3D space)
# Result is a dim_physical × dim_parametric matrix
J_data = zeros(dim_physical, dim_parametric)
@inbounds for k in 1:length(X)
for i in 1:dim_physical
for j in 1:dim_parametric
J_data[i,j] += X[k][i] * dB[k][j]
end
end
end
# Convert to Tensor (note: Tensor{2,N} is N×N, but we need dim_physical×dim_parametric)
# For now, return as Matrix
return J_data
end
# Consolidate FEMBase.jl into src/ (Phase 1 continued)
# Order matters: fields → types → sparse → elements → integrate → problems → assembly
+12 -5
View File
@@ -641,9 +641,13 @@ function (element::Element)(ip, time, ::Type{Val{:Jacobian}})
X_dict = element("geometry", time)
# Convert to Vector{Vec} for Tensors.jl compatibility
X = [Vec(x...) for x in X_dict]
# Convert ip.coords (Tuple) to Vec
xi = Vec(ip.coords)
J = jacobian(element.properties, X, xi)
# Convert ip to Vec - handle both Tuple and IntegrationPoint
if isa(ip, Tuple)
xi = Vec(ip)
else
xi = Vec(ip.coords)
end
J = jacobian(element.basis, X, xi)
return J
end
@@ -653,10 +657,13 @@ function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}})
if n == m # volume element
return det(J)
end
# For embedded elements (1D in 2D/3D, 2D in 3D):
# detJ = || ∂X/∂ξ || for 1D elements
# detJ = || ∂X/∂ξ₁ × ∂X/∂ξ₂ || for 2D elements
JT = transpose(J)
if size(JT, 2) == 1 # boundary of 2d problem, || ∂X/∂ξ ||
if m == 1 # 1D element (boundary of 2D or 3D), J is n×1, JT is 1×n
return norm(JT)
else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
else # 2D element (manifold on 3D problem), J is 3×2, JT is 2×3
return norm(cross(JT[:, 1], JT[:, 2]))
end
end
+19 -7
View File
@@ -5,15 +5,27 @@ using JuliaFEM, Test
# 1d strain
X = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
u = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
element = Element(Seg2, (1, 2))
update!(element, "geometry", X)
update!(element, "displacement", u)
# Global node data (Dict format for backward compatibility in tests)
X_global = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
u_global = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
# Convert to element-local format (extract data for element nodes)
connectivity = (1, 2)
X = tuple([X_global[i] for i in connectivity]...)
u = tuple([u_global[i] for i in connectivity]...)
# Wrap in field objects (DVTI = Discrete, Variable, Time-Invariant)
X_field = JuliaFEM.DVTI(X)
u_field = JuliaFEM.DVTI(u)
# Create element with fields at construction (immutable pattern)
element = Element(Seg2, connectivity; fields=(geometry=X_field, displacement=u_field))
xi, time = (0.0,), 0.0
detJ = element(xi, time, Val{:detJ})
J = element(xi, time, Val{:Jacobian})
# gradu = element("displacement", xi, time, Val{:Grad})
@debug("1d seg2 info", xi ,time, detJ, J)
@debug("1d seg2 info", xi, time, detJ, J)
@test isapprox(detJ, sqrt(3)/2)
@test isapprox(J, [0.5 0.5 0.5])
# Jacobian is 3×1 (physical_dim × parametric_dim) for 1D element in 3D
@test isapprox(J, [0.5; 0.5; 0.5]) # column vector