fix: Add deprecated FEMBase methods for backward compatibility

Added deprecated_fembase.jl with legacy methods that tests and user code depend on:
- length(element): Returns number of nodes in element
- size(element): Returns (dim, nnodes)
- getproperty override: Maps element.fields → element.dfields

Bug fix:
- Changed sym == :fields to sym === :fields in getproperty
- Reason: fields.jl overrides == operator, breaking normal Symbol comparisons
- This is a known issue (Code Smell documented in Phase 4 plan)

Result:
- ✅ Element length() works correctly
- ✅ Basic element operations functional
- ⚠️  Test suite still has 44 errors (investigating other API mismatches)

Next: Investigate remaining test failures, likely more API incompatibilities

src/JuliaFEM.jl

@@ -152,6 +152,7 @@ include("assembly/problems.jl")      # Problem types
 include("assembly/assembly.jl")      # Assembly framework
 include("solvers/solvers_base.jl")   # Base solver types
 include("analysis.jl")               # Analysis and AbstractResultsWriter
+include("deprecated_fembase.jl")     # Deprecated/legacy methods from FEMBase (length, size, etc.)
 
 using TimerOutputs
 export @timeit, print_timer

src/deprecated_fembase.jl

@@ -0,0 +1,135 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+"""
+    length(element)
+
+Return the length of basis (number of nodes).
+"""
+function length(element::AbstractElement)
+    return length(element.properties)
+end
+
+"""
+    size(element)
+
+Return the size of basis (dim, nnodes).
+"""
+function size(element::AbstractElement)
+    return size(element.properties)
+end
+
+function getindex(element::AbstractElement, field_name::String)
+    return get_field(element, Symbol(field_name))
+end
+
+function setindex!(element::AbstractElement, data::T, field_name) where T<:AbstractField
+    element.fields[field_name] = data
+end
+
+function setindex!(element::AbstractElement, data::Function, field_name)
+    if hasmethod(data, Tuple{AbstractElement, Vector, Float64})
+        # create enclosure to pass element as argument
+        element.fields[field_name] = field((ip,time) -> data(element,ip,time))
+    else
+        element.fields[field_name] = field(data)
+    end
+end
+
+function setindex!(element::AbstractElement, data, field_name)
+    element.dfields[field_name] = field(data)
+end
+
+function setindex!(fieldset::Dict{Symbol, AbstractField}, field_data, field_name::String)
+    setindex!(fieldset, field_data, Symbol(field_name))
+end
+
+function has_field(element, field_name::String)
+    return has_field(element, Symbol(field_name))
+end
+
+#""" Return a Field object from element.
+#Examples
+#--------
+#>>> element = Element(Seg2, [1, 2])
+#>>> data = Dict(1 => 1.0, 2 => 2.0)
+#>>> update!(element, "my field", data)
+#>>> element("my field")
+#"""
+function (element::Element)(field_name::String)
+    return element[field_name]
+end
+
+function size(element::AbstractElement, dim)
+    return size(element)[dim]
+end
+
+""" Check existence of field. """
+function haskey(element::AbstractElement, field_name)
+    return has_field(element, field_name)
+end
+
+# will be deprecated
+function assemble!(::Assembly, ::Problem{P}, ::AbstractElement, ::Any) where P
+    @warn("One must define assemble! function for problem of type $P. " *
+          "Not doing anything.")
+    return nothing
+end
+
+# will be deprecated
+function assemble!(assembly::Assembly, problem::Problem{P},
+                   elements::Vector{Element}, time) where P
+    @warn("This is default assemble! function. Decreased performance can be " *
+          "expected without preallocation of memory. One should implement " *
+          "`assemble_elements!(problem, assembly, elements, time)` function.")
+    for element in elements
+        assemble!(assembly, problem, element, time)
+    end
+    return nothing
+end
+
+# generally a bad idea to have functions like update! and interpolate, which
+# are not explicitly giving targe?
+
+function update!(element::AbstractElement, field_name, field_data)
+    update_field!(element, field_name, field_data)
+end
+
+function update!(field::AbstractField, data)
+    update_field!(field, data)
+end
+
+function update!(element::AbstractElement, field_name::String, field_data)
+    update_field!(element, Symbol(field_name), field_data)
+end
+
+function update!(elements::Vector{Element}, field_name::String, field_data)
+    update_field!(elements, Symbol(field_name), field_data)
+end
+
+function update!(elements::Vector{Element{T}}, field_name::String, field_data) where T
+    update_field!(convert(Vector{Element}, elements), field_name, field_data)
+end
+
+function interpolate(element::AbstractElement, field_name::String, time)
+    interpolate_field(element, Symbol(field_name), time)
+end
+
+function interpolate(element::AbstractElement, field_name::String, ip, time)
+    interpolate(element, Symbol(field_name), ip, time)
+end
+
+# fields is now dfields
+function Base.getproperty(element::Element, sym::Symbol)
+    if sym === :fields  # Use === instead of == (which is overridden by fields.jl)
+        return getfield(element, :dfields)
+    else
+        return getfield(element, sym)
+    end
+end
+
+# getindex, when someone asks key element.fields["f"] => element.fields[:f]
+function Base.getindex(fieldset::Dict{Symbol,AbstractField}, field_name::String)
+    return getindex(fieldset, Symbol(field_name))
+end
+

src/preprocess.jl

@@ -13,14 +13,14 @@
 =#
 
 mutable struct Mesh
-    nodes :: Dict{Int, Vector{Float64}}
-    node_sets :: Dict{Symbol, Set{Int}}
-    elements :: Dict{Int, Vector{Int}}
-    element_types :: Dict{Int, Symbol}
-    element_codes :: Dict{Int, Symbol}
-    element_sets :: Dict{Symbol, Set{Int}}
-    surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}}
-    surface_types :: Dict{Symbol, Symbol}
+    nodes::Dict{Int,Vector{Float64}}
+    node_sets::Dict{Symbol,Set{Int}}
+    elements::Dict{Int,Vector{Int}}
+    element_types::Dict{Int,Symbol}
+    element_codes::Dict{Int,Symbol}
+    element_sets::Dict{Symbol,Set{Int}}
+    surface_sets::Dict{Symbol,Vector{Tuple{Int,Symbol}}}
+    surface_types::Dict{Symbol,Symbol}
 end
 
 function Mesh()
@@ -68,7 +68,7 @@ end
 
 Add nodes into the mesh.
 """
-function add_nodes!(mesh::Mesh, nodes::Dict{Int, Vector{Float64}})
+function add_nodes!(mesh::Mesh, nodes::Dict{Int,Vector{Float64}})
     for (nid, ncoords) in nodes
         add_node!(mesh, nid, ncoords)
     end
@@ -133,7 +133,7 @@ end
 
 Add elements into the mesh.
 """
-function add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
+function add_elements!(mesh::Mesh, elements::Dict{Int,Tuple{Symbol,Vector{Int}}})
     for (elid, (eltype, elcon)) in elements
         add_element!(mesh, elid, eltype, elcon)
     end
@@ -233,7 +233,7 @@ function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothi
     elements = [create_element(mesh, id) for id in element_ids]
 
     nelements = length(elements)
-    content = Dict{Symbol, Int}()
+    content = Dict{Symbol,Int}()
     for elid in element_ids
         eltype = mesh.element_types[elid]
         content[eltype] = get(content, eltype, 0) + 1
@@ -268,14 +268,14 @@ end
 find npts nearest nodes from the mesh and return their id numbers as a list.
 """
 function find_nearest_nodes(mesh::Mesh, coords::Vector{Float64}, npts::Int=1; node_set=nothing)
-    dist = Dict{Int, Float64}()
+    dist = Dict{Int,Float64}()
     for (nid, c) in mesh.nodes
         if node_set != nothing && !(nid in mesh.node_sets[Symbol(node_set)])
             continue
         end
-        dist[nid] = norm(coords-c)
+        dist[nid] = norm(coords - c)
     end
-    s = sort(collect(dist), by=x->x[2])
+    s = sort(collect(dist), by=x -> x[2])
     nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)]
     node_ids = [n[1] for n in nd]
     return node_ids
@@ -298,7 +298,7 @@ mapping :: Dict{Symbol, Vector{Int}}
     e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
 
 """
-function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{Int}})
+function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol,Vector{Int}})
     for (elid, eltype) in mesh.element_types
         haskey(mapping, eltype) || continue
         new_order = mapping[eltype]
@@ -330,7 +330,7 @@ It is safe to assemble elements with the same color in parallel
 """
 function create_coloring(mesh::Mesh)
     # Contains the elements that each node contain
-    cell_containing_node = Dict{Int, Set{Int}}()
+    cell_containing_node = Dict{Int,Set{Int}}()
     for (cellid, nodes) in mesh.elements
         for v in nodes
             if !haskey(cell_containing_node, v)
@@ -355,7 +355,7 @@ function create_coloring(mesh::Mesh)
 
     incidence_matrix = sparse(I, J, V)
     # cell -> color of cell
-    cell_colors = Dict{Int, Int}()
+    cell_colors = Dict{Int,Int}()
     # color -> list of cells
     final_colors = Set{Int}[]
     occupied_colors = Set{Int}()

src/solvers_modal.jl

@@ -4,24 +4,24 @@
 using SparseArrays, Arpack
 
 mutable struct Modal <: AbstractSolver
-    time :: Float64
-    geometric_stiffness :: Bool
-    eigvals :: Vector
-    eigvecs :: Matrix
-    nev :: Int
-    which :: Symbol
-    bc_invertible :: Bool
-    P :: Vector{SparseMatrixCSC}
-    symmetric :: Bool
-    empty_assemblies_before_solution :: Bool
-    dense :: Bool
-    info_matrices :: Bool
-    sigma :: Float64
+    time::Float64
+    geometric_stiffness::Bool
+    eigvals::Vector
+    eigvecs::Matrix
+    nev::Int
+    which::Symbol
+    bc_invertible::Bool
+    P::Vector{SparseMatrixCSC}
+    symmetric::Bool
+    empty_assemblies_before_solution::Bool
+    dense::Bool
+    info_matrices::Bool
+    sigma::Float64
 end
 
 function Modal(nev=10, which=:SM)
-    solver = Modal(0.0, false, [], Matrix{Float64}(undef,0,0), nev, which,
-                   false, [], true, true, false, false, 0.0)
+    solver = Modal(0.0, false, [], Matrix{Float64}(undef, 0, 0), nev, which,
+        false, [], true, true, false, false, 0.0)
 end
 
 """
@@ -66,10 +66,10 @@ function eliminate_boundary_conditions!(problem::P, K, M, f) where {P}
     isdiag(C1) || error("Cannot eliminate boundary condition $P: C is not diagonal")
     @info("Eliminating boundary condition $(problem.name) from global system.")
     fixed_dofs = get_nonzero_rows(C1)
-    K[fixed_dofs,:] .= 0.0
-    K[:,fixed_dofs] .= 0.0
-    M[fixed_dofs,:] .= 0.0
-    M[:,fixed_dofs] .= 0.0
+    K[fixed_dofs, :] .= 0.0
+    K[:, fixed_dofs] .= 0.0
+    M[fixed_dofs, :] .= 0.0
+    M[:, fixed_dofs] .= 0.0
     dropzeros!(K)
     dropzeros!(M)
     return nothing
@@ -118,8 +118,8 @@ function run!(solver::Solver{Modal})
 
     for P in properties.P
         @info("Using P to make transformation K_red = P'*K*P and M_red = P'*M*P")
-        K[:,:] .= P'*K*P
-        M[:,:] .= P'*M*P
+        K[:, :] .= P' * K * P
+        M[:, :] .= P' * M * P
     end
 
     for problem in get_problems(solver)
@@ -136,8 +136,8 @@ function run!(solver::Solver{Modal})
     SparseArrays.droptol!(K, 1.0e-9)
     SparseArrays.droptol!(M, 1.0e-9)
     nz = get_nonzero_rows(K)
-    K = K[nz,nz]
-    M = M[nz,nz]
+    K = K[nz, nz]
+    M = M[nz, nz]
 
     sigma = 0.0
     if properties.sigma != 0.0
@@ -172,12 +172,12 @@ function run!(solver::Solver{Modal})
 
     try
         @timeit "solve eigenvalue problem using `eigs`" begin
-            om2, X = eigs(K + sigma*I, M; nev=props.nev, which=props.which)
+            om2, X = eigs(K + sigma * I, M; nev=props.nev, which=props.which)
         end
         passed = true
     catch
         @info("Failed to calculate eigenvalues for problem.",
-              issymmetric(K), issymmetric(M), isposdef(K), isposdef(M))
+            issymmetric(K), issymmetric(M), isposdef(K), isposdef(M))
         if !isapprox(properties.sigma, 0.0)
             @info("Stiffness matrix is not positive definite and Cholesky " *
                   "factorization is failing. Model is not supported enough " *
@@ -195,7 +195,7 @@ function run!(solver::Solver{Modal})
               "positive definite and Cholesky factorization is failing. Trying " *
               "again by adjusting problem.properties.sigma to $sigma.")
         try
-            om2, X = eigs(K + sigma*I, M; nev=props.nev, which=props.which)
+            om2, X = eigs(K + sigma * I, M; nev=props.nev, which=props.which)
             passed = true
         catch
             @info("Failed to calculate eigenvalues with sigma value $sigma. " *
@@ -209,15 +209,15 @@ function run!(solver::Solver{Modal})
     props.eigvals = om2
     neigvals = length(om2)
     props.eigvecs = zeros(ndofs, neigvals)
-    for i=1:neigvals
-        props.eigvecs[nz,i] = X[:,i]
+    for i = 1:neigvals
+        props.eigvecs[nz, i] = X[:, i]
         for problem in get_boundary_problems(solver)
             isa(problem, Problem{Mortar}) || continue
             s, m, P = calc_projection(problem)
             # FIXME: store projection to boundary problem, i.e.
             # update!(problem, "master-slave projection", time => P)
             # us = P*um
-            props.eigvecs[s,i] = P*props.eigvecs[m,i]
+            props.eigvecs[s, i] = P * props.eigvecs[m, i]
         end
     end
 
@@ -245,13 +245,13 @@ function update_xdmf!(solver::Solver{Modal})
 
     @timeit "fetch geometry" X_ = solver("geometry", solver.properties.time)
     node_ids = keys(X_)
-    @timeit "create node permutation" P = Dict(j=>i for (i, j) in enumerate(node_ids))
+    @timeit "create node permutation" P = Dict(j => i for (i, j) in enumerate(node_ids))
     nnodes = length(X_)
-    ndofs = round(Int, size(solver.properties.eigvecs, 1)/nnodes)
+    ndofs = round(Int, size(solver.properties.eigvecs, 1) / nnodes)
     ndim = length(X_[first(node_ids)])
     @info("Number of nodes: $nnodes. ",
-          "Number of dofs/node: $ndofs. ",
-          "Dimension of geometry: $ndim.")
+        "Number of dofs/node: $ndofs. ",
+        "Dimension of geometry: $ndim.")
     @timeit "create ncoords array" begin
         X = zeros(ndim, nnodes)
         for j in node_ids
@@ -263,7 +263,7 @@ function update_xdmf!(solver::Solver{Modal})
     element_types = unique(map(get_element_type, all_elements))
     nelements = length(all_elements)
 
-    elcon_arrays  = Dict()
+    elcon_arrays = Dict()
     @timeit "create topology arrays" for element_type in element_types
         elements = collect(filter_by_element_type(element_type, all_elements))
         nelements = length(elements)
@@ -271,7 +271,7 @@ function update_xdmf!(solver::Solver{Modal})
         element_conn = zeros(Int, eldim, nelements)
         for (i, element) in enumerate(elements)
             for (j, conn) in enumerate(get_connectivity(element))
-                element_conn[j,i] = P[conn]-1
+                element_conn[j, i] = P[conn] - 1
             end
         end
         elcon_arrays[element_type] = element_conn
@@ -305,7 +305,7 @@ function update_xdmf!(solver::Solver{Modal})
             @warn("negative real eigenvalue found, om2=$eigval, setting to zero.")
             eigval = 0.0
         end
-        freq = sqrt(eigval)/(2.0*pi)
+        freq = sqrt(eigval) / (2.0 * pi)
         path = "/Results/Natural Frequency Analysis/$unknown_field_name/Mode $j"
         @info("Creating frequency frame f=$(round(freq; digits=3)), path=$path")
 
@@ -340,11 +340,11 @@ function update_xdmf!(solver::Solver{Modal})
         end
 
         @timeit "store eigenmode" begin
-            mode_ = reshape(solver.properties.eigvecs[:,j], ndofs, nnodes)
+            mode_ = reshape(solver.properties.eigvecs[:, j], ndofs, nnodes)
             @timeit "create mode array" begin
                 mode = zeros(ndofs, nnodes)
                 for nid in node_ids
-                    mode[:,P[nid]] = mode_[:,nid]
+                    mode[:, P[nid]] = mode_[:, nid]
                 end
             end
             if ndofs == 1