JuliaFEM.jl/src/elements/elements.jl

# 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{F, B}

Abstract supertype for all elements.

# Type Parameters
- `F`: Field type (any type-stable container: NamedTuple, struct, etc.)
- `B`: Basis function type (e.g., Lagrange{Triangle,1})
"""
abstract type AbstractElement{F,B} end

"""
    Element{N,NIP,F,B} <: AbstractElement{F,B}

Immutable finite element with type-stable fields.

# Type Parameters
- `N`: Number of nodes (compile-time constant)
- `NIP`: Number of integration points (compile-time constant)
- `F`: Field container type (NamedTuple, struct, etc.) - **must be type-stable!**
- `B`: Basis function type (e.g., Lagrange{Triangle,1})

# Fields
- `id`: Element identifier
- `connectivity`: Node IDs as NTuple (zero-cost, immutable)
- `integration_points`: Integration points as NTuple (zero-cost, immutable)
- `fields`: Type-stable field container (can be empty tuple `()` if no fields)
- `basis`: Basis function type instance

# Design Philosophy
- **Type-stable fields**: `F` parameter ensures compile-time type knowledge
- **Immutable**: All fields immutable (GPU-compatible, thread-safe)
- **Zero-allocation**: NTuple for connectivity and IPs
- **GPU-ready**: Immutable fields can be transferred to GPU without copying

# Examples
```julia
# With fields (NamedTuple):
fields = (E = 210e3, ν = 0.3)
el = Element(UInt(1), (1,2,3), (), fields, Lagrange{Triangle,1}())

# Without fields (empty tuple):
el = Element(UInt(1), (1,2,3), (), (), Lagrange{Triangle,1}())

# Access fields (type-stable!):
E = el.fields.E  # Float64, known at compile time
```
"""
struct Element{N,NIP,F,B} <: AbstractElement{F,B}
    id::UInt
    connectivity::NTuple{N,UInt}  # Tuple for zero-cost, compile-time known size
    integration_points::NTuple{NIP,IP}  # Tuple for zero-cost
    fields::F  # Type-stable field container (NamedTuple, struct, or ())
    basis::B
end

"""
    Element(basis_type, connectivity; fields=(), id=UInt(0))

Construct element from basis type with optional fields.

# Arguments
- `basis_type`: Basis function type (e.g., Lagrange{Triangle,1})
- `connectivity`: Node IDs as tuple or vector
- `fields`: Optional field container (NamedTuple, struct, or empty tuple)
- `id`: Element identifier (default: 0)

# Examples
```julia
# Simple element without fields:
element = Element(Lagrange{Triangle,1}, (1, 2, 3))

# Element with material properties:
element = Element(Lagrange{Triangle,1}, (1, 2, 3), 
                  fields=(E=210e3, ν=0.3))

# Element with complete specification:
element = Element(Lagrange{Quadrilateral,1}, (1,2,3,4),
                  fields=(E=210e3, ν=0.3, thickness=0.01),
                  id=UInt(42))
```
"""
function Element(::Type{B}, connectivity::NTuple{N,<:Integer};
    fields::F=(), id::UInt=UInt(0)) where {N,B<:AbstractBasis,F}
    connectivity_uint = UInt.(connectivity)
    integration_points = ntuple(i -> IP(UInt(0), 0.0, ()), 0)  # Empty initially
    basis = B()
    return Element{N,0,F,B}(id, connectivity_uint, integration_points, fields, basis)
end

function Element(::Type{B}, connectivity::Vector{<:Integer}; kwargs...) where B<:AbstractBasis
    return Element(B, (connectivity...,); kwargs...)
end

# ============================================================================
# Topology → Basis constructors: Accept topology types, create Lagrange basis
# ============================================================================

"""
    Element(::Type{<:AbstractTopology}, connectivity; fields=(), id=UInt(0))

Construct element from topology type. Creates linear Lagrange basis automatically.

# Examples
```julia
element = Element(Triangle, (1, 2, 3))        # → Lagrange{Triangle,1}
element = Element(Quadrilateral, (1,2,3,4),   # → Lagrange{Quadrilateral,1}
                  fields=(E=210e3, ν=0.3))
element = Element(Tet10, (1,2,3,4,5,6,7,8,9,10))  # → Lagrange{Tet10,2}
```
"""
function Element(::Type{T}, connectivity::NTuple{N,<:Integer}; kwargs...) where {N,T<:AbstractTopology}
    # Determine order from number of nodes
    order = infer_lagrange_order(T, N)
    BasisType = Lagrange{T,order}
    return Element(BasisType, connectivity; kwargs...)
end

function Element(::Type{T}, connectivity::Vector{<:Integer}; kwargs...) where T<:AbstractTopology
    return Element(T, (connectivity...,); kwargs...)
end

"""Infer Lagrange order from topology type and number of nodes"""
function infer_lagrange_order(::Type{T}, n::Int) where T<:AbstractTopology
    # Get the type name for pattern matching
    tname = string(nameof(T))

    # Check explicit higher-order topologies first (Tet10, Hex20, Tri6, etc.)
    if tname == "Tet10" && n == 10
        return 2
    elseif tname == "Tri6" && n == 6
        return 2
    elseif tname == "Tri7" && n == 7
        return 3
    elseif tname == "Quad8" && n == 8
        return 2  # Serendipity
    elseif tname == "Quad9" && n == 9
        return 2  # Full quadratic
    elseif tname == "Hex20" && n == 20
        return 2  # Serendipity
    elseif tname == "Hex27" && n == 27
        return 2  # Full quadratic
    elseif tname == "Wedge15" && n == 15
        return 2
    elseif tname == "Seg3" && n == 3
        return 2
    end

    # For base topologies (and their aliases), infer order from number of nodes
    if T === Segment || T === Seg2
        n == 2 && return 1
        n == 3 && return 2
    elseif T === Triangle || T === Tri3
        n == 3 && return 1
        n == 6 && return 2
        n == 7 && return 3
    elseif T === Quadrilateral || T === Quad4
        n == 4 && return 1
        n == 8 && return 2  # Serendipity
        n == 9 && return 2  # Full
    elseif T === Tetrahedron || T === Tet4
        n == 4 && return 1
        n == 10 && return 2
    elseif T === Hexahedron || T === Hex8
        n == 8 && return 1
        n == 20 && return 2  # Serendipity
        n == 27 && return 2  # Full
    elseif T === Pyramid || T === Pyr5
        n == 5 && return 1
    elseif T === Wedge || T === Wedge6
        n == 6 && return 1
        n == 15 && return 2
    end

    # Default to order 1 (linear)
    return 1
end
# Commenting out until new architecture is fully implemented.

# # Helper: Map topology type to basis type
# topology_to_basis(::Type{Poi1}) = Poi1Basis
# topology_to_basis(::Type{Seg2}) = Seg2Basis
# topology_to_basis(::Type{Seg3}) = Seg3Basis
# topology_to_basis(::Type{Tri3}) = Tri3Basis
# topology_to_basis(::Type{Tri6}) = Tri6Basis
# topology_to_basis(::Type{Tri7}) = Tri7Basis
# topology_to_basis(::Type{Quad4}) = Quad4Basis
# topology_to_basis(::Type{Quad8}) = Quad8Basis
# topology_to_basis(::Type{Quad9}) = Quad9Basis
# topology_to_basis(::Type{Tet4}) = Tet4Basis
# topology_to_basis(::Type{Tet10}) = Tet10Basis
# topology_to_basis(::Type{Hex8}) = Hex8Basis
# topology_to_basis(::Type{Hex20}) = Hex20Basis
# topology_to_basis(::Type{Hex27}) = Hex27Basis
# topology_to_basis(::Type{Pyr5}) = Pyr5Basis
# topology_to_basis(::Type{Wedge6}) = Wedge6Basis
# topology_to_basis(::Type{Wedge15}) = Wedge15Basis

# # Shim constructors that forward to basis-based constructors
# function Element(::Type{T}, connectivity::NTuple{N,<:Integer}) where {N,T<:AbstractTopology}
#     return Element(topology_to_basis(T), connectivity)
# end

# function Element(::Type{T}, ::Type{M}, connectivity::NTuple{N,<:Integer}) where {N,M<:AbstractFieldSet,T<:AbstractTopology}
#     return Element(topology_to_basis(T), M, connectivity)
# end

# function Element(::Type{T}, connectivity::Vector{<:Integer}) where T<:AbstractTopology
#     return Element(topology_to_basis(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

# ============================================================================
# COMPATIBILITY SHIM: Old dfields API → New fields API
# ============================================================================
# NOTE: This is temporary until field system is redesigned (Phase 3)
# Old API: element.dfields (mutable Dict for dynamic fields)
# New API: element.fields (type-stable NamedTuple/struct, immutable)
#
# For now, we map dfields to the same place as sfields (element.fields)
# since the new API doesn't distinguish between static and dynamic fields.

function has_dfield(element, field_name)
    # In new API, dfields don't exist as separate from sfields
    # Check if field_name is in element.fields
    F = typeof(element.fields)
    if F === Tuple{}
        return false  # Empty fields
    elseif F <: NamedTuple
        return haskey(element.fields, field_name)
    else
        return isdefined(element.fields, field_name)
    end
end

function get_dfield(element, field_name)
    # Return the field value from element.fields
    F = typeof(element.fields)
    if F === Tuple{}
        error("Element has no fields (empty tuple)")
    end
    return getfield(element.fields, field_name)
end

function create_dfield!(element, field_name, field_::AbstractField)
    # WARNING: In new API, fields are IMMUTABLE!
    # This function cannot actually create the field in-place
    # TODO Phase 3: Redesign field system to support mutable time-varying fields
    @warn "create_dfield!() called but fields are immutable in new API. Field not created." maxlog = 1
    return
end

function create_dfield!(element, field_name, field_data)
    # WARNING: In new API, fields are IMMUTABLE!
    @warn "create_dfield!() called but fields are immutable in new API. Field not created." maxlog = 1
    return
end

function update_dfield!(element, field_name, field_data)
    # WARNING: In new API, fields are IMMUTABLE!
    # This function cannot actually update the field in-place
    # TODO Phase 3: Redesign field system to support mutable time-varying fields
    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
    return
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
    # WARNING: In new API, fields are IMMUTABLE!
    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
    return
end

function update_dfield!(element, field_name, field_data::Dict{Int,V}) where V
    # WARNING: In new API, fields are IMMUTABLE!
    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
    return
end

function update_dfield!(element, field_name, field_data::Function)
    # WARNING: In new API, fields are IMMUTABLE!
    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
    return
end

function interpolate_dfield(element, field_name, time)
    # In new API, fields are static (no time variation)
    # Just return the field value
    return get_dfield(element, field_name)
end

### sfields statically defined fields

# ============================================================================
# COMPATIBILITY SHIM: Old sfields/dfields API → New fields API
# ============================================================================
# NOTE: This is temporary until field system is redesigned (Phase 3)
# Old API: element.sfields (static fields) and element.dfields (dynamic fields)
# New API: element.fields (type-stable NamedTuple/struct, immutable)
#
# Problem: Tests use update_field!() to mutate fields, but new fields are immutable
# Solution: Make has_sfield/get_sfield work with new API, but warn about mutations

function has_sfield(element, field_name)
    # In new API, sfields don't exist - but we can check if field_name is in element.fields
    F = typeof(element.fields)
    if F === Tuple{}
        return false  # Empty fields
    elseif F <: NamedTuple
        return haskey(element.fields, field_name)
    else
        return isdefined(element.fields, field_name)
    end
end

function get_sfield(element, field_name)
    # Return the field value from element.fields
    F = typeof(element.fields)
    if F === Tuple{}
        error("Element has no fields (empty tuple)")
    end
    return getfield(element.fields, field_name)
end

function update_sfield!(element, field_name, field_data)
    # WARNING: In new API, fields are IMMUTABLE!
    # This function cannot actually update the field in-place
    # TODO Phase 3: Redesign field system to support mutable time-varying fields
    @warn "update_sfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
    return nothing
end

function interpolate_sfield(element, field_name, time)
    # In new API, fields are static (no time variation)
    # Just return the field value
    return get_sfield(element, field_name)
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)
    # In new API: fields are immutable, so we can't actually update them
    # For now, just check if field exists and warn
    if has_sfield(element, field_name)
        update_sfield!(element, field_name, field_data)
    elseif has_dfield(element, field_name)
        update_dfield!(element, field_name, field_data)
    else
        # Field doesn't exist - this is OK in new API where fields are optional
        # Silently ignore (many tests add fields dynamically that we don't need)
    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}
    # Handle both raw coordinates (Tuple) and IP struct
    coords = isa(ip, IP) ? ip.coords : ip
    T = typeof(first(coords))
    # Convert to Vec for Tensors.jl compatibility
    xi = Vec{length(coords),T}(coords)
    # eval_basis! now returns a tuple directly - zero allocations!
    N_tuple = eval_basis!(B, T, xi)
    # Return as row matrix for compatibility with old code
    # This still allocates, but only at the API boundary
    return reshape(collect(N_tuple), 1, length(element))
end

function get_dbasis(element::AbstractElement{M,B}, ip, ::Any) where {M,B}
    # Handle both raw coordinates (Tuple) and IP struct
    coords = isa(ip, IP) ? ip.coords : ip
    T = typeof(first(coords))
    # Convert to Vec for Tensors.jl compatibility
    xi = Vec{length(coords),T}(coords)
    # eval_dbasis! now returns NTuple{N,Vec{D}} directly - zero allocations!
    dN_tuple = eval_dbasis!(B, xi)
    # Return as Vector for compatibility with old code
    # This still allocates, but only at the API boundary
    return collect(dN_tuple)
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_dict = element("geometry", time)
    # Convert to Vector{Vec} for Tensors.jl compatibility
    X = [Vec(x...) for x in X_dict]
    # 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

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
    # 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 m == 1  # 1D element (boundary of 2D or 3D), J is n×1, JT is 1×n
        return norm(JT)
    else  # 2D element (manifold on 3D problem), J is 3×2, JT is 2×3
        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::Element{N,NIP,M,B}) where {N,NIP,M,B}
    # If integration points already set, return them
    if NIP > 0
        return element.integration_points
    end
    # Otherwise get default integration points for this element type
    ips = get_integration_points(element.properties)
    return tuple([IP(UInt(i), w, xi) for (i, (w, xi)) in enumerate(ips)]...)
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 tuple([IP(UInt(i), w, xi) for (i, (w, xi)) in enumerate(ips)]...)
end

""" 
    with_integration_points(element, integration_points_tuple) -> Element

Create a new element with the given integration points. Since Element is immutable,
this returns a new instance with updated integration points.
"""
function with_integration_points(element::Element{N,NIP,M,B}, ips::NTuple{NNEW,IP}) where {N,NIP,M,B,NNEW}
    return Element{N,NNEW,M,B}(element.id, element.connectivity, ips,
        element.dfields, element.sfields, element.properties)
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

# OLD: Uses BasisInfo which was defined in math.jl (commented out)
# 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