introduced new way of handling fields, now including time.

src/JuliaFEM.jl

@@ -2,8 +2,8 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 module JuliaFEM
 
-VERSION < v"0.4-" && using Docile
 using Lexicon
+
 using Logging
 @Logging.configure(level=DEBUG)
 
@@ -12,8 +12,9 @@ include("elements.jl") # elements
 include("equations.jl") # formulations
 include("problems.jl") # problems
 
-include("math.jl") # basic mathematical operations -- obsolete ..?
 
+
+include("math.jl") # basic mathematical operations -- obsolete ..?
 include("elasticity_solver.jl")
 include("xdmf.jl")
 include("abaqus_reader.jl")

src/types.jl

@@ -1,17 +1,104 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-type Assembly
-    # LHS
-    I :: Array{Int64, 1}
-    J :: Array{Int64, 1}
-    A :: Array{Float64, 1}
-    # RHS
-    i :: Array{Int64, 1}
-    b :: Array{Float64, 1}
-    # global dofs for each element
-    gdofs :: Dict{Int64, Array{Int64, 1}}
-end
-Assembly() = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())
-Assembly(gdofs::Dict{Int64,Array{Int64,1}}) = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)
+# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
+
+using ForwardDiff
+
+
+""" Field. """
+type Field{T}
+    time :: Float64
+    increment :: Int64
+    values :: Array{T, 1}
+end
+
+""" Initialize field. """
+function Field(time, values)
+    Field(time, 1, values)
+end
+
+""" Get length of a field (number of basis functions in practice). """
+Base.length(f::Field) = length(f.values)
+
+""" Get field discrete value at point i. """
+Base.getindex(f::Field, i::Int64) = f.values[i]
+
+""" Interpolate field h(ξ)*f = x*f """
+Base.(:*)(x::Array{Float64, 1}, f::Field) = sum(x .* f.values)
+Base.(:*)(x::Array{Float64, 2}, f::Field) = sum([f[i]*x[i,:] for i in 1:length(f)])
+
+""" Interpolate field (h*f)(ξ) """
+Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
+
+""" Multiply field with some constant k. """
+Base.(:*)(k::Float64, f::Field) = Field(f.time, k*f.values)
+
+""" Sum two fields. """
+function Base.(:+)(f1::Field, f2::Field)
+    @assert(f1.time == f2.time, "Cannot add fields: time mismatch, $(f1.time) != $(f2.time)")
+    Field(f1.time, f1.values + f2.values)
+end
+
+""" Interpolate from set of fields x*[f1, f2] where x is evaluated basis. """
+Base.(:*)(x::Array{Float64, 1}, f::Array{Field}) = sum(x .* f)
+
+""" Interpolate from set of fields with basis b, i.e. f(t) = b(t)*[f1, f2] """
+Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
+
+"""
+Interpolate a field from finite set of fields some time t ∈ R.
+"""
+function call(fields :: Array{Field, 1}, t::Float64)
+    if t <= fields[1].time
+        return fields[1]
+    end
+    if t >= fields[end].time
+        return fields[end]
+    end
+    i = length(fields)
+    while fields[i].time >= t
+        i -= 1
+    end
+    if fields[i].time == t
+        return fields[i]
+    end
+    #Logging.debug("doing linear interpolation between fields $i and $(i+1)")
+    f1 = fields[i]
+    t1 = f1.time
+    f2 = fields[i+1]
+    t2 = f2.time
+    dt = t2 - t1
+    nw = (t2-t)/dt*f1.values + (t-t1)/dt*f2.values
+    f = Field(t, nw)
+    return f
+end
+function call(field::Field, t::Float64)
+    Field(t, field.increment, field.values)
+end
+
+
+
+""" Basis function. """
+type Basis
+    basis :: Function
+    dbasisdxi :: Function
+end
+
+""" Constructor of basis function. """
+function Basis(basis)
+    Basis(basis, ForwardDiff.jacobian(basis))
+end
+
+""" Interpolate field f using basis b. """
+Base.(:*)(b::Basis, f::Field) = (x) -> b(x)*f
+Base.(:*)(b::Basis, f::Array{Field}) = (t) -> b(t)*f
+
+""" Evaluate basis function in point ξ. """
+call(b::Basis, xi) = b.basis(xi)
+
+""" Get partial derivative of basis function. """
+∂(h::Basis) = h.dbasisdxi
+diff(h::Basis) = h.dbasisdxi
+derivative(h::Basis) = h.dbasisdxi
 

test/test_types.jl

@@ -0,0 +1,60 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM: Basis, Field, get_field, diff
+using FactCheck
+
+facts("test fields and interpolation") do
+
+    # simple interpolation in domain [-1, 1]
+    N = Basis((ξ) -> [0.5*(1.0-ξ[1]), 0.5*(1.0+ξ[1])])
+    u = Field(0.0, [0.0, 1.0])
+    @fact N([0.0])*u --> 0.5
+    @fact (N*u)([0.0]) --> 0.5
+
+    # multiply of field with constant
+    u1 = Field(0.0, [0.0, 1.0])
+    u2 = 3.0*u1
+    @fact u1.time --> 0.0
+    @fact u2.time --> 0.0
+    @fact u2.values --> [0.0, 3.0]
+
+    # addition of fields together
+    u1 = Field(0.0, [0.0, 1.0])
+    u2 = Field(0.0, [1.0, 2.0])
+    u3 = u1 + u2
+    @fact u3.values --> [1.0, 3.0]
+
+    # interpolation between two fields in time domain
+    u1 = Field(0.0, [0.0, 1.0])
+    u2 = Field(0.0, [1.0, 2.0])
+    t = Basis((t) -> [1-t, t])
+    u = Field[u1, u2]
+    u2 = (t*u)(0.5)
+    @fact u2.values --> [0.5, 1.5]
+
+    # interpolation between set of fields
+    u1 = Field(0.0, [0.0, 0.0])
+    u2 = Field(1.0, [1.0, 2.0])
+    u3 = Field(2.0, [0.5, 1.5])
+    u = Field[u1, u2, u3]
+    @fact u(-1.0).values --> [0.0, 0.0]
+    @fact u(0.0).values --> [0.0, 0.0]
+    @fact u(1.0).values --> [1.0, 2.0]
+    @fact u(2.0).values --> [0.5, 1.5]
+    @fact u(3.0).values --> [0.5, 1.5]
+    @fact u(0.5).values --> [0.5, 1.0]
+    @fact u(1.5).values --> [0.75, 1.75]
+
+    X = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    h = Basis((xi) ->
+    [(1-xi[1])*(1-xi[2])/4
+     (1+xi[1])*(1-xi[2])/4
+     (1+xi[1])*(1+xi[2])/4
+     (1-xi[1])*(1+xi[2])/4])
+   @fact (h*X)([0.0, 0.0]) --> [0.5, 0.5]
+   @fact h([0.0, 0.0])*X --> [0.5, 0.5]
+   @fact diff(h)([0.0, 0.0])*X --> [0.5 0.0; 0.0 0.5]
+   @fact (diff(h)*X)([0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
+
+end