From aab8b7d6ce72adc0e337311ed56411b49dcc5def Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Sun, 9 Nov 2025 18:42:56 +0200
Subject: [PATCH] feat(test): First test rewritten for immutable elements
 (test_elasticity_1d)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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
---
 src/JuliaFEM.jl            | 24 ++++++++++++++++++++++++
 src/elements/elements.jl   | 17 ++++++++++++-----
 test/test_elasticity_1d.jl | 26 +++++++++++++++++++-------
 3 files changed, 55 insertions(+), 12 deletions(-)

diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 5504fd7..2a60d64 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -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
diff --git a/src/elements/elements.jl b/src/elements/elements.jl
index c4ab03a..5506113 100644
--- a/src/elements/elements.jl
+++ b/src/elements/elements.jl
@@ -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
diff --git a/test/test_elasticity_1d.jl b/test/test_elasticity_1d.jl
index 02df34d..8e57ba6 100644
--- a/test/test_elasticity_1d.jl
+++ b/test/test_elasticity_1d.jl
@@ -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