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

@@ -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

src/FEMBase.jl

@@ -26,7 +26,10 @@ include("elements.jl")
 export element_info!
 include("elements_lagrange.jl")
 include("integrate.jl")
+
 include("problems.jl")
+export set_gdofs!, get_gdofs
+
 include("assembly.jl")
 export assemble_prehook!, assemble_posthook!
 
@@ -40,7 +43,7 @@ export AbstractAnalysis, Analysis, add_problems!,
 export FieldProblem, BoundaryProblem, Problem, Element, Assembly
 export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
        Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27
-export add_elements!, add!, get_gdofs, group_by_element_type,
+export add_elements!, add!, group_by_element_type,
        get_unknown_field_name, get_unknown_field_dimension,
        get_integration_points, initialize!, assemble!
 export SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo,

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 i in conn for j=1:dim]
+    return gdofs
 end

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
@@ -259,3 +257,13 @@ end
     assemble!(problem.assembly, problem, problem.elements, 0.0)
     assemble_elements!(problem, problem.assembly, elements, 0.0)
 end
+
+@testset "test set and 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)
+end