Use DOFMap in Assembly

Also, losen data type definitions in `add!` so that generators works 
also when adding data do sparse matricse.

Author: ahojukka5

src/FEMBase.jl

@@ -27,10 +27,12 @@ export element_info!
 include("elements_lagrange.jl")
 include("integrate.jl")
 
+include("dofmap.jl")
+export get_gdofs
+
 include("problems.jl")
-export set_gdofs!, get_gdofs, get_formulation_type
+export get_formulation_type
 
-include("dofmap.jl")
 
 include("assembly.jl")
 export assemble_prehook!, assemble_posthook!, assemble_elements!

src/problems.jl

@@ -100,7 +100,7 @@ type Problem{P<:AbstractProblem}
     dimension :: Int                 # degrees of freedom per node
     parent_field_name :: AbstractString # (optional) name of the parent field e.g. "displacement"
     elements :: Vector{Element}
-    dofmap :: Dict{Element, Vector{Int64}} # connects the element local dofs to the global dofs
+    dofmap :: AbstractDOFMap         # connects the element local dofs to the global dofs
     assembly :: Assembly
     fields :: Dict{String, AbstractField}
     postprocess_fields :: Vector{String}
@@ -125,7 +125,13 @@ prob1 = Problem(Elasticity, "this is my problem", 3)
 
 """
 function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
-    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P())
+    dofmap = DOFMap()
+    assembly = Assembly()
+    local_dof_indices = Dict(Symbol("u$i") => i for i=1:dimension)
+    set_local_dof_indices!(dofmap, local_dof_indices)
+    problem = Problem{P}(name, dimension, "none", [], dofmap, assembly,
+                         Dict(), Vector(), P())
+    return problem
 end
 
 """
@@ -139,14 +145,20 @@ julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
 solver.
 """
 function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name)
-    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), P())
+    dofmap = DOFMap()
+    assembly = Assembly()
+    local_dof_indices = Dict(Symbol("u$i") => i for i=1:dimension)
+    set_local_dof_indices!(dofmap, local_dof_indices)
+    problem = Problem{P}(name, dimension, parent_field_name, [], dofmap,
+                         assembly, Dict(), Vector(), P())
+    return problem
 end
 
 function get_formulation_type(::Problem)
     return :incremental
 end
 
-""" 
+"""
     get_unknown_field_dimension(problem)
 
 Return the dimension of the unknown field of this problem.
@@ -419,35 +431,47 @@ function push!(problem::Problem, elements_::Vector...)
 end
 
 """
-    set_gdofs!(problem, element)
+    get_dof_names(problem, element)
+
+Return an iterable containing the abbreviations of dof names for this element
+and problem. This could be e.g. `(:u1, :u2, :T)` for two-dimensional problem
+having two displacement dofs and temperature dof.
 
-Set element global degrees of freedom.
+Each new Problem should define this function. If now defined, assumption is that
+dofs for problem are (:u1, :u2, ..., :uN) where N is the unknown field dimension.
 """
-function set_gdofs!(problem, element, dofs)
-    problem.dofmap[element] = dofs
+function get_dof_names(problem::Problem{P}, ::Element) where {P}
+    dim = get_unknown_field_dimension(problem)
+    dof_names = (Symbol("u$i") for i=1:dim)
+    warn("Define function FEMBase.get_dof_names(::Problem{$P}, ::Element).")
+    warn("Assuming that dofs for problem $P are $(collect(dof_names)).")
+    code = quote
+        function FEMBase.get_dof_names(::Problem{$P}, ::Element)
+            return $dof_names
+        end
+    end
+    eval(code)
+    return dof_names
+end
+
+"""
+    get_dofmap(problem)
+
+Return `DOFMap` for problem.
+"""
+function get_dofmap(problem)
+    return problem.dofmap
 end
 
 """
     get_gdofs(problem, element)
 
 Return the global degrees of freedom for element.
 
-First make lookup from problem dofmap. If not defined there, make implicit
-assumption that dofs follow formula `gdofs = [dim*(nid-1)+j for j=1:dim]`,
-where `nid` is node id and `dim` is the dimension of problem. This formula
-arranges dofs so that first comes all dofs of node 1, then node 2 and so on:
-(u11, u12, u13, u21, u22, u23, ..., un1, un2, un3) for 3 dofs/node setting.
 """
-function get_gdofs(problem::Problem, element::Element)
-    if haskey(problem.dofmap, element)
-        return problem.dofmap[element]
-    end
-    conn = get_connectivity(element)
-    if length(conn) == 0
-        error("element connectivity not defined, cannot determine global ",
-              "degrees of freedom for element #: $(element.id)")
-    end
-    dim = get_unknown_field_dimension(problem)
-    gdofs = [dim*(i-1)+j for i in conn for j=1:dim]
-    return gdofs
+function get_gdofs(problem::Problem{P}, element::Element) where {P}
+    nodes = get_connectivity(element)
+    dof_names = get_dof_names(problem, element)
+    dofmap = get_dofmap(problem)
+    return get_gdofs(dofmap, nodes, dof_names)
 end

src/sparse.jl

@@ -116,11 +116,10 @@ Example
  0.0  0.0  0.0  0.0  0.0  8.0  9.0  10.0
 
 """
-function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix)
-    n, m = size(data)
-    for j=1:m
-        for i=1:n
-            add!(A, dofs1[i], dofs2[j], data[i,j])
+function add!(A::SparseMatrixCOO, dofs1, dofs2, data)
+    for (j, d1) in enumerate(dofs1)
+        for (i, d2) in enumerate(dofs2)
+            add!(A, d1, d2, data[i,j])
         end
     end
     return nothing
@@ -134,15 +133,11 @@ function add!(A::SparseMatrixCOO, B::SparseMatrixCSC)
 end
 
 """ Add new data to COO Sparse vector. """
-function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1)
-    if length(dofs) != length(data)
-        info("dofs = $dofs")
-        info("data = $(vec(data))")
-        error("when adding to sparse vector dimension mismatch!")
+function add!(A::SparseMatrixCOO, dofs, data::Array{Float64}, dim::Int=1)
+    for (d, v) in zip(dofs, data)
+        add!(A, d, dim, v)
     end
-    append!(A.I, dofs)
-    append!(A.J, dim*ones(Int, length(dofs)))
-    append!(A.V, vec(data))
+    return nothing
 end
 
 """ Add SparseVector to SparseVectorCOO. """

src/test.jl

@@ -112,7 +112,7 @@ function assemble_elements!{E}(problem::Problem{Dirichlet},
     for element in elements
         for i=1:dim
             haskey(element, "$name $i") || continue
-            gdofs = get_gdofs(problem, element)
+            gdofs = collect(get_gdofs(problem, element))
             ldofs = gdofs[i:dim:end]
             xis = get_reference_element_coordinates(E)
             for (ldof, xi) in zip(ldofs, xis)
@@ -263,7 +263,7 @@ end
     test_resource(resource::String) -> String
 
 `@test_resource(resource)` expands to a string containing full path to the some
-`resource` what is needed by a test. 
+`resource` what is needed by a test.
 
 # Example
 
@@ -274,7 +274,7 @@ One needs a test mesh file `mesh.inp`. Then, function
 mesh_file = @test_resource("mesh.inp")
 ```
 
-expands to 
+expands to
 
     `~/.julia/v0.6/Models/test/test_run_model/mesh.inp`
 

test/test_assembly.jl

@@ -59,7 +59,7 @@ end
     M = full(p.assembly.M)
 
     M_expected = [
-        6   1   1   1  -4  -6  -4  -4  -6  -6 
+        6   1   1   1  -4  -6  -4  -4  -6  -6
         1   6   1   1  -4  -4  -6  -6  -4  -6
         1   1   6   1  -6  -4  -4  -6  -6  -4
         1   1   1   6  -6  -6  -6  -4  -4  -4
@@ -83,8 +83,8 @@ end
     A = Assembly()
     B = Assembly()
     @test isempty(A)
-    add!(A.K, [1,2,3,4], [1,2,3,4], [1.0 2.0; 3.0 4.0])
-    add!(B.K, [1,2,3,4], [1,2,3,4], [1.0 2.0; 3.0 4.0])
+    add!(A.K, [1,2], [1,2], [1.0 2.0; 3.0 4.0])
+    add!(B.K, [1,2], [1,2], [1.0 2.0; 3.0 4.0])
     @test isapprox(A, B)
     @test !isempty(A)
     empty!(A)

test/test_problems.jl

@@ -98,7 +98,7 @@ end
     println("elgeom is ", el("geometry", 0.0))
     @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0])
     @test get_parent_field_name(p2) == "p"
-    @test get_gdofs(p1, el) == [1, 2]
+    @test collect(get_gdofs(p1, el)) == [1, 2]
 
     p3 = Problem(P2, "P3", 1, "p")
     as = get_assembly(p3)
@@ -203,7 +203,7 @@ function assemble_elements!{E}(problem::Problem{DirBC},
     for element in elements
         for i=1:dim
             haskey(element, "$name $dim") || continue
-            gdofs = get_gdofs(problem, element)
+            gdofs = collect(get_gdofs(problem, element))
             ldofs = gdofs[i:dim:end]
             xis = get_reference_element_coordinates(E)
             for (ldof, xi) in zip(ldofs, xis)
@@ -258,12 +258,8 @@ end
     assemble_elements!(problem, problem.assembly, elements, 0.0)
 end
 
-@testset "test set and get global dofs for element" begin
+@testset "test get global dofs for element" begin
     element = Element(Seg2, [1, 2])
     problem = Problem(P3, "P3", 2)
-    @test get_gdofs(problem, element) == [1, 2, 3, 4]
-    set_gdofs!(problem, element, [2, 3, 4, 5])
-    @test get_gdofs(problem, element) == [2, 3, 4, 5]
-    element2 = Element(Seg2, Int[])
-    @test_throws ErrorException get_gdofs(problem, element2)
+    @test collect(get_gdofs(problem, element)) == [1, 2, 3, 4]
 end

test/test_sparse.jl

@@ -21,11 +21,6 @@ end
     @test isapprox(full(b), full(b2))
 end
 
-@testset "Failure to add data to sparse vector due dimension mismatch" begin
-    b = SparseVectorCOO()
-    @test_throws ErrorException add!(b, [1, 2], [1.0, 2.0, 3.0])
-end
-
 @testset "Test combining of SparseMatrixCOO" begin
     k = convert(Matrix{Float64}, reshape(collect(1:9), 3, 3))
     dofs1 = [1, 2, 3]