Elements

In JuliaFEM, elements can be considered as "containers", combining fields and basis functions described above. Among that, element has information about topology (connectivity) and numerical integration rule. These fundamentals forms a finite element, the backbone of finite element method.

New elements are constructed by passing basis type (e.g. Seg2, Quad4, Tet10, ...) to Element and list of id numbers to where element is topologically connected.

julia> el = Element(Quad4, [1, 2, 3, 4])
Element{FEMBase.EmptyFieldSet{4},Quad4}(-1, [1, 2, 3, 4], FEMBase.Point{IntegrationPoint}[], Dict{Symbol,AbstractField}(), FEMBase.EmptyFieldSet{4}(), Quad4())

Setting fields to element is done using a command update!, which either creates a new field if field with that name does not already exist, or updates the old one. Typically, at least field called geometry needs to be defined to element as it's used to calculate Jacobian of element. Fields can be discrete, continuous, time invariant or variant, variable or constant, or anything else that is subclassed from AbstractField.

julia> X = Dict(1 => [0.0,0.0], 2=>[1.0,0.0], 3=>[1.0,1.0], 4=>[0.0,1.0])
Dict{Int64,Array{Float64,1}} with 4 entries:
  4 => [0.0, 1.0]
  2 => [1.0, 0.0]
  3 => [1.0, 1.0]
  1 => [0.0, 0.0]

julia> update!(el, "geometry", X)

Internally, fields are stored in a Dict:

julia> el.fields
Dict{Symbol,AbstractField} with 1 entry:
  :geometry => DVTI{4,Array{Float64,1}}(([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0…

The command update! is automatically inspecting field type based in input. For example, to define temporal field varying spatially:

julia> u0 = ([0.0,0.0], [0.0,0.0], [0.0,0.0], [0.0,0.0])
([0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0])

julia> u1 = ([0.0,0.0], [0.0,0.0], [0.5,0.0], [0.0,0.0])
([0.0, 0.0], [0.0, 0.0], [0.5, 0.0], [0.0, 0.0])

julia> update!(el, "displacement", 0.0 => u0)

julia> update!(el, "displacement", 1.0 => u1)
2-element Array{Pair{Float64,NTuple{4,Array{Float64,1}}},1}:
 0.0 => ([0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0])
 1.0 => ([0.0, 0.0], [0.0, 0.0], [0.5, 0.0], [0.0, 0.0])

julia> el.fields
Dict{Symbol,AbstractField} with 2 entries:
  :geometry     => DVTI{4,Array{Float64,1}}(([0.0, 0.0], [1.0, 0.0], [1.0, 1.0]…
  :displacement => DVTV{4,Array{Float64,1}}(Pair{Float64,NTuple{4,Array{Float64…

Interpolating of fields can be done using syntax Element(field_name, xi, time). For example, position of material particle $X$ in initial configuration and deformed configuration in the middle of the element at time $t=1$ can be found as

julia> xi = (0.0, 0.0)
(0.0, 0.0)

julia> time = 1.0
1.0

julia> X = el("geometry", xi, time)
2-element Array{Float64,1}:
 0.5
 0.5

julia> u = el("displacement", xi, time)
2-element Array{Float64,1}:
 0.125
 0.0

julia> x = X+u
2-element Array{Float64,1}:
 0.625
 0.5

Jacobian, determinant of Jacobian and gradient of field can be calculated adding extra argument Val{:Jacobian}, Val{:detJ}, Val{:Grad} to the above command and not passing field name, i.e.

julia> el(xi, time, Val{:Jacobian})
2×2 Array{Float64,2}:
 0.5  0.0
 0.0  0.5

julia> el(xi, time, Val{:detJ})
0.25

julia> el(xi, time, Val{:Grad})
2×4 Array{Float64,2}:
 -0.5   0.5  0.5  -0.5
 -0.5  -0.5  0.5   0.5

Usually what is needed when calculating local stiffness matrices is a gradient of some field. For example, displacement gradient and temperature gradient can be obtained following way:

julia> gradu = el("displacement", xi, time, Val{:Grad})
2×2 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
 0.25  0.25
 0.0   0.0

julia> update!(el, "temperature", (1.0, 2.0, 3.0, 4.0))

julia> gradT = el("temperature", xi, time, Val{:Grad})
1×2 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 0.0  2.0

Accessing integration points of element is done using function get_integration_points. Combining interpolation and integration one can already calculate local matrices of a single element or, for example area and strain energy:

update!(el, "lambda", 96.0)
update!(el, "mu", 48.0)

A = 0.0
W = 0.0
for ip in get_integration_points(el)
    detJ = el(ip, time, Val{:detJ})
    A += ip.weight * detJ
    ∇u = el("displacement", ip, time, Val{:Grad})
    E = 1/2*(∇u + ∇u')
    λ = el("lambda", ip, time)
    μ = el("mu", ip, time)
    W += ip.weight * ( λ/2*trace(E*E') + μ*trace(E)^2) * detJ
end

println("Area: $A")
println("Strain energy: $W")

To calculate local stiffness matrix for Poisson problem:

K = zeros(4,4)
update!(el, "coefficient", 36.0)
for ip in get_integration_points(el)
    dN = el(ip, time, Val{:Grad})
    detJ = el(ip, time, Val{:detJ})
    c = el("coefficient", ip, time)
    K += ip.weight * c*dN'*dN * detJ
end
K

For more memory efficient code it's a good idea to use BasisInfo and element_info! which allocates working memory to calculate all "basic stuff" for a single integration point, like Jacobian, determinant of Jacobian, basis and it's partial derivatives with respect to reference configuration $X$.

bi = BasisInfo(Quad4)
fill!(K, 0.0)
for ip in get_integration_points(el)
    J, detJ, N, dN = element_info!(bi, el, ip, time)
    c = el("coefficient", ip, time)
    K += ip.weight * c*dN'*dN * detJ
end
K

Using analytical fields

Creating fields depending from other fields