multiple dirichlet boundary conditions for vector valued functions. direct solver design.

src/JuliaFEM.jl

@@ -10,6 +10,10 @@ module JuliaFEM
 #@Logging.configure(level=DEBUG)
 #using Lexicon
 
+import Base: +, -, /, *, push!, convert, getindex, length, similar, call, vec, endof
+
+#importall Base
+
 """
 A very simple debugging macro. It prints debug message if environment variable
 JULIAFEM_DEBUG is found.
@@ -78,6 +82,7 @@ include("elasticity.jl")
 ### ASSEMBLY + SOLVE ###
 include("assembly.jl")
 include("solvers.jl")
+include("directsolver.jl") # parallel sparse direct solver for non-lniear problems
 
 # PRE AND POSTPROCESS
 include("xdmf.jl")

src/directsolver.jl

@@ -0,0 +1,114 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+## Direct solver
+
+type DirectSolver <: Solver
+    field_problems :: Vector{FieldProblem}
+    boundary_problems :: Vector{BoundaryProblem}
+    nonlinear_problem :: Bool
+    max_iterations :: Int64
+    tol :: Float64
+end
+
+function push!(solver::DirectSolver, problem::FieldProblem)
+    push!(solver.field_problems, problem)
+end
+
+function push!(solver::DirectSolver, problem::BoundaryProblem)
+    push!(solver.boundary_problems, problem)
+end
+
+""" Default initializer. """
+function DirectSolver()
+    DirectSolver([], [], true, 10, 1.0e-6)
+end
+
+""" Call solver to solve a set of problems. """
+function call(solver::DirectSolver, time::Number=0.0)
+    @assert length(solver.field_problems) == 1
+    @assert length(solver.boundary_problems) == 1
+    @assert solver.nonlinear_problem == true
+
+    problem1 = solver.field_problems[1]
+    problem2 = solver.boundary_problems[1]
+
+    x = zeros(3)
+    dx = zeros(3)
+    dims = nothing
+
+    for iter=1:solver.max_iterations
+        tic()
+        info("Starting iteration $iter")
+        assembly1 = Assembly()
+        assemble!(assembly1, problem1, time)
+        assembly2 = Assembly()
+        assemble!(assembly2, problem2, time)
+
+        A1 = sparse(assembly1.stiffness_matrix)
+        dims = size(A1)
+        b1 = sparse(assembly1.force_vector, dims[1], 1)
+        A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2])
+        b2 = sparse(assembly2.force_vector, dims[1], 1)
+    
+        # create a saddle point problem
+        A = [A1 A2; A2' zeros(A2)]
+        b = [b1; b2]
+
+        if length(b) != length(x)
+            info("iter $iter: resizing solution vector")
+            resize!(x, length(b))
+            resize!(dx, length(b))
+            fill!(x, 0.0)
+            fill!(dx, 0.0)
+        end
+
+        # solve problem, update solution vector
+        nz = unique(rowvals(A))  # take only non-zero rows
+        dx[nz] = lufact(A[nz,nz]) \ full(b[nz])
+        x += dx
+
+        # get "problem-wise" solution vectors
+        x1 = x[1:dims[1]]
+        x2 = x[dims[1]+1:end]
+
+        # update field for elements in problem 1
+        for equation in get_equations(problem1)
+            element = get_element(equation)
+            field_name = get_unknown_field_name(problem1)
+            gdofs = get_gdofs(problem1, equation)
+            local_sol = vec(full(x1[gdofs]))
+            eqsize = size(equation)
+            if eqsize[1] != 1
+                local_sol = reshape(local_sol, eqsize)
+            end
+            #info("problem1: pushing to $field_name")
+            push!(element[field_name], time => local_sol)
+        end
+
+        # update field for elements in problem 2 (Dirichlet boundary)
+        for equation in get_equations(problem2)
+            element = get_element(equation)
+            field_name = "reaction force" #get_unknown_field_name(problem2)
+            gdofs = get_gdofs(problem2, equation)
+            local_sol = vec(full(x1[gdofs]))
+            eqsize = size(equation)
+            if eqsize[1] != 1
+                local_sol = reshape(local_sol, eqsize)
+            end
+            #info("problem2: pushing to $field_name")
+            push!(element[field_name], time => local_sol)
+        end
+        
+        if norm(dx[1:dims[1]]) < solver.tol
+            return (iter, true)
+        end
+
+        info("Iteration took $(toq()) seconds")
+    end
+
+    info("Warning: did not coverge in $(solver.max_iterations) iterations!")
+    return (solver.max_iterations, false)
+
+end
+

src/dirichlet.jl

@@ -3,11 +3,7 @@
 
 # Dirichlet boundary conditions in weak form
 
-abstract DirichletEquation <: Equation
-
-function get_unknown_field_name(equation::DirichletEquation)
-    return "reaction force"
-end
+abstract DirichletEquation <: BoundaryEquation
 
 ### Dirichlet problem + equations
 
@@ -15,8 +11,6 @@ type DirichletProblem <: BoundaryProblem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
     equations :: Vector{DirichletEquation}
-#   element_mapping :: Dict{DataType, DataType}
-    field_value :: Function
 end
 
 """ Initialize new Dirichlet boundary condition.
@@ -24,27 +18,14 @@ end
 Parameters
 ----------
 dimension
-    dimension of unknown field
-field_value
-    boundary function
+    dimension of unknown field (scalar, vector, ...)
 
 Examples
 --------
 
-Create u(X) = 0.0 boundary condition for three-dimensional elasticity problem:
-
->>> u(X) = [0.0, 0.0, 0.0]
->>> bc = DirichletProblem(3, u)
 """
-function DirichletProblem(dimension::Int=1, field_value::Function=(X)->[0.0,0.0,0.0], equations=[])
-#    element_mapping = nothing
-#    if dimension == 1
-#        element_mapping = Dict(
-#            Seg2 => DBC2D2
-#            )
-#    end
-    DirichletProblem("reaction force", dimension, equations, field_value)
-#    DirichletProblem("reaction force", dimension, [], element_mapping, field_value)
+function DirichletProblem(unknown_field_name::ASCIIString, dimension::Int=1)
+    DirichletProblem(unknown_field_name, dimension, [])
 end
 
 """ Dirichlet boundary condition element for 2 node line segment """
@@ -57,31 +38,45 @@ function Base.size(equation::DBC2D2)
     return (1, 2)
 end
 
-#function DBC2D2(element::Seg2)
 function Base.convert(::Type{DirichletEquation}, element::Seg2)
     integration_points = line3()
-    haskey(element, "reaction force") || (element["reaction force"] = zeros(1, 2))
+    haskey(element, "reaction force") || (element["reaction force"] = 0.0 => zeros(2))
     DBC2D2(element, integration_points)
 end
 
 
 function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number=0.0, problem=nothing)
-    gdofs = get_gdofs(equation)
-#   info("gdofs = $gdofs")
+    isa(problem, Void) && error("Dicihlet boundary condition needs problem defined")
+    field_dim = problem.unknown_field_dimension
+    field_name = problem.unknown_field_name
     element = get_element(equation)
+    gdofs = get_gdofs(element, field_dim)
     basis = get_basis(element)
     detJ = det(basis)
     for ip in get_integration_points(equation)
         w = ip.weight * detJ(ip)
         N = basis(ip, time)
-        add!(assembly.stiffness_matrix, gdofs, gdofs, w*N'*N)
-#       info("added $(w*N'*N)")
-#        if !isa(problem, Void)
-#            X = basis("geometry", ip, time)
-#            u = problem.field_value(X)[1:length(gdofs)]
-#            add!(assembly.force_vector, gdofs, w*N'*u)
-#        end
+        A = w*N'*N
+
+        if haskey(element, field_name)
+            # add all dimensions at once
+            for i=1:field_dim
+                g = element(field_name, ip, time)
+                ldofs = gdofs[i:field_dim:end]
+                add!(assembly.stiffness_matrix, ldofs, ldofs, A)
+                add!(assembly.force_vector, ldofs, w*g*N')
+            end
+        end
+  
+        for i=1:field_dim
+            if haskey(element, field_name*" $i")
+                # add single component
+                g = element(field_name*" $i", ip, time)
+                ldofs = gdofs[i:field_dim:end]
+                add!(assembly.stiffness_matrix, ldofs, ldofs, A)
+                add!(assembly.force_vector, ldofs, w*g*N')
+            end
+        end
     end
-#   info("assembly done")
 end
 

src/elasticity.jl

@@ -3,7 +3,7 @@
 
 # Elasticity problems
 
-abstract ElasticityProblem <: Problem
+abstract ElasticityProblem <: FieldProblem
 abstract ElasticityEquation <: Equation
 
 function get_unknown_field_name(equation::ElasticityEquation)

src/elements.jl

@@ -67,6 +67,10 @@ function Base.getindex(element::Element, field_name)
     return element.fields[field_name]
 end
 
+function get_integration_points(element)
+    return get_default_integration_points(element)
+end
+
 """Add new Field to element.
 
 Examples
@@ -78,6 +82,9 @@ Examples
 function Base.setindex!(element::Element, data, name::ASCIIString)
     element.fields[name] = Field(data)
 end
+function Base.setindex!(element::Element, field::Field, name::ASCIIString)
+    element.fields[name] = field
+end
 function Base.setindex!(element::Element, data::Tuple, name::ASCIIString)
     element.fields[name] = Field(data...)
 end
@@ -143,6 +150,39 @@ function call(gradu::GradientFunctionSpace, field_name, xi::Union{Vector, Integr
     gradu.basis(geometry, field, xi, Val{:grad})
 end
 
+typealias VecOrIP Union{Vector, IntegrationPoint}
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number)
+    return element.basis(element[field_name](time), xi)
+end
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}})
+    return element.basis(element["geometry"](time), element[field_name](time), xi, Val{:grad})
+end
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP)
+    return element.basis(element[field_name], xi)
+end
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}})
+    return element.basis(element["geometry"], element[field_name], xi, Val{:grad})
+end
+
+function Base.call(element::Element, field_name::ASCIIString, time::Number)
+    return element[field_name](time)
+end
+
+function Base.call(element::Element, xi::VecOrIP)
+    element.basis(xi)
+end
+
+function Base.call(element::Element, xi::VecOrIP, ::Type{Val{:grad}})
+    element.basis(element["geometry"], xi, Val{:grad})
+end
+
+function Base.call(element::Element, field_name::ASCIIString)
+    return element[field_name]
+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)
@@ -176,12 +216,19 @@ function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
     m, n = size(J)
     return m == n ? det(J) : norm(J)
 end
+
 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
+
+function LinAlg.det(element::Element)
+    return det(get_basis(element))
+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...)

src/equations.jl

@@ -4,6 +4,8 @@
 # Functions to handle element level things -- integration, assembly, ...
 
 abstract Equation
+abstract FieldEquation <: Equation
+abstract BoundaryEquation <: Equation
 
 type Assembly
     mass_matrix :: SparseMatrixIJV

src/fields2.jl

@@ -12,6 +12,7 @@ abstract Variable <: AbstractField
 abstract TimeVariant <: AbstractField
 abstract TimeInvariant <: AbstractField
 
+
 type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}}
     data
 end
@@ -81,6 +82,10 @@ typealias TimeVariantField   Union{DCTV, DVTV, CCTV, CVTV}
 
 ### Convenient functions to create fields
 
+#function Base.convert(::Type{Field}, data)
+#    return Field(data)
+#end
+
 function Field(data)
     return DCTI(data)
 end
@@ -101,6 +106,20 @@ function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...)
     return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
+function Field(func::Function)
+    if method_exists(func, Tuple{})
+        return CCTI(func)
+    elseif method_exists(func, Tuple{Float64})
+        return CCTV(func)
+    elseif method_exists(func, Tuple{Vector})
+        return CVTI(func)
+    elseif method_exists(func, Tuple{Vector, Number})
+        return CVTV(func)
+    else
+        error("no proper definition found for function: check methods.")
+    end
+end
+
 function CVTI(basis::Function, dbasis::Function)
     return CVTI(Basis(basis, dbasis))
 end
@@ -203,28 +222,55 @@ function Base.convert(::Type{Basis}, field::CVTI)
     return field.data
 end
 
+function Base.call(field::CCTV, time::Number)
+    return field.data(time)
+end
+
 ### Interpolation 
 
+""" Interpolate time-invariant field in time direction. """
 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(field::CVTI, time::Float64)
+    return field.data()
+end
+function Base.call(field::CCTI, time::Float64)
+    return field.data()
+end
 
+""" Interpolate time-variant field in time direction. """
+function Base.call(field::DCTV, time::Float64)
+    for i=reverse(1:length(field))
+        if isapprox(field[i].time, time)
+            return DCTI(field[i].data)
+        end
+    end
+    info(field.data)
+    info(time)
+    error("interpolate DCTV: not implemented yet")
+end
+
+function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
+    for i=reverse(1:length(field))
+        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
+
+""" Interpolate constant field in spatial dimension. """
 function Base.call(basis::CVTI, field::DCTI, xi::Vector)
-    # try to interpolate constant value -> no effect
-    return field
+    return field.data
 end
 
-#function Base.call(basis::Basis, field::DCTI, xi::Vector)
-    # calling constant field with basis -> no effect
-#    return field
-#end
-
+""" Interpolate variable field in spatial dimension. """
 function Base.call(basis::CVTI, values::DVTI, xi::Vector)
     N = basis(xi)
     return sum([N[i]*values[i] for i=1:length(N)])
@@ -244,27 +290,8 @@ function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type
     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")
+function Base.call(basis::CVTI, xi::Vector, time::Number)
+    call(basis, xi)
 end
 
 ### FIELDSET ###

src/solvers.jl

@@ -9,7 +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; 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, callback=nothing)
     unknown_field_name = get_unknown_field_name(equation)
     element = get_element(equation)
     x0 = element[unknown_field_name](0.0)
@@ -31,6 +31,9 @@ function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_it
         data = eqsize[1] != 1 ? reshape(x, eqsize) : x
         push!(element[unknown_field_name], time => data)
         norm(dx) < tolerance && return
+        if !isa(callback, Void)
+            callback(x)
+        end
     end
     error("Did not converge in $max_iterations iterations")
 end
@@ -41,7 +44,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::Float64; 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, callback=nothing)
     info("start solver")
     assembly = Assembly()
 #   x = zeros(ga.ndofs)
@@ -66,6 +69,9 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite
         dx[free_dofs] = lufact(A[free_dofs,free_dofs]) \ full(b)[free_dofs]
         info("Difference in solution norm: $(norm(dx))")
         x += dx
+        if !(isa(callback, Void))
+            callback(x)
+        end
         for equation in get_equations(problem)
             element = get_element(equation)
             gdofs = get_gdofs(equation)
@@ -81,11 +87,6 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite
     error("Did not converge in $max_iterations iterations")
 end
 
-""" Add new problem to solver. """
-function add_problem!(solver::Solver, problem::Problem)
-    push!(solver.problems, problem)
-end
-
 function Base.push!(solver::Solver, problem::Problem)
     push!(solver.problems, problem)
 end
@@ -126,9 +127,10 @@ function call(solver::SimpleSolver, time::Number=0.0)
 
 #   info("Creating sparse matrices")
     A1 = sparse(assembly1.stiffness_matrix)
-    b1 = sparse(assembly1.force_vector, size(A1, 1), 1)
-    A2 = sparse(assembly2.stiffness_matrix)
-    b2 = sparse(assembly2.force_vector, size(A2, 1), 1)
+    dims = size(A1)
+    b1 = sparse(assembly1.force_vector, dims[1], 1)
+    A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2])
+    b2 = sparse(assembly2.force_vector, dims[1], 1)
     
     # create a saddle point problem
     A = [A1 A2; A2' zeros(A2)]
@@ -149,16 +151,28 @@ function call(solver::SimpleSolver, time::Number=0.0)
         field_name = get_unknown_field_name(problem1)
         gdofs = get_gdofs(problem1, equation)
         local_sol = vec(full(x1[gdofs]))
+        eqsize = size(equation)
+        if eqsize[1] != 1
+            local_sol = reshape(local_sol, eqsize)
+        end
+        #info("problem1: pushing to $field_name")
         push!(element[field_name], time => local_sol)
     end
 
     # update field for elements in problem 2 (Dirichlet boundary)
     for equation in get_equations(problem2)
         element = get_element(equation)
-        field_name = get_unknown_field_name(problem2)
+        field_name = "reaction force" #get_unknown_field_name(problem2)
         gdofs = get_gdofs(problem2, equation)
         local_sol = vec(full(x1[gdofs]))
+        eqsize = size(equation)
+        if eqsize[1] != 1
+            local_sol = reshape(local_sol, eqsize)
+        end
+        #info("problem2: pushing to $field_name")
         push!(element[field_name], time => local_sol)
     end
+
+    return norm(x1)
 end
 

test/test_dirichlet.jl

@@ -1,18 +1,60 @@
-
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-module TestAutoDiffWeakForm
+module TestDirichletBoundaryCondition
 
 using JuliaFEM.Test
 using JuliaFEM
-using JuliaFEM: Seg2, DirichletProblem
+using JuliaFEM: Seg2, DirichletProblem, Assembly, assemble!
 
-function test_dirichlet_problem()
-    element = Seg2([3, 4])
+function test_dirichlet_problem_1_dim()
+    element = Seg2([1, 2])
     element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
-    problem = DirichletProblem(1)
+    element["temperature"] = 0.0
+    problem = DirichletProblem("temperature", 1)
     push!(problem, element)
+    assembly = Assembly()
+    assemble!(assembly, problem)
+    A = full(assembly.stiffness_matrix)
+    b = full(assembly.force_vector)
+    @test isapprox(A, 1/6*[2 1; 1 2])
+    @test isapprox(b, [0.0, 0.0])
+end
+
+function test_dirichlet_problem_2_dim()
+    element = Seg2([1, 2])
+    element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
+    element["displacement"] = 0.0
+    problem = DirichletProblem("displacement", 2)
+    push!(problem, element)
+    assembly = Assembly()
+    assemble!(assembly, problem)
+    A = full(assembly.stiffness_matrix)
+    b = full(assembly.force_vector)
+    A_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2]
+    @test isapprox(A, A_expected)
+    @test isapprox(b, [0.0, 0.0, 0.0, 0.0])
+end
+
+function test_dirichlet_problem_2_dim_single_dof_fixed()
+    element = Seg2([1, 2])
+    element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
+    element["displacement 2"] = 0.0
+    problem = DirichletProblem("displacement", 2)
+    push!(problem, element)
+    assembly = Assembly()
+    assemble!(assembly, problem)
+    A = full(assembly.stiffness_matrix)
+    b = full(assembly.force_vector)
+    info(b)
+    info("A = \n$A")
+    A_expected = 1/6*[
+        0 0 0 0
+        0 2 0 1
+        0 0 0 0
+        0 1 0 2]
+    @test isapprox(A, A_expected)
+    @test isapprox(b, [0.0, 0.0, 0.0, 0.0])
 end
 
 end

test/test_elasticity.jl

@@ -25,27 +25,28 @@ function test_elasticity_volume_load()
 end
 
 function test_elasticity_surface_load()
-    N = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
 
-    element1 = Quad4([1, 2, 3, 4])
-    element1["geometry"] = Vector[N[1], N[2], N[3], N[4]]
-    element1["youngs modulus"] = 500.0
-    element1["poissons ratio"] = 0.3
+    element1 = Quad4([1, 2, 4, 3])
+    element1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    element1["youngs modulus"] = 900.0
+    element1["poissons ratio"] = 0.25
     element2 = Seg2([3, 4])
     element2["geometry"] = Vector[N[3], N[4]]
-    element2["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]]
+    element2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
 
-    free_dofs = [3, 4, 5, 6]
+    #free_dofs = [3, 5, 6, 8]
+    free_dofs = [3, 6, 7, 8]
     problem = PlaneStressElasticityProblem()
     push!(problem, element1)
     push!(problem, element2)
-    solve!(problem, free_dofs, 1.0; max_iterations=10)
-    disp = get_basis(element1)("displacement", [1.0, 1.0], 1.0)[2]
+    solve!(problem, free_dofs, 0.0; max_iterations=10)
+    #disp = get_basis(element1)("displacement", [1.0, 1.0], 1.0)[2]
+    info(last(element1["displacement"]))
+    disp = element1("displacement", [1.0, 1.0], 0.0)
     info("displacement at tip: $disp")
     # verified using Code Aster.
-    @test isapprox(disp, -9.33106637611714)
+    @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
 end
 
-#test_elasticity_volume_load()
-
 end

test/test_elements.jl

@@ -5,7 +5,7 @@ module ElementTests
 
 using JuliaFEM.Test
 
-using JuliaFEM: Element, Basis, Field, FieldSet, FunctionSpace, test_element
+using JuliaFEM: Element, Field, FieldSet, test_element
 
 """ Prototype element
 
@@ -13,11 +13,11 @@ This should always pass test_element if everything is ok.
 """
 type MockElement <: Element
     connectivity :: Vector{Int}
-    basis :: Basis
+    basis :: Field
     fields :: FieldSet
 end
 
-function MockElement(connectivity)
+function MockElement(connectivity, fields...)
 
     h(xi) = 1/4*[
         (1-xi[1])*(1-xi[2])
@@ -29,95 +29,51 @@ function MockElement(connectivity)
         -(1-xi[2])    (1-xi[2])   (1+xi[2])  -(1+xi[2])
         -(1-xi[1])   -(1+xi[1])   (1+xi[1])   (1-xi[1])]
 
-    basis = Basis(h, dh)
-    MockElement(connectivity, basis, Dict())
+    MockElement(connectivity, Field(h, dh), FieldSet(fields...))
 end
 
 Base.size(element::Type{MockElement}) = (2, 4)
 
-"""test test_element against mock element"""
-function test_mockelement()
+""" Return test element with some fields. """
+function get_element()
+    el = MockElement([1, 2, 3, 4])
+    el["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
+    el["temperature"] = (
+        0.0 => [0.0, 0.0, 0.0, 0.0],
+        1.0 => [1.0, 2.0, 3.0, 4.0])
+    el["displacement"] = (
+        0.0 => Vector{Float64}[[0.0,0.0], [0.0, 0.0], [0.0,0.0], [0.0,0.0]],
+        1.0 => Vector{Float64}[[0.0,0.0], [1.0,-1.0], [2.0,3.0], [0.0,0.0]])
+    return el
+end
+
+function test_mock_element()
     test_element(MockElement)
 end
 
-""" test adding fieldsets and fields to element"""
 function test_add_fields_to_element()
-    el = MockElement([1, 2, 3, 4])
-    #geometry = Field([0.0, 0.0, 0.0, 0.0])
-    el["geometry"] = Field([0.0, 0.0, 0.0, 0.0])
-    @test el["geometry"][1].time == 0.0
-    @test last(el["geometry"]) == [0.0, 0.0, 0.0, 0.0]
-    el["geometry"] = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
-    @test last(el["geometry"])[3] == [1.0, 1.0]
-    el["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
-    @test last(el["geometry"])[3] == [1.0, 1.0]
-    el["geometry"] = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]'
-    @test last(el["geometry"])[3] == [1.0, 1.0]
-    el["geometry"] = (0.0, [0.0, 0.0, 0.0, 0.0]), (1.0, [1.0, 1.0, 1.0, 1.0])
-    field = el["geometry"]
-    @test length(field) == 2  # two time steps
-    el["boundary flux"] = (0.0, 0.0), (1.0, 6.0)
+    el = get_element()
+    info(el.fields)
 end
 
-function test_add_fields_to_element_2()
-    el = MockElement([1, 2, 3, 4])
-    el["data"] = (0.0 => [1, 2], 1.0 => [2, 3])
-    @test length(el["data"]) == 2
-    @test el["data"][1].time == 0.0
-    @test el["data"][2].time == 1.0
-    @test last(el["data"][1]) == [1, 2]
-    @test last(el["data"][2]) == [2, 3]
+function test_interpolate()
+    el = get_element()
+    @test isapprox(el("geometry", [0.0, 0.0]), [0.5, 0.5])
+    @test isapprox(el("geometry", [0.0, 0.0], 0.0), [0.5, 0.5])
+    @test isapprox(el([0.0, 0.0]), [0.25 0.25 0.25 0.25])
+    @test isapprox(el([0.0, 0.0], Val{:grad}), [-0.5 0.5 0.5 -0.5; -0.5 -0.5 0.5 0.5])
+    gradT = el("temperature", [0.0, 0.0], 1.0, Val{:grad})
+    info("gradT = $gradT")
+    X = [0.5, 0.5]
+    gradT_expected = [1-2*X[2] 3-2*X[1]]
+    info("gradT(expected) = $gradT_expected")
+    @test isapprox(gradT, gradT_expected)
+
+    @test isapprox(el("temperature", [0.0, 0.0], 0.5), 1/2*gradT_expected)
+
+    gradT = el("temperature", [0.0, 0.0], 0.5, Val{:grad})
+    info("gradT = $gradT")
+    @test isapprox(gradT, 1/2*gradT_expected)
 end
 
-function test_add_data_to_element_using_push()
-    el = MockElement([1, 2, 3, 4])
-    el["data"] = [0, 0, 0, 0]
-   
-    push!(el["data"], [1, 2, 3, 4])
-    @test length(el["data"]) == 1
-    @test length(el["data"][1]) == 2
-    @test el["data"][1].time == 0.0
-
-    push!(el["data"], 1.0 => [2, 3, 4, 5])  # creates new timestep at t=1.0
-    push!(el["data"], [3, 4, 5, 6])  # adds new increment data to last timestep
-    @test length(el["data"]) == 2
-    @test length(el["data"][2]) == 2
-    @test el["data"][2].time == 1.0
-    
-end
-
-#=
-facts("interpolation of fields in some function space") do
-
-    el = MockElement([1, 2, 3, 4])
-    fieldset1 = FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])
-    fieldset2 = FieldSet("constant scalar field", [Field(0.0, 1.0)])
-    fieldset3 = FieldSet("scalar field", [Field(0.0, [1.0, 2.0, 3.0, 4.0])])
-    fieldset4 = FieldSet("vector field 1", [Field(0.0, Vector[[1.0], [2.0], [3.0], [4.0]])])
-    fieldset5 = FieldSet("vector field 2", [Field(0.0, Vector[[1.0, 5.0], [2.0, 6.0], [3.0, 7.0], [4.0, 8.0]])])
-    fieldset6 = FieldSet("vector field 3", [Field(0.0, Vector[[1.0, 5.0, 9.0], [2.0, 6.0, 10.0], [3.0, 7.0, 11.0], [4.0, 8.0, 12.0]])])
-    fieldset7 = FieldSet("tensor field 1", [Field(0.0, Matrix[[1.0 5.0; 9.0 13.0], [2.0 6.0; 10.0 14.0], [3.0 7.0; 11.0 15.0], [4.0 8.0; 12.0 16.0]])])
-
-    element["geometry"] = fieldset1
-    element["constant scalar field"] = fieldset2
-    element["scalar field"] = fieldset3
-    element["vector field 1"] = fieldset4
-    element["vector field 2"] = fieldset5
-    element["vector field 3"] = fieldset6
-    element["tensor field 1"] = fieldset7
-
-    xi = [0.0, 0.0]
-    t = 0.0
-    u = FunctionSpace(element)
-    v = FunctionSpace(element)
-
-    @fact v("constant scalar field", xi, t) --> 1.0
-    @fact v("scalar field", xi, t) --> 1/4*(1+2+3+4)
-    @fact v("vector field 1", xi, t) --> [1/4*(1+2+3+4)]
-    @fact v("vector field 2", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8]
-    @fact v("vector field 3", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8, 9+10+11+12]
-    @fact v("tensor field 1", xi, t) --> 1/4*[1+2+3+4 5+6+7+8; 9+10+11+12 13+14+15+16]
-end
-=#
-
 end

test/test_solver.jl

@@ -6,7 +6,7 @@ module SolverTests
 using JuliaFEM.Test
 using JuliaFEM
 
-using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis, MortarElement, MortarProblem, DirectSolver, PlaneStressElasticityProblem
+using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis, MortarElement, MortarProblem, PlaneStressElasticityProblem, solve!, DirectSolver
 
 """ Define Problem 1:
 
@@ -38,7 +38,8 @@ end
 function get_boundaryproblem()
     el3 = Seg2([3, 4])
     el3["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
-    problem2 = DirichletProblem(1)
+    el3["temperature"] = 0.0
+    problem2 = DirichletProblem("temperature", 1)
     push!(problem2, el3)
     return problem2
 end
@@ -67,7 +68,7 @@ function test_simplesolver()
     @test isapprox(T, 100.0)
 end
 
-function test_direct_solver()
+function atest_direct_solver()
 
     N = Dict{Int, Vector}(
         1 => [0.0, 0.0],
@@ -87,7 +88,7 @@ function test_direct_solver()
 
    # volume elements
    e1 = Quad4([1, 2, 5, 4])
-   e1["geometry"] = Vector[N[1], N[2], N[3], N[4]]
+   e1["geometry"] = Vector[N[1], N[2], N[5], N[4]]
    e2 = Quad4([2, 3, 6, 5])
    e2["geometry"] = Vector[N[2], N[3], N[6], N[5]]
    e3 = Quad4([7, 8, 12, 11])
@@ -178,4 +179,56 @@ function test_direct_solver()
 
 end
 
+function test_solver_multiple_dirichlet_bc()
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
+
+    e1 = Quad4([1, 2, 4, 3])
+    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    e1["youngs modulus"] = 900.0
+    e1["poissons ratio"] = 0.25
+    b1 = Seg2([3, 4])
+    b1["geometry"] = Vector[N[3], N[4]]
+    b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
+
+    #free_dofs = [3, 5, 6, 8]
+    free_dofs = [3, 6, 7, 8]
+    problem = PlaneStressElasticityProblem()
+    push!(problem, e1)
+    push!(problem, b1)
+
+    # manually solve problem 1
+    #solve!(problem, free_dofs, 0.0; max_iterations=10)
+    #disp = e1("displacement", [1.0, 1.0], 0.0)
+    #info("displacement at tip: $disp")
+    #@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
+
+    # boundary elements for dirichlet dx=0
+    dx = Seg2([1, 3])
+    dx["geometry"] = Vector[N[1], N[3]]
+    dx["displacement 1"] = 0.0
+
+    # boundary elements for dirichlet dy=0
+    dy = Seg2([1, 2])
+    dy["geometry"] = Vector[N[1], N[2]]
+    dy["displacement 2"] = 0.0
+  
+    problem2 = DirichletProblem("displacement", 2)
+    push!(problem2, dx)
+    push!(problem2, dy)
+
+    solver = DirectSolver()
+    push!(solver, problem)
+    push!(solver, problem2)
+
+    # launch solver
+    norm = solver(0.0)
+
+    disp = e1("displacement", [1.0, 1.0], 0.0)
+    info("displacement at tip: $disp")
+    @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
+
+end
+
+# test_solver_multiple_dirichlet_bc()
+
 end