[WIP] Implement new type for more flexible mapping of degrees of freedom

src/FEMBase.jl

@@ -27,8 +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("assembly.jl")
 export assemble_prehook!, assemble_posthook!, assemble_elements!

src/dofmap.jl

@@ -0,0 +1,71 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+abstract type AbstractDOFMap end
+
+"""
+    SimpleDOFMap
+
+DOF maps makes conversions between node id/dof pair and global dofs. The basic
+question to anwwer is that given some node `j` and its local dof index `i`, what
+is the corresponding row number in global matrix assembly.
+
+In the simplest case, each node has same number of degrees of freedom and nodes
+starts from 1. For example, given continuum mechanics model having ``d=3``
+displacement degrees of freedom per dof, we have ``(u_{11}, u_{12}, u_{13},
+u_{21}, u_{22}, u_{23}, \ldots, u_{N1}, u_{N2}, u_{N3})``, where, in general,
+``u_{ij}`` maps to row d*(i-1)+j in matrix assembly level.
+"""
+type SimpleDOFMap <: AbstractDOFMap
+    local_dof_indices :: Dict{Symbol, Int64}
+end
+
+function SimpleDOFMap()
+    return SimpleDOFMap(Dict())
+end
+
+"""
+    set_local_dof_indices!(dofmap, local_dof_indices)
+
+Set local dof indices for dofmap, connecting together dof abbreviation and its
+order in node.
+
+# Example
+
+```julia
+local_dof_indices = Dict(1 => :u1, 2 => :u2, 3 => :ur3, 4 => :T)
+set_local_dof_indices!(dofmap, local_dof_indices)
+
+# output
+
+nothing
+```
+"""
+function set_local_dof_indices!(dofmap, local_dof_indices)
+    dofmap.local_dof_indices = local_dof_indices
+    return nothing
+end
+
+"""
+    get_gdofs(dofmap, nodes, dof_names)
+
+Return global dofs for some nodes and dofs. With `SimpleDOFMap`, this is a
+trivial calculation operation `gdof(i,j) = d(j-1)+i`, where ``j`` is node id,
+``i`` is the local dof index and ``d`` is number of dofs per node.
+"""
+function get_gdofs(dofmap::SimpleDOFMap, nodes, dof_names)
+    ldi = dofmap.local_dof_indices
+    max_dim = length(ldi)
+    return (max_dim*(j-1)+ldi[n] for j in nodes for n in dof_names)
+end
+
+function get_gdofs!(gdofs, dofmap::SimpleDOFMap, nodes, dof_names)
+    for (k,j) in enumerate(get_gdofs(dofmap, nodes, dof_names))
+        gdofs[k] = j
+    end
+    return gdofs
+end
+
+function DOFMap()
+    return SimpleDOFMap()
+end

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,45 @@ 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.
+
+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 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
+        FEMBase.get_dof_names(::Problem{$P}, ::Element) = $dof_names
+    end
+    eval(code)
+    return dof_names
+end
+
+"""
+    get_dofmap(problem)
 
-Set element global degrees of freedom.
+Return `DOFMap` for problem.
 """
-function set_gdofs!(problem, element, dofs)
-    problem.dofmap[element] = dofs
+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/runtests.jl

@@ -16,6 +16,7 @@ push!(test_files, "test_solvers.jl")
 push!(test_files, "test_analysis.jl")
 push!(test_files, "test_test.jl")
 push!(test_files, "test_types.jl")
+push!(test_files, "test_dofmap.jl")
 
 @testset "FEMBase.jl" begin
     for fn in test_files

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_dofmap.jl

@@ -0,0 +1,34 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+# # Mapping nodal degrees of freedom to matrix rows (global dofs)
+
+using FEMBase
+using FEMBase.Test
+using FEMBase: DOFMap, set_local_dof_indices!
+
+# With `SimpleDOFMap`, matrix row indices, what we call to global dofs, are
+# calculated based on the number of node ``j`` and maximum number of dofs per
+# node ``d``. So for example, if each node has two dofs ``u_1``, ``u_2``, we
+# have ``(u_{11}, u_{12}, u_{21}, u_{22}, \ldots, u_{N1}, u_{N2})``. In general,
+# this formula is then `gdof(i,j) = d*(1-j)+i`. `SimpleDOFMap` mode does not
+# need any storage of mapping because matrix row is explicitly determined based
+# on node id number and local dof index.
+
+dm = DOFMap()
+set_local_dof_indices!(dm, Dict(:u1=>1, :u2=>2, :T=>3))
+
+# Accessing of data is done using `get_gdofs` or in-place version `get_gdofs!`:
+
+@test collect(get_gdofs(dm, (1, 3), (:u1, :u2))) == [1, 2, 7, 8]
+
+# In general, `get_gdofs` takes two iterables, one describing nodes and another
+# describing name of the dofs. Order does matter:
+
+@test collect(get_gdofs(dm, (3, 1), (:u1, :u2))) == [7, 8, 1, 2]
+@test collect(get_gdofs(dm, (1, 3), (:u2, :u1))) == [2, 1, 8, 7]
+
+# There is also in-place version to get dofs:
+
+gdofs = zeros(Int64, 4)
+@test collect(get_gdofs!(dm, (1, 2), (:u1, :u2))) == [1, 2, 3, 4]

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]