From d466bd6275be9aa48b6c6e0d51799bcdfb6ba735 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Tue, 30 Jan 2018 07:49:13 +0200
Subject: [PATCH] Modifications to dof mapping routines

- refactor get_gdofs, removed potentially unnecessary functions
- implement set_gdofs to set dofs for an element manually
- in documentation, clarify how dofs are automatically calculated
---
 docs/src/problems.md  | 13 ++++++++++++
 src/problems.jl       | 48 +++++++++++++++++++++++--------------------
 test/test_problems.jl |  4 +---
 3 files changed, 40 insertions(+), 25 deletions(-)

diff --git a/docs/src/problems.md b/docs/src/problems.md
index f23b2f5..5470906 100644
--- a/docs/src/problems.md
+++ b/docs/src/problems.md
@@ -121,6 +121,18 @@ matrix can be parallelized (at least almost) perfectly and is not considered
 as a bottleneck when models get big enough. It's anyway a good idea to pay
 attention to the memory allocations.
 
+### Setting and getting global degrees of freedom for element
+
+[`get_gdofs`](@ref) is returning the global degrees of freedom for element. They
+can be set manually using [`set_gdofs(problem, element, dofs)`](@ref). Otherwise
+they are calculated automatically based on the problem dimension using formula
+`g(i,j,d) = d*(i-1)+j`, where `i` is local node number, `j` is the number of
+degree of freedom and `d` is the maximum number of degrees of freedom for node.
+With this formula dofs are ordered so that first comes all dofs for node 1,
+then for node 2 and so on. For 3 dofs/node we get ``(u_{11}, u_{12}, u_{13},
+u_{21}, u_{22}, u_{23}, \ldots, u_{N1}, u_{N2}, u_{N3})``, where the first index
+is node number and second index is dof number.
+
 Let's create some test problem and test our implementation:
 
 ```@example FEMBase
@@ -131,6 +143,7 @@ update!(el1, "thermal conductivity", 6.0)
 elements = [el1]
 assembly = Assembly()
 problem = Problem(Heat, "test problem", 1)
+nothing # hide
 ```
 
 Now the struct `Heat` is empty. If we need to store some problem-wide settings
diff --git a/src/problems.jl b/src/problems.jl
index cc337c1..490f5ff 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -418,32 +418,36 @@ function push!(problem::Problem, elements_::Vector...)
     end
 end
 
-function get_gdofs(element::Element, dim::Int)
-    conn = get_connectivity(element)
-    if length(conn) == 0
-        error("element connectivity not defined, cannot determine global dofs for element: $element")
-    end
-    gdofs = vec([dim*(i-1)+j for j=1:dim, i in conn])
-    return gdofs
-end
+"""
+    set_gdofs(problem, element)
 
-function empty!(problem::Problem)
-    empty!(problem.assembly)
+Set element global degrees of freedom.
+"""
+function set_gdofs(problem, element, dofs)
+    problem.dofmap[element] = dofs
 end
 
-""" Return global degrees of freedom for element.
+"""
+    get_gdofs(problem, element)
 
-Notes
------
-First look dofs from problem.dofmap, it not found, update dofmap from
-element.element connectivity using formula gdofs = [dim*(nid-1)+j for j=1:dim]
-1. look element dofs from problem.dofmap
-2. if not found, use element.connectivity to update dofmap and 1.
+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)
-        dim = get_unknown_field_dimension(problem)
-        problem.dofmap[element] = get_gdofs(element, dim)
+function get_gdofs(problem, element)
+    if haskey(problem.dofmap, element)
+        return problem.dofmap[element]
     end
-    return problem.dofmap[element]
+    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 j=1:dim for i in conn]
+    return gdofs
 end
diff --git a/test/test_problems.jl b/test/test_problems.jl
index 7700ff4..cfd28ce 100644
--- a/test/test_problems.jl
+++ b/test/test_problems.jl
@@ -98,9 +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(el, 1) == [1, 2]
     @test get_gdofs(p1, el) == [1, 2]
-    empty!(p2)
 
     p3 = Problem(P2, "P3", 1, "p")
     as = get_assembly(p3)
@@ -117,7 +115,7 @@ end
     @test haskey(p3, "f")
     @test isapprox(p3["f"].data, 1.0)
     f = p3("f", 0.0)
-    @test_throws Exception get_gdofs(e3, 1)
+    @test_throws Exception get_gdofs(problem, e3)
 end
 
 @testset "test initializing new problems" begin