issue #67

src/JuliaFEM.jl

@@ -52,8 +52,8 @@ function Base.linspace{T<:Array}(X1::T, X2::T, n)
 end
 
 # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
-include("fields.jl")
-include("basis.jl")  # interpolation of discrete fields
+include("fields2.jl")
+#include("basis.jl")  # interpolation of discrete fields
 include("symbolic.jl") # a thin symbolic layer for fields
 include("types.jl")  # type definitions
 

src/elasticity.jl

@@ -114,7 +114,7 @@ end
 """ 4-node plane stress element. """
 type CPS4 <: PlaneStressElasticityEquation
     element :: Quad4
-    integration_points :: Array{IntegrationPoint, 1}
+    integration_points :: Vector{IntegrationPoint}
 end
 
 function Base.size(equation::CPS4)
@@ -123,7 +123,9 @@ end
 
 function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4)
     integration_points = get_default_integration_points(element)
-    haskey(element, "displacement") || (element["displacement"] = zeros(2, 4))
+    if !haskey(element, "displacement")
+        element["displacement"] = 0.0 => [zeros(2) for i=1:4]
+    end
     CPS4(element, integration_points)
 end
 
@@ -139,7 +141,9 @@ end
 
 function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2)
     integration_points = get_default_integration_points(element)
-    haskey(element, "displacement") || (element["displacement"] = zeros(2, 2))
+    if !haskey(element, "displacement")
+        element["displacement"] = 0.0 => [zeros(2) for i=1:2]
+    end
     CPS2(element, integration_points)
 end
 

src/elements.jl

@@ -75,18 +75,21 @@ Examples
 >>> element["temperature"] = (0.0, [0, 0, 0, 0]), (1.0, [1, 2, 3, 4])
 >>> element["temperature"] = (0.0 => [0, 0, 0, 0], 1.0 => [1, 2, 3, 4])
 """
-function Base.setindex!(element::Element, field_data, field_name)
-    setindex!(element.fields, field_data, field_name)
+function Base.setindex!(element::Element, data, name::ASCIIString)
+    element.fields[name] = Field(data)
+end
+function Base.setindex!(element::Element, data::Tuple, name::ASCIIString)
+    element.fields[name] = Field(data...)
 end
 
-function Base.setindex!(element::Element, field_data::Tuple, field_name)
-    field = Field()
-    for (time, data) in field_data
-        ts = TimeStep(time, Increment[Increment(data)])
-        push!(field, ts)
-    end
-    element[field_name] = field
-end
+#function Base.setindex!(element::Element, field_data::Tuple, field_name)
+#    field = Field()
+#    for (time, data) in field_data
+#        ts = TimeStep(time, Increment[Increment(data)])
+#        push!(field, ts)
+#    end
+#    element[field_name] = field
+#end
 
 function get_connectivity(el::Element)
     return el.connectivity
@@ -95,12 +98,12 @@ end
 abstract AbstractFunctionSpace
 
 type FunctionSpace <: AbstractFunctionSpace
-    basis :: Basis
+    basis :: CVTI
     fields :: FieldSet
 end
 
 type GradientFunctionSpace <: AbstractFunctionSpace
-    basis :: Basis
+    basis :: CVTI
     fields :: FieldSet
 end
 
@@ -130,9 +133,6 @@ end
 """ Evaluate field on element function space. """
 function call(u::FunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
     field = !isa(variation, Void) ? variation : u.fields[field_name](t)
-    if length(field) == 1
-        return field.data[1]
-    end
     u.basis(field, xi)
 end
 
@@ -171,7 +171,7 @@ end
 """ Get a determinant of element in point ξ. """
 function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
     X = u.fields["geometry"](time)
-    dN = u.basis.dbasisdxi(xi)
+    dN = u.basis(xi, Val{:grad})
     J = sum([dN[:,i]*X[i]' for i=1:length(X)])
     m, n = size(J)
     return m == n ? det(J) : norm(J)

src/equations.jl

@@ -161,7 +161,8 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
 
     # 3. virtual work -- user has defined some residual r = p - f = 0
     if has_residual_vector(equation)
-        field = element[unknown_field_name](time)
+
+        field = DVTI(last(element[unknown_field_name]).data)
 
         """ Wrapper for virtual work for ForwardDiff. """
         function calc_R(data::Vector)
@@ -177,9 +178,12 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
             if haskey(element, "$unknown_field_name nodal load")
                 R -= vec(element["$unknown_field_name nodal load"](time))
             end
+            #info("return = $R")
             return R
         end
-       
+
+        #info("field = $field")
+        #info("vec(field) = $(vec(field))")
         jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
         add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
         add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults))

src/fields2.jl

@@ -16,7 +16,48 @@ type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{
     data
 end
 
-# Different field combinations
+### Basic data structure for discrete field
+
+type Increment{T}
+    time :: Float64
+    data :: T
+end
+
+function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T})
+    return Increment{T}(data[1], data[2])
+end
+
+function Base.convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}})
+    time = data[1]
+    content = data[2]
+    return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)])
+end
+
+function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
+    return increment.data[i]
+end
+
+function Base.(:*)(d, increment::Increment)
+    return d*increment.data
+end
+
+### Basic data structure for continuous field
+
+type Basis
+    basis :: Function
+    dbasis :: Function
+end
+
+function Base.call(basis::Basis, xi::Vector)
+    basis.basis(xi)
+end
+
+function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
+    basis.dbasis(xi)
+end
+
+### Different field combinations and other typealiases
+
 typealias DCTI Field{Discrete,   Constant, TimeInvariant}
 typealias DVTI Field{Discrete,   Variable, TimeInvariant}
 typealias DCTV Field{Discrete,   Constant, TimeVariant}
@@ -26,71 +67,207 @@ typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to inter
 typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time
 typealias CVTV Field{Continuous, Variable, TimeVariant}
 
-# Basic data structure for discrete field
-type Increment{T}
-    time :: Float64
-    data :: T
-end
+typealias ScalarIncrement{T} Increment{T}
+typealias VectorIncrement{T} Increment{Vector{T}}
+typealias TensorIncrement{T} Increment{Matrix{T}}
 
-typealias VectorIncrement Increment{Vector}
+typealias DiscreteField      Union{DCTI, DVTI, DCTV, DVTV}
+typealias ContinuousField    Union{CCTI, CVTI, CCTV, CVTV}
+typealias ConstantField      Union{DCTI, DCTV, CCTI, CCTV}
+typealias VariableField      Union{DVTI, DVTV, CVTI, CVTV}
+typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI}
+typealias TimeVariantField   Union{DCTV, DVTV, CCTV, CVTV}
 
-function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
-    return increment.data[i]
-end
 
-# Basic data structure for continuous field
-type Basis
-    basis :: Function
-    dbasis :: Function
-end
+### Convenient functions to create fields
 
-# Functions simplifying definition of fields.
-
-"""
-All other data than vectors are considered as constant time invariant fields.
-"""
 function Field(data)
-    DCTI(data)
+    return DCTI(data)
 end
 
-"""
-Vector data is considered as variable field time invariant field.
-"""
 function Field(data::Vector)
-    DVTI(data)
+    return DVTI(data)
 end
 
-"""
-Data given in (time, value) pairs, where value is not vector, is considered as
-constant time variant field.
-"""
 function Field{T}(data::Pair{Float64, T}...)
-    increments = [Increment{T}(d[1], d[2]) for d in data]
-    DCTV(increments)
+    return DCTV([Increment{T}(d[1], d[2]) for d in data])
 end
 
-"""
-Data given in (time, value) pairs, where value is a vector, is considered as
-variable time variant field.
-"""
 function Field{T}(data::Pair{Float64, Vector{T}}...)
-    increments = [Increment{Vector{T}}(d[1], d[2]) for d in data]
-    DVTV(increments)
+    return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
-""" Special case, constant time-variant vector, converted automatically. """
 function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...)
-    increments = [Increment(d[1], d[2]) for d in data]
-    DCTV(increments)
+    return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
-## Other field related functions
+function CVTI(basis::Function, dbasis::Function)
+    return CVTI(Basis(basis, dbasis))
+end
+
+function Field(basis::Function, dbasis::Function)
+    return CVTI(basis, dbasis)
+end
+
+### Accessing and manipulating discrete fields
 
 function Base.getindex(field::DVTV, i::Int64)
     return field.data[i]
 end
 
+function Base.push!(field::DCTV, data::Pair)
+    push!(field.data, data)
+end
+
+function Base.push!(field::DVTV, data::Pair)
+#    info("field.data = \n$(field.data)")
+#    info("data = \n$data")
+    push!(field.data, data)
+end
+
+function Base.getindex(field::DVTV, i::Int64)
+    return field.data[i]
+end
+
+function Base.getindex(field::DVTI, i::Int64)
+    return field.data[i]
+end
+
+function Base.getindex(field::DCTV, i::Int64)
+    return field.data[i]
+end
+
+function Base.getindex(field::Field, i::Int64)
+    return field.data[i]
+end
+
+function Base.length(field::DVTI)
+    return length(field.data)
+end
+
+function Base.length(field::DCTI)
+    return 1
+end
+
+function Base.length(field::DVTV)
+    return length(field.data)
+end
+
+function Base.length(field::DCTV)
+    return length(field.data)
+end
+
+for op = (:+, :*, :/, :-)
+    @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data)
+    @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data)
+    @eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data)
+    @eval ($op)(field::DCTI, k) = ($op)(field.data, k)
+    @eval ($op)(k, field::DCTI) = ($op)(field.data, k)
+end
+
+function Base.vec(field::DVTI)
+    return [field.data...;]
+end
+
+function Base.vec(field::DCTV)
+    info("trying to vectorize $field")
+    error("does not make sense")
+end
+
+function Base.endof(field::Field)
+    return endof(field.data)
+end
+
+#function Base.similar{T}(field::DVTI, data::Vector{T})
+#    return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
+#end
+
+function Base.similar{T}(field::DVTI, data::Vector{T})
+    n = length(field.data)
+    data = reshape(data, round(Int, length(data)/n), n)
+    newdata = Vector[data[:,i] for i=1:n]
+    return typeof(field)(newdata)
+end
+
+### Accessing continuous fields
+
+function Base.call(field::CVTI, xi::Vector)
+    field.data(xi)
+end
+
+function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}})
+    field.data(xi, Val{:grad})
+end
+
+function Base.convert(::Type{Basis}, field::CVTI)
+    return field.data
+end
+
+### Interpolation 
+
+function Base.call(field::DVTI, time::Float64)
+    # interpolating time-invariant field in time direction -> no effect
+    return field
+end
+
+function Base.call(field::DCTI, time::Float64)
+    # interpolating time-invariant field in time direction -> no effect
+    return field
+end
+
+function Base.call(basis::CVTI, field::DCTI, xi::Vector)
+    # try to interpolate constant value -> no effect
+    return field
+end
+
+#function Base.call(basis::Basis, field::DCTI, xi::Vector)
+    # calling constant field with basis -> no effect
+#    return field
+#end
+
+function Base.call(basis::CVTI, values::DVTI, xi::Vector)
+    N = basis(xi)
+    return sum([N[i]*values[i] for i=1:length(N)])
+end
+
+function Base.call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
+    dbasis = basis(xi, Val{:grad})
+    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
+    invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
+    grad = invJ * dbasis
+    return grad
+end
+
+function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
+    grad = call(basis, geometry, xi, Val{:grad})
+    gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
+    return length(gradf) == 1 ? gradf[1] : gradf
+end
+
+function Base.call(field::DCTV, time::Float64)
+    for i in length(field)
+        if isapprox(field[i].time, time)
+            return DCTI(field[i].data)
+        end
+    end
+    error("interpolate DCTV: not implemented yet")
+end
+
+function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
+#    info("length of field DVTV: $(length(field))")
+    for i=reverse(1:length(field))
+        res = isapprox(field[i].time, time)
+        #info("isapprox $(field[i].time) to $time ? $res")
+        if isapprox(field[i].time, time)
+            return DVTI(field[i].data)
+        end
+    end
+    info(field.data)
+    info(time)
+    error("interpolate DVTV: not implemented yet")
+end
 
 ### FIELDSET ###
 
 typealias FieldSet Dict{ASCIIString, Field}
+

src/heat.jl

@@ -60,6 +60,7 @@ function assemble!(assembly::Assembly, equation::HeatEquation, time::Number=0.0,
             add!(assembly.force_vector, gdofs, w*N'*f)
         end
         if haskey(element, "temperature flux")
+            info("assemble boundary flux")
             g = basis("temperature flux", ip, time)
             add!(assembly.force_vector, gdofs, w*N'*g)
         end
@@ -93,13 +94,13 @@ end
 
 function Base.convert(::Type{HeatEquation}, element::Quad4)
     integration_points = get_default_integration_points(element)
-    haskey(element, "temperature") || (element["temperature"] = zeros(4))
+    haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4))
     DC2D4(element, integration_points)
 end
 
 function Base.convert(::Type{HeatEquation}, element::Seg2)
     integration_points = get_default_integration_points(element)
-    haskey(element, "temperature") || (element["temperature"] = zeros(2))
+    haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2))
     DC2D2(element, integration_points)
 end
 

src/lagrange.jl

@@ -36,11 +36,11 @@ macro create_lagrange_element(element_name, element_description, X, P)
         basis, dbasisdxi = calculate_lagrange_basis($P, $X)
         type $eltype <: CG
             connectivity :: Array{Int, 1}
-            basis :: Basis
+            basis :: CVTI
             fields :: FieldSet
         end
         function $eltype(connectivity, args...)
-            $eltype(connectivity, Basis(basis, dbasisdxi), FieldSet())
+            $eltype(connectivity, CVTI(basis, dbasisdxi), FieldSet())
         end
         get_element_description(el::Type{$eltype}) = $element_description
         Base.size(element::Type{$eltype}) = Base.size($X)

src/solvers.jl

@@ -9,8 +9,7 @@ abstract Solver
 Solve field equations for single element with some dofs fixed. This can be used
 to test nonlinear element formulations.
 """
-function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number=0.0; 
-                 max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
     unknown_field_name = get_unknown_field_name(equation)
     element = get_element(equation)
     x0 = element[unknown_field_name](0.0)
@@ -28,7 +27,9 @@ function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number=0.0;
         end
         dx[free_dofs] = A \ b
         x += dx
-        push!(element[unknown_field_name], reshape(x, size(equation)))
+        eqsize = size(equation)
+        data = eqsize[1] != 1 ? reshape(x, eqsize) : x
+        push!(element[unknown_field_name], time => data)
         norm(dx) < tolerance && return
     end
     error("Did not converge in $max_iterations iterations")
@@ -40,7 +41,7 @@ to test nonlinear element formulations. Dirichlet boundary is assumed to be homo
 and degrees of freedom are eliminated. So if boundary condition is known in nodal
 points and everything is zero this should be quite good.
 """
-function solve!(problem::Problem, free_dofs::Vector{Int}, time::Number=1.0; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
     info("start solver")
     assembly = Assembly()
 #   x = zeros(ga.ndofs)
@@ -68,8 +69,12 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Number=1.0; max_
         for equation in get_equations(problem)
             element = get_element(equation)
             gdofs = get_gdofs(equation)
-            data = reshape(full(x[gdofs]), size(equation))
-            push!(element[field_name], data)
+            data = full(x[gdofs])
+            eqsize = size(equation)
+            if eqsize[1] != 1
+                data = reshape(data, eqsize)
+            end
+            push!(element[field_name], time => data)
         end
         norm(dx) < tolerance && return
     end

src/symbolic.jl

@@ -3,7 +3,7 @@
 
 # https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
 
-immutable SymbolicField <: Field
+immutable SymbolicField <: AbstractField
     name :: ASCIIString
 end
 

src/types.jl

@@ -25,25 +25,45 @@ type IntegrationPoint
 end
 
 function IntegrationPoint(xi, weight)
-    IntegrationPoint(xi, weight, Dict())
+    return IntegrationPoint(xi, weight, FieldSet())
 end
 
 function Base.convert(::Type{Number}, ip::IntegrationPoint)
     return ip.xi
 end
 
-function Base.call(basis::Basis, ip::IntegrationPoint)
-    return basis(ip.xi)
+function Base.call(field::CVTI, ip::IntegrationPoint)
+    return call(field, ip.xi)
 end
 
-function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint)
-    return call(basis, increment, ip.xi)
+function Base.call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
+    call(basis, field, ip.xi)
 end
 
-function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
-    return call(basis, increment, ip.xi, Val{:grad})
+function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+    call(basis, field, ip.xi, Val{:grad})
 end
 
-function Base.call(basis::Basis, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
-    return call(basis, geometry, field, ip.xi, Val{:grad})
+function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint)
+    call(basis, field, ip.xi)
 end
+
+function Base.call(basis::CVTI, geometry::DVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+    call(basis, geometry, field, ip.xi, Val{:grad})
+end
+
+#function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint)
+#    return call(basis, increment, ip.xi)
+#end
+#function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+#    return call(basis, increment, ip.xi, Val{:grad})
+#end
+#function Base.call(basis::Basis, field::Field, ip::IntegrationPoint, ::Type{Val{:grad}})
+#    return call(basis, field, ip.xi, Val{:grad})
+#end
+#function Base.call(basis::Basis, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+#    return call(basis, geometry, field, ip.xi, Val{:grad})
+#end
+#function Base.call(basis::Basis, field::Field, ip::IntegrationPoint)
+#    return call(basis, field, ip.xi)
+#end