feat: Consolidate FEMBase.jl into JuliaFEM (Phase 1 complete)

MAJOR MILESTONE: FEMBase + FEMBasis fully consolidated, JuliaFEM loads! Consolidated files: - src/elements/ (3 files): elements.jl, elements_lagrange.jl, integrate.jl - src/fields/ (1 file): fields.jl (DCTI, DVTI, DCTV, DVTV, etc.) - src/sparse/ (1 file): sparse.jl (SparseMatrixCOO, SparseVectorCOO) - src/assembly/ (2 files): problems.jl, assembly.jl - src/solvers/ (1 file): solvers_base.jl - src/analysis.jl, src/core_types.jl (Node, IP, IntegrationPoint) Changes to JuliaFEM.jl: - Added dependencies: Tensors, Calculus - Removed @reexport using FEMBase (now consolidated) - Added 20+ include statements for consolidated files - Include order: fields → core_types → fembase_compat → sparse → elements Compatibility layer: - Created fembase_compat.jl: Minimal FEMBase submodule for vendor packages - Temporarily disabled vendor-specific Mortar2D functions in solvers_modal.jl Bug fixes: - Changed i == 1 → isequal(i, 1) in integrate.jl (== operator overridden by fields) - Resolved all FEMBasis. namespace references throughout codebase Result: - ✅ JuliaFEM loads successfully on Julia 1.12.1 - ✅ 134 exported symbols (was 171 with separate FEMBase) - ✅ Core types accessible: Seg2, Quad4, Problem, AbstractProblem, etc. - ⚠️ Vendor packages show FEMBase cache warnings (expected, harmless) TODO: - Re-enable Mortar2D functions after vendor consolidation - Field system == operator override needs redesign (Phase 4) - Continue Phase 2: Consolidate remaining vendor packages
2026-10-03 22:57:57 +00:00 · 2025-11-08 09:39:16 +02:00
parent d3fc55f13e
commit 73ec910589
31 changed files with 2228 additions and 55 deletions
@@ -0,0 +1,527 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+"""
+    AbstractFieldSet{N<:Int}
+
+Abstract supertype for all field sets, where `N` is the length of the discrete
+fields (typically is the number of the nodes in element).
+"""
+abstract type AbstractFieldSet{N} end
+
+"""
+    EmptyFieldSet{N} <: AbstractFieldSet{N}
+
+Empty field set used as a default for all elements.
+"""
+struct EmptyFieldSet{N} <: AbstractFieldSet{N}
+end
+
+const DefaultFieldSet = EmptyFieldSet
+
+"""
+    AbstractElement{M<:AbstractFieldSet, B<:AbstractBasis}
+
+Abstract supertype for all elements.
+"""
+abstract type AbstractElement{M<:AbstractFieldSet, B<:AbstractBasis} end
+
+mutable struct Element{M,B} <: AbstractElement{M,B}
+    id :: Int
+    connectivity :: Vector{Int}
+    integration_points :: Vector{IP}
+    dfields :: Dict{Symbol, AbstractField}
+    sfields :: M
+    properties :: B
+end
+
+"""
+    Element(topology, connectivity)
+
+Construct a new element where `topology` is the topological type of the element
+and connectivity contains node numbers where element is connected.
+
+# Topological types
+
+## 1d elements
+- `Seg2`
+- `Seg3`
+
+## 2d elements
+- `Tri3`
+- `Tri6`
+- `Tri7`
+- `Quad4`
+- `Quad8`
+- `Quad9`
+
+## 3d elements
+- `Tet4`
+- `Tet10`
+- `Hex8`
+- `Hex20`
+- `Hex27`
+- `Pyr5`
+- `Wedge6`
+- `Wedge15`
+
+# Examples
+
+```julia
+element = Element(Tri3, (1, 2, 3))
+```
+"""
+function Element(::Type{T}, connectivity::NTuple{N, Int}) where {N, T<:AbstractBasis}
+    return Element(T, DefaultFieldSet, connectivity)
+end
+
+function Element(::Type{T}, ::Type{M}, connectivity::NTuple{N, Int}) where {N, M<:AbstractFieldSet, T<:AbstractBasis}
+    element_id = -1
+    topology = T()
+    integration_points = Point{IntegrationPoint}[]
+    dfields = Dict{Symbol,AbstractField}()
+    sfields = M{N}()
+    element = Element(element_id, collect(connectivity), integration_points,
+                      dfields, sfields, topology)
+    return element
+end
+
+function Element(::Type{T}, connectivity::Vector{Int}) where T<:AbstractBasis
+    return Element(T, (connectivity...,))
+end
+
+function get_element_id(element::AbstractElement)
+    return element.id
+end
+
+function get_element_type(::AbstractElement{M,T}) where {M,T}
+    return T
+end
+
+function is_element_type(::AbstractElement{M,T}, element_type) where {M,T}
+    return T === element_type
+end
+
+function filter_by_element_type(element_type, elements)
+    return Iterators.filter(element -> is_element_type(element, element_type), elements)
+end
+
+function get_connectivity(element::AbstractElement)
+    return element.connectivity
+end
+
+"""
+    group_by_element_type(elements)
+
+Given a vector of elements, group elements by element type to several vectors.
+Returns a dictionary, where key is the element type and value is a vector
+containing all elements of type `element_type`.
+"""
+function group_by_element_type(elements)
+    eltypes = map(T -> typeof(T), elements)
+    elgroups = Dict(T => T[] for T in eltypes)
+    for element in elements
+        T = typeof(element)
+        push!(elgroups[T], element)
+    end
+    return elgroups
+end
+
+### dfields - dynamically defined fields
+
+# This is the "old" field system, where fields are defined to dictionary.
+# It is known that this approach is having a performance issue caused by
+# type instability.
+
+function has_dfield(element, field_name)
+    return haskey(element.dfields, field_name)
+end
+
+function get_dfield(element, field_name)
+    return getindex(element.dfields, field_name)
+end
+
+function create_dfield!(element, field_name, field_::AbstractField)
+    T = typeof(field_)
+    if has_dfield(element, field_name)
+        @debug("Replacing the content of a field $field_name with a new field of type $T.")
+    else
+        @debug("Creating a new dfield $field_name of type $T")
+    end
+    element.dfields[field_name] = field_
+    return
+end
+
+function create_dfield!(element, field_name, field_data)
+    create_dfield!(element, field_name, field(field_data))
+end
+
+function update_dfield!(element, field_name, field_data)
+    if has_dfield(element, field_name)
+        field = get_dfield(element, field_name)
+        @debug("Update $field_name with data $field_data")
+        update_field!(field, field_data)
+    else
+        create_dfield!(element, field_name, field_data)
+    end
+end
+
+# A helper function to pick element data from dictionary
+function pick_data_(element, field_data)
+    connectivity = get_connectivity(element)
+    N = length(connectivity)
+    picked_data = ntuple(i -> getindex(field_data, connectivity[i]), N)
+    return picked_data
+end
+
+function update_dfield!(element, field_name, (time, field_data)::Pair{Float64, Dict{Int,V}}) where V
+    update_dfield!(element, field_name, time => pick_data_(element, field_data))
+end
+
+function update_dfield!(element, field_name, field_data::Dict{Int,V}) where V
+    update_dfield!(element, field_name, pick_data_(element, field_data))
+end
+
+function update_dfield!(element, field_name, field_data::Function)
+    if hasmethod(field_data, Tuple{Element, Any, Any})
+        element.dfields[field_name] = field((ip, time) -> field_data(element, ip, time))
+    else
+        element.dfields[field_name] = field(field_data)
+    end
+end
+
+function interpolate_dfield(element, field_name, time)
+    field = get_dfield(element, field_name)
+    return interpolate(field, time)
+end
+
+### sfields statically defined fields
+
+# A new-style field system, where fields are defined in sfields <: AbstractFieldSet
+# during the initialization of element.
+
+function has_sfield(element, field_name)
+    return isdefined(element.sfields, field_name)
+end
+
+function get_sfield(element, field_name)
+    return getfield(element.sfields, field_name)
+end
+
+function update_sfield!(element, field_name, field_data)
+    field = get_sfield(element, field_name)
+    update!(field, field_data)
+end
+
+function interpolate_sfield(element, field_name, time)
+    field = get_sfield(element, field_name)
+    return interpolate(field, time)
+end
+
+### dfield & sfield -- common routines
+
+function has_field(element, field_name)
+    return has_sfield(element, field_name) || has_dfield(element, field_name)
+end
+
+function get_field(element, field_name)
+    if has_sfield(element, field_name)
+        return get_sfield(element, field_name)
+    else
+        return get_dfield(element, field_name)
+    end
+end
+
+function update_field!(element, field_name, field_data)
+    if has_sfield(element, field_name)
+        update_sfield!(element, field_name, field_data)
+    else
+        update_dfield!(element, field_name, field_data)
+    end
+end
+
+function interpolate_field(element, field_name::Symbol, time)
+    if has_sfield(element, field_name)
+        return interpolate_sfield(element, field_name, time)
+    elseif has_dfield(element, field_name)
+        return interpolate_dfield(element, field_name, time)
+    else
+        error("Cannot interpolate from field $field_name: no such field.")
+    end
+end
+
+function interpolate(element::AbstractElement, field_name, time)
+    return interpolate_field(element, field_name, time)
+end
+
+function update_field!(elements::Vector{Element}, field_name, field_data)
+    for element in elements
+        update_field!(element, field_name, field_data)
+    end
+end
+
+# Update fields when given a dictionary or time => dictionary:
+# pick data from dictionary diven by the connectivity information of element
+#=
+function update_field!(element::AbstractElement, field::F,
+                       data::Dict{T,V}) where {F<:DVTI,T,V}
+    connectivity = get_connectivity(element)
+    N = length(connectivity)
+    picked_data = ntuple(i -> data[connectivity[i]], N)
+    update_field!(field, picked_data)
+end
+
+function update_field!(element::AbstractElement, field::F,
+                       ddata::Pair{Float64, Dict{T,V}}) where {F<:DVTV,T,V}
+    time, data = ddata
+    connectivity = get_connectivity(element)
+    N = length(connectivity)
+    picked_data = ntuple(i -> data[connectivity[i]], N)
+    update_field!(field, time => picked_data)
+end
+=#
+
+"""
+    interpolate(element, field_name, time)
+
+Interpolate field `field_name` from element at given `time`.
+
+# Example
+```
+element = Element(Seg2, [1, 2])
+data1 = Dict(1 => 1.0, 2 => 2.0)
+data2 = Dict(1 => 2.0, 2 => 3.0)
+update!(element, "my field", 0.0 => data1)
+update!(element, "my field", 1.0 => data2)
+interpolate(element, "my field", 0.5)
+
+# output
+
+(1.5, 2.5)
+
+```
+"""
+function interpolate(element::AbstractElement, field_name::String, time::Float64)
+    field = element[field_name]
+    result = interpolate(field, time)
+    if isa(result, Dict)
+        connectivity = get_connectivity(element)
+        return tuple((result[i] for i in connectivity)...)
+    else
+        return result
+    end
+end
+
+function info_update_field(elements, field_name, data)
+    nelements = length(elements)
+    @info("Updating field `$field_name` for $nelements elements.")
+end
+
+function info_update_field(elements, field_name, data::Float64)
+    nelements = length(elements)
+    @info("Updating field `$field_name` => $data for $nelements elements.")
+end
+
+"""
+    update!(elements, field_name, data)
+
+Given a list of elements, field name and data, update field to elements. Data
+is passed directly to the `field`-function.
+
+# Examples
+
+Create two elements with topology `Seg2`, one is connecting to nodes (1, 2) and
+the other is connecting to (2, 3). Some examples of updating fields:
+
+```julia
+elements = [Element(Seg2, [1, 2]), Element(Seg2, [2, 3])]
+X = Dict(1 => 0.0, 2 => 1.0, 3 => 2.0)
+u = Dict(1 => 0.0, 2 => 0.0, 3 => 0.0)
+update!(elements, "geometry", X)
+update!(elements, "displacement", 0.0 => u)
+update!(elements, "youngs modulus", 210.0e9)
+update!(elements, "time-dependent force", 0.0 => 0.0)
+update!(elements, "time-dependent force", 1.0 => 100.0)
+```
+
+When using dictionaries in definition of fields, key of dictionary corresponds
+to node id, that is, updating field `geometry` in the example above is updating
+values `(0.0, 1.0)` for the first elements and values `(1.0, 2.0)` to the second
+element. For time dependent field, syntax `time => data` is used. If field is
+initialized without time-dependency, it cannot be changed to be time-dependent
+afterwards. If unsure, it's better to initialize field with time dependency.
+
+"""
+function update!(elements, field_name, data)
+    info_update_field(elements, field_name, data)
+    for element in elements
+        update!(element, field_name, data)
+    end
+end
+
+
+## Interpolate fields in spatial direction
+
+const ConstantField = Union{DCTI, DCTV}
+const VariableFields = Union{DVTV, DVTI}
+const DictionaryFields = Union{DVTVd, DVTId}
+
+function interpolate_field(::AbstractElement, field::ConstantField, ip, time)
+    return interpolate_field(field, time)
+end
+
+function interpolate_field(element::AbstractElement, field::VariableFields, ip, time)
+    data = interpolate_field(field, time)
+    basis = get_basis(element, ip, time)
+    N = length(basis)
+    return sum(data[i]*basis[i] for i=1:N)
+end
+
+function interpolate_field(element::AbstractElement, field::DictionaryFields, ip, time)
+    data = interpolate_field(field, time)
+    basis = element(ip, time)
+    N = length(element)
+    c = get_connectivity(element)
+    return sum(data[c[i]]*basis[i] for i=1:N)
+end
+
+function interpolate_field(::AbstractElement, field::CVTV, ip, time)
+    return field(ip, time)
+end
+
+function interpolate(element::AbstractElement, field_name, ip, time)
+    field = get_field(element, field_name)
+    interpolate_field(element, field, ip, time)
+end
+
+
+## Other stuff
+
+function get_basis(element::AbstractElement{M,B}, ip, ::Any) where {M,B}
+    T = typeof(first(ip))
+    N = zeros(T, 1, length(element))
+    eval_basis!(B, N, tuple(ip...))
+    return N
+end
+
+function get_dbasis(element::AbstractElement{M,B}, ip, ::Any) where {M,B}
+    T = typeof(first(ip))
+    dN = zeros(T, size(element)...)
+    eval_dbasis!(B, dN, tuple(ip...))
+    return dN
+end
+
+function (element::Element)(ip, time::Float64=0.0)
+    return get_basis(element, ip, time)
+end
+
+#"""
+#Examples
+#julia> el = Element(Quad4, [1, 2, 3, 4]);
+#julia> el([0.0, 0.0], 0.0, 1)
+#1x4 Array{Float64,2}:
+# 0.25  0.25  0.25  0.25
+#julia> el([0.0, 0.0], 0.0, 2)
+#2x8 Array{Float64,2}:
+# 0.25  0.0   0.25  0.0   0.25  0.0   0.25  0.0
+# 0.0   0.25  0.0   0.25  0.0   0.25  0.0   0.25
+#"""
+function (element::Element)(ip, time::Float64, dim::Int)
+    dim == 1 && return get_basis(element, ip, time)
+    Ni = vec(get_basis(element, ip, time))
+    N = zeros(dim, length(element)*dim)
+    for i=1:dim
+        N[i,i:dim:end] += Ni
+    end
+    return N
+end
+
+function (element::Element)(ip, time, ::Type{Val{:Jacobian}})
+    X = element("geometry", time)
+    J = jacobian(element.properties, X, ip)
+    return J
+end
+
+function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}})
+    J = element(ip, time, Val{:Jacobian})
+    n, m = size(J)
+    if n == m  # volume element
+        return det(J)
+    end
+    JT = transpose(J)
+    if size(JT, 2) == 1  # boundary of 2d problem, || ∂X/∂ξ ||
+        return norm(JT)
+    else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
+        return norm(cross(JT[:,1], JT[:,2]))
+    end
+end
+
+function (element::Element)(ip, time::Float64, ::Type{Val{:Grad}})
+    J = element(ip, time, Val{:Jacobian})
+    return inv(J)*get_dbasis(element, ip, time)
+end
+
+function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:Grad}})
+    X = element("geometry", time)
+    u = element(field_name, time)
+    return grad(element.properties, u, X, ip)
+end
+
+
+function get_integration_points(element::AbstractElement{E}) where E
+    # first time initialize default integration points
+    if length(element.integration_points) == 0
+        ips = get_integration_points(element.properties)
+        element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
+    end
+    return element.integration_points
+end
+
+""" This is a special case, temporarily change order
+of integration scheme mainly for mass matrix.
+"""
+function get_integration_points(element::AbstractElement{E}, change_order::Int) where E
+    ips = get_integration_points(element.properties, Val{change_order})
+    return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
+end
+
+""" Find inverse isoparametric mapping of element. """
+function get_local_coordinates(element::AbstractElement, X::Vector, time::Float64; max_iterations=10, tolerance=1.0e-6)
+    haskey(element, "geometry") || error("element geometry not defined, cannot calculate inverse isoparametric mapping")
+    dim = size(element, 1)
+    dim == length(X) || error("manifolds not supported.")
+    xi = zeros(dim)
+    dX = element("geometry", xi, time) - X
+    for i=1:max_iterations
+        J = element(xi, time, Val{:Jacobian})'
+        xi -= J \ dX
+        dX = element("geometry", xi, time) - X
+        norm(dX) < tolerance && return xi
+    end
+    debug("get_local_coordinates", X, dX, xi)
+    error("Unable to find inverse isoparametric mapping for element $element for X = $X")
+end
+
+""" Test is X inside element. """
+function inside(element::AbstractElement{M,B}, X, time) where {M,B}
+    xi = get_local_coordinates(element, X, time)
+    return inside(B, xi)
+end
+
+## Convenience functions
+
+# element("displacement", 0.0)
+function (element::Element)(field_name::String, time::Float64)
+    return interpolate(element, field_name, time)
+end
+
+# element("displacement", (0.0, 0.0), 0.0)
+function (element::Element)(field_name::String, ip, time::Float64)
+    return interpolate(element, field_name, ip, time)
+end
+
+function element_info!(bi::BasisInfo{T}, element::AbstractElement{M,T}, ip, time) where {M,T}
+    X = interpolate(element, "geometry", time)
+    eval_basis!(bi, X, ip)
+    return bi.J, bi.detJ, bi.N, bi.grad
+end
@@ -0,0 +1,53 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+struct Poi1 <: AbstractBasis{0} end
+
+function get_basis(::E, ::Any, ::Any) where E<:AbstractElement{M,Poi1} where M
+    return [1]
+end
+
+function get_dbasis(::E, ::Any, ::Any) where E<:AbstractElement{M,Poi1} where M
+    return [0]
+end
+
+function (::Element{M,Poi1})(::Any, ::Float64, ::Type{Val{:detJ}}) where M
+    return 1.0
+end
+
+function get_integration_order(::Poi1)
+    return 1
+end
+
+function get_integration_points(::Poi1, ::Int)
+    return [(1.0, (0.0,))]
+end
+
+function size(::Type{Poi1})
+    return (0, 1)
+end
+
+function length(::Type{Poi1})
+    return 1
+end
+
+function get_reference_element_coordinates(::Type{Poi1})
+    Vector{Float64}[[0.0]]
+end
+
+function inside(::Union{Type{Seg2},Type{Seg3},Type{Quad4},
+        Type{Quad8},Type{Quad9},Type{Pyr5},
+        Type{Hex8},Type{Hex20},
+        Type{Hex27}}, xi)
+    return all(-1.0 .<= xi .<= 1.0)
+end
+
+function inside(::Union{Type{Tri3},Type{Tri6},Type{Tri7},
+        Type{Tet4},Type{Tet10}}, xi)
+    return all(xi .>= 0.0) && (sum(xi) <= 1.0)
+end
+
+function get_reference_coordinates(::E) where E<:AbstractElement{M,B} where {M,B}
+    return get_reference_element_coordinates(B)
+end
+
@@ -0,0 +1,57 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+# Default number of integration points for each element. First rule is the
+# default integration rule returned by `get_integration_points(element)`.
+# Sometimes we want to increase integration order, e.g. when integrating mass
+# matrix or boundary conditions. For that reason, additional rules are provied
+# in list, so e.g. `get_integration_points(element, 1)` returns the second rule,
+# `get_integration_points(element, 2)` third rule and so on. Rules should be
+# ordered so that picking next one integrates more accurately.
+integration_rule_mapping = (
+    :Seg2 => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
+    :Seg3 => (:GLSEG3, :GLSEG4, :GLSEG5),
+    :NSeg => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
+    :Quad4 => (:GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :Quad8 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :Quad9 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :NSurf => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :Hex8 => (:GLHEX8, :GLHEX27, :GLHEX64, :GLHEX125),
+    :Hex20 => (:GLHEX27, :GLHEX64, :GLHEX125),
+    :Hex27 => (:GLHEX27, :GLHEX64, :GLHEX125),
+    :NSolid => (:GLHEX27, :GLHEX64, :GLHEX125),
+    :Tri3 => (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
+    :Tri6 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
+    :Tri7 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
+    :Tet4 => (:GLTET1, :GLTET4, :GLTET5, :GLTET15),
+    :Tet10 => (:GLTET4, :GLTET5, :GLTET15),
+    :Pyr5 => (:GLPYR5,),
+    :Wedge6 => (:GLWED6, :GLWED21),
+    :Wedge15 => (:GLWED21,))
+
+for (E, R) in integration_rule_mapping
+    for i in 1:length(R)
+        P = Val{R[i]}
+        order = Val{i - 1}
+        local code  # Explicitly declare as local to avoid warning
+        if isequal(i, 1)
+            code = quote
+                function get_integration_points(element::$E)
+                    return FEMQuad.get_quadrature_points($P)
+                end
+            end
+        else
+            code = quote
+                function get_integration_points(element::$E, ::Type{$order})
+                    return FEMQuad.get_quadrature_points($P)
+                end
+            end
+        end
+        eval(code)
+    end
+end
+
+# All good codes needs a special case. Here we have it: Poi1
+function get_integration_points(::Poi1)
+    [(1.0, (0.0,))]
+end