diff --git a/src/FEMBase.jl b/src/FEMBase.jl index c3f4492..adb46fd 100644 --- a/src/FEMBase.jl +++ b/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! diff --git a/src/dofmap.jl b/src/dofmap.jl new file mode 100644 index 0000000..1e333b1 --- /dev/null +++ b/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 diff --git a/src/problems.jl b/src/problems.jl index 25128c1..98793c9 100644 --- a/src/problems.jl +++ b/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,12 +431,34 @@ 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 """ @@ -432,22 +466,10 @@ end 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 diff --git a/src/sparse.jl b/src/sparse.jl index e6002ff..463584c 100644 --- a/src/sparse.jl +++ b/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. """ diff --git a/src/test.jl b/src/test.jl index d6b25a5..0e2267e 100644 --- a/src/test.jl +++ b/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` diff --git a/test/runtests.jl b/test/runtests.jl index c8162ff..08e7b10 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_assembly.jl b/test/test_assembly.jl index b76fed5..75b1da9 100644 --- a/test/test_assembly.jl +++ b/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) diff --git a/test/test_dofmap.jl b/test/test_dofmap.jl new file mode 100644 index 0000000..d75e7c7 --- /dev/null +++ b/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] diff --git a/test/test_problems.jl b/test/test_problems.jl index 62a1e74..1587e3a 100644 --- a/test/test_problems.jl +++ b/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 diff --git a/test/test_sparse.jl b/test/test_sparse.jl index 7270cc9..95de4c6 100644 --- a/test/test_sparse.jl +++ b/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]