updated developers guide + tests

src/basis.jl

@@ -45,12 +45,14 @@ function Base.call(field::DiscreteField, time::Number,
 
     # special cases, only 1 timestep defined or time = -Inf -> return first ts
     if (length(field) == 1) || (time == -Inf)
-        return field[1][end]
+        #return field[1][end]
+        return first(field)
     end
 
     # special case, time = +Inf -> return last ts
     if time == +Inf
-        return field[end][end]
+        #return field[end][end]
+        return last(field)
     end
 
     # very likely we are always near some defined timestep, usually field

src/dirichlet.jl

@@ -35,7 +35,7 @@ function DBC2D2(element::Seg2)
         IntegrationPoint([-sqrt(1/3)], 1.0),
         IntegrationPoint([+sqrt(1/3)], 1.0)]
     if !haskey(element, "reaction force")
-        element["reaction force"] = FieldSet()
+        element["reaction force"] = zeros(1, 2)
     end
     DBC2D2(element, integration_points)
 end

src/elements.jl

@@ -42,9 +42,9 @@ function test_element(element_type)
     end
 
     # try to interpolate some scalar field
-    element["field1"] = Field(0.0, collect(1:n))
+    element["field1"] = Field(collect(1:n))
     # TODO: how to parametrize this?
-    element["geometry"] = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    element["geometry"] = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
 
     # evaluate basis functions at middle point of element
     basis = get_basis(element)
@@ -55,7 +55,7 @@ function test_element(element_type)
     val2 = basis("field1", mid, 0.0)
     info("field val at $mid: $val2")
     val3 = dbasis(mid, 0.0)
-    info("derivative of basis at $mid: $val3")
+    info("derivative of basis at $mid:\n$val3")
     val4 = dbasis("field1", mid, 0.0)
     info("field val at $mid: $val4")
 
@@ -64,135 +64,130 @@ end
 
 """ Get FieldSet from element. """
 function Base.getindex(element::Element, field_name)
-    element.fields[field_name]
+    return element.fields[field_name]
 end
 
-"""Add new FieldSet to element.
+"""Add new Field to element.
 
 Examples
 --------
->>> element["geometry"] = [1, 2, 3, 4]
-JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
+>>> element["temperature"] = [1, 2, 3, 4]
+>>> 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)
-    #element.fields[field_name] = field_data
     setindex!(element.fields, field_data, field_name)
 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)
-    el.connectivity
+    return el.connectivity
 end
 
 abstract AbstractFunctionSpace
 
 type FunctionSpace <: AbstractFunctionSpace
-    element :: Element
+    basis :: Basis
+    fields :: FieldSet
 end
 
 type GradientFunctionSpace <: AbstractFunctionSpace
-    element :: Element
-end
-
-type MixedFunctionSpace <: AbstractFunctionSpace
-    element1 :: Element
-    element2 :: Element
+    basis :: Basis
+    fields :: FieldSet
 end
 
 function get_basis(element::Element)
-    return FunctionSpace(element)
+    return FunctionSpace(element.basis, element.fields)
 end
 
 function get_dbasis(element::Element)
-    return GradientFunctionSpace(element)
+    return GradientFunctionSpace(element.basis, element.fields)
 end
 
 function grad(u::FunctionSpace)
-    return GradientFunctionSpace(u.element)
-end
-
-""" Evaluate field on element function space. """
-function call(u::FunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
-    f = !isa(variation, Void) ? variation : u.element[field_name](t)
-    if length(f) == 1
-        return f.data[1]
-    end
-    h = u.element.basis.basis(xi)
-    #@debug("vec(h) = $(vec(h)), size(h) = $(size(vec(h)))")
-    #@debug("f = $f, size(f) = $(size(f))")
-    #return dot(vec(h), f)
-    return sum(vec(h).*f)
+    return GradientFunctionSpace(u.basis, u.fields)
 end
 
 """ If basis is called without a field, return basis functions evaluated at that point. """
-function call(u::FunctionSpace, xi::Vector, t::Number=Inf)
-    return u.element.basis.basis(xi)
-end
-
-""" Evaluate gradient of field on element function space. """
-function call(gradu::GradientFunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
-    f = !isa(variation, Void) ? variation : gradu.element[field_name](t)
-    X = gradu.element["geometry"](t)
-    dN = gradu.element.basis.dbasisdxi(xi)
-    J = sum([dN[:,i]*X[i]' for i=1:length(X)])
-    grad = inv(J)*dN
-    gradf = sum([grad[:,i]*f[i]' for i=1:length(f)])'
-    return gradf
+function call(u::FunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0)
+    return u.basis(xi)
 end
 
 """ If gradient of basis is called without a field, return "empty" gradient evaluated at that point. """
-function call(gradu::GradientFunctionSpace, xi::Vector, t::Number=Inf)
-    X = gradu.element["geometry"](t)
-    dN = gradu.element.basis.dbasisdxi(xi)
-    J = sum([dN[:,i]*X[i]' for i=1:length(X)])
-    grad = inv(J)*dN
-    return grad
+function call(gradu::GradientFunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0)
+    geometry = gradu.fields["geometry"](t)
+    gradu.basis(geometry, xi, Val{:grad})
 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
+
+""" Evaluate gradient of field on element function space. """
+function call(gradu::GradientFunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
+    field = !isa(variation, Void) ? variation : gradu.fields[field_name](t)
+    geometry = gradu.fields["geometry"](t)
+    gradu.basis(geometry, field, xi, Val{:grad})
+end
+
+
 # on-line functions to get api more easy to use, ip -> xi.ip
-call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
-call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
+#call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
+#call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
 # i think these will be the most called functions.
-call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation)
-call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation)
-call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
-call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
+#call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+#call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+#call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
+#call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
 
 """ Return a field from function space. """
-function get_field(u::FunctionSpace, field_name, time=Inf)
-    return u.element[field_name](time)
+function get_field(u::FunctionSpace, field_name, time::Number=0.0)
+    return u.fields[field_name](time)
 end
 
 """ Return a field from function space. """
-function get_field(u::FunctionSpace, field_name, time=Inf, variation=nothing)
-    return !isa(variation, Void) ? variation : u.element[field_name](time)
+function get_field(u::FunctionSpace, field_name, time::Number=0.0, variation=nothing)
+    return !isa(variation, Void) ? variation : u.fields[field_name](time)
 end
 
-""" Return a fieldset from function space. """
+""" Return a field from function space. """
 function get_fieldset(u::FunctionSpace, field_name)
-    return u.element[field_name]
+    return u.fields[field_name]
 end
 
 """ Get a determinant of element in point ξ. """
-function LinAlg.det(u::FunctionSpace, xi::Vector, t::Number=Inf)
-    X = u.element["geometry"](t)
-    dN = u.element.basis.dbasisdxi(xi)
+function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
+    X = u.fields["geometry"](time)
+    dN = u.basis.dbasisdxi(xi)
     J = sum([dN[:,i]*X[i]' for i=1:length(X)])
     m, n = size(J)
     return m == n ? det(J) : norm(J)
 end
-function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf)
-    LinAlg.det(u, ip.xi, t)
+function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, time::Number=0.0)
+    LinAlg.det(u, ip.xi, time)
 end
 function LinAlg.det(u::FunctionSpace)
     return (args...) -> det(u, args...)
 end
+#Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
+#Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
+#Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)
+#Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...)
 
-Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
-Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
-Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)
-Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...)
-
-""" Check does fieldset exist. """
+""" Check does field exist. """
 function Base.haskey(element::Element, what)
     haskey(element.fields, what)
 end

src/equations.jl

@@ -57,10 +57,10 @@ function initialize_local_assembly!(assembly::LocalAssembly, equation::Equation)
 end
 
 has_mass_matrix(equation::Equation) = false
-function get_mass_matrix(equation::Equation, ip, time=Inf, problem=nothing)
+function get_mass_matrix(equation::Equation, ip, time=0.0, problem=nothing)
     get_mass_matrix(equation, ip, time)
 end
-function get_mass_matrix(equation::Equation, ip, time=Inf)
+function get_mass_matrix(equation::Equation, ip, time=0.0)
     get_mass_matrix(equation, ip)
 end
 function get_mass_matrix(equation::Equation, ip)
@@ -68,10 +68,10 @@ function get_mass_matrix(equation::Equation, ip)
 end
 
 has_stiffness_matrix(equation::Equation) = false
-function get_stiffness_matrix(equation::Equation, ip, time=Inf, problem=nothing)
+function get_stiffness_matrix(equation::Equation, ip, time=0.0, problem=nothing)
     get_stiffness_matrix(equation, ip, time)
 end
-function get_stiffness_matrix(equation::Equation, ip, time=Inf)
+function get_stiffness_matrix(equation::Equation, ip, time=0.0)
     get_stiffness_matrix(equation, ip)
 end
 function get_stiffness_matrix(equation::Equation, ip)
@@ -79,10 +79,10 @@ function get_stiffness_matrix(equation::Equation, ip)
 end
 
 has_force_vector(equation::Equation) = false
-function get_force_vector(equation::Equation, ip, time=Inf, problem=nothing)
+function get_force_vector(equation::Equation, ip, time=0.0, problem=nothing)
     get_force_vector(equation, ip, time)
 end
-function get_force_vector(equation::Equation, ip, time=Inf)
+function get_force_vector(equation::Equation, ip, time=0.0)
     get_force_vector(equation, ip)
 end
 function get_force_vector(equation::Equation, ip)
@@ -90,10 +90,10 @@ function get_force_vector(equation::Equation, ip)
 end
 
 has_residual_vector(equation::Equation) = false
-function get_residual_vector(equation::Equation, ip, time=Inf, problem=nothing)
+function get_residual_vector(equation::Equation, ip, time=0.0, problem=nothing)
     get_residual_vector(equation, ip, time)
 end
-function get_residual_vector(equation::Equation, ip, time=Inf)
+function get_residual_vector(equation::Equation, ip, time=0.0)
     get_residual_vector(equation, ip)
 end
 function get_residual_vector(equation::Equation, ip)
@@ -101,10 +101,10 @@ function get_residual_vector(equation::Equation, ip)
 end
 
 has_potential_energy(equation::Equation) = false
-function get_potential_energy(equation::Equation, ip, time=Inf, problem=nothing)
+function get_potential_energy(equation::Equation, ip, time=0.0, problem=nothing)
     get_potential_energy(equation, ip, time)
 end
-function get_potential_energy(equation::Equation, ip, time=Inf)
+function get_potential_energy(equation::Equation, ip, time=0.0)
     get_potential_energy(equation, ip)
 end
 function get_potential_energy(equation::Equation, ip)
@@ -117,7 +117,7 @@ get_integration_points(equation::Equation) = equation.integration_points
 
 """ Return a local assembly for element. """
 function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
-                                   unknown_field_name::ASCIIString, time::Number=Inf,
+                                   unknown_field_name::ASCIIString, time::Number=0.0,
                                    problem=nothing)
 
     initialize_local_assembly!(assembly, equation) # zero all
@@ -174,6 +174,7 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
         assembly.stiffness_matrix += hessian
         assembly.force_vector -= ForwardDiff.gradient(allresults)  # <--- minus explained in tutorial
         assembly.potential_energy = ForwardDiff.value(allresults)
+        #info("potential energy of system: $(assembly.potential_energy)")
     end
 
     # 3. virtual work form - user has defined residual vector δW_int(u,δu) + δW_ext(u,δu) = 0 ∀ v

src/fields.jl

@@ -214,6 +214,15 @@ function Base.push!(field::DefaultDiscreteField, timestep::TimeStep)
     push!(field.timesteps, timestep)
 end
 
+function Base.push!(field::DefaultDiscreteField, data::Union{Vector, Matrix})
+    push!(field[end], Increment(data))
+end
+
+function Base.push!(field::DefaultDiscreteField, data::Pair)
+    ts = TimeStep(data[1], Increment(data[2]))
+    push!(field, ts)
+end
+
 """Quickly create fields.
 
 Examples
@@ -287,3 +296,7 @@ function Base.convert(::Type{ContinuousField}, data::Function)
     return convert(DefaultContinuousField, data)
 end
 
+function Base.length(::Field)
+    return 1
+end
+

src/integrate.jl

@@ -16,3 +16,48 @@ function get_default_integration_points(element::Seg2)
         IntegrationPoint([0.0], 2.0)
     ]
 end
+
+function line3()
+    [
+        IntegrationPoint([0.0], 8/9),
+        IntegrationPoint([-sqrt(3/5)], 5/9),
+        IntegrationPoint([+sqrt(3/5)], 5/9)
+    ]
+end
+
+function line5()
+    [
+        IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
+        IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
+        IntegrationPoint([0.0], 128/225),
+        IntegrationPoint([ 1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
+        IntegrationPoint([ 1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900)
+    ]
+end
+
+    #integration_points = [
+    #    IntegrationPoint([ 0.0000000000000000], 0.5688888888888889),
+    #    IntegrationPoint([-0.5384693101056831], 0.4786286704993665),
+    #    IntegrationPoint([ 0.5384693101056831], 0.4786286704993665),
+    #    IntegrationPoint([-0.9061798459386640], 0.2369268850561891),
+    #    IntegrationPoint([ 0.9061798459386640], 0.2369268850561891)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
+    #    IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
+    #    IntegrationPoint([+sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
+    #    IntegrationPoint([-sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([0.0], 8/9),
+    #    IntegrationPoint([-sqrt(3/5)], 5/9),
+    #    IntegrationPoint([+sqrt(3/5)], 5/9)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([-sqrt(1/3)], 1)
+    #    IntegrationPoint([+sqrt(1/3)], 1)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([0.0], 2)
+    #]
+

src/solvers.jl

@@ -10,10 +10,10 @@ Solve field equations for single element with some dofs fixed. This can be used
 to test nonlinear element formulations.
 """
 function solve!(equation::Equation, unknown_field_name::ASCIIString,
-                 free_dofs::Array{Int, 1}, time::Number=Inf; 
+                 free_dofs::Array{Int, 1}, time::Number=0.0; 
                  max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
     element = get_element(equation)
-    x0 = element[unknown_field_name](-Inf)
+    x0 = element[unknown_field_name](0.0)
     x = zeros(prod(size(equation)))
     dx = fill!(similar(x), 0.0)
     la = initialize_local_assembly()
@@ -27,15 +27,10 @@ function solve!(equation::Equation, unknown_field_name::ASCIIString,
         end
         dx[free_dofs] = A \ b
         x += dx
-        new_field = similar(x0, x)
-        new_field.time = time
-        new_field.increment = i
-        push!(element[unknown_field_name], new_field)
-        if norm(dx) < tolerance
-            return
-        end
+        push!(element[unknown_field_name], reshape(x, size(equation)))
+        norm(dx) < tolerance && return
     end
-    Logging.err("Did not converge in $max_iterations iterations")
+    error("Did not converge in $max_iterations iterations")
 end
 
 """
@@ -44,15 +39,18 @@ 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::Array{Int, 1}, time::Number=Inf;
+function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=1.0;
                 max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+    info("start solver")
     ga = initialize_global_assembly(problem)
     x = zeros(ga.ndofs)
     dx = fill!(similar(x), 0.0)
     field_name = get_unknown_field_name(problem)
     dim = get_unknown_field_dimension(problem)
     for i=1:max_iterations
+        info("calculate global assembly")
         calculate_global_assembly!(ga, problem)
+        info("done")
         A = ga.stiffness_matrix[free_dofs, free_dofs]
         b = ga.force_vector[free_dofs]
         if dump_matrices
@@ -60,20 +58,17 @@ function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=Inf;
             dump(full(b)')
         end
         dx[free_dofs] = lufact(A) \ full(b)
+        info("Difference in solution norm: $(norm(dx))")
         x += dx
         for equation in get_equations(problem)
             element = get_element(equation)
-            conn = get_connectivity(element)
-            gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
-            old_field = element[field_name](Inf)
-            new_field = similar(old_field, full(x[gdofs]))
-            push!(element[field_name][end], new_field)
-        end
-        if norm(dx) < tolerance
-            return
+            gdofs = get_gdofs(problem, equation)
+            data = reshape(full(x[gdofs]), size(equation))
+            push!(element[field_name], data)
         end
+        norm(dx) < tolerance && return
     end
-    Logging.err("Did not converge in $max_iterations iterations")
+    error("Did not converge in $max_iterations iterations")
 end
 
 """ Add new problem to solver. """

src/types.jl

@@ -31,3 +31,19 @@ end
 function Base.convert(::Type{Number}, ip::IntegrationPoint)
     return ip.xi
 end
+
+function Base.call(basis::Basis, ip::IntegrationPoint)
+    return basis(ip.xi)
+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, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+    return call(basis, geometry, field, ip.xi, Val{:grad})
+end