Merge 6551f2a into 110acde

src/FEMBase.jl

@@ -36,5 +36,6 @@ export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
 export update!, add_elements!, add!, get_gdofs, group_by_element_type,
        get_unknown_field_name, get_unknown_field_dimension,
        get_integration_points, initialize!, assemble!
+export DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV
 
 end

src/assembly.jl

@@ -44,6 +44,11 @@ function assemble!(problem::Problem, time=0.0; auto_initialize=true)
     return true
 end
 
+function assemble!{P}(assembly::Assembly, problem::Problem{P}, element::Element, time)
+    warn("One must define assemble! function for problem of type $P. Not doing anything.")
+    return nothing
+end
+
 function assemble!(assembly::Assembly, problem::Problem, elements::Vector{Element}, time)
     warn("assemble!() this is default assemble operation, decreased performance can be expected without preallocation of memory!")
     for element in elements

src/problems.jl

@@ -395,6 +395,11 @@ function get_unknown_field_dimension(problem::Problem)
     return problem.dimension
 end
 
+function get_unknown_field_name{P<:AbstractProblem}(::Type{P})
+    warn("The name of unknown field (e.g. displacement, temperature, ...) of the problem type must be given by defining function `get_unknown_field_name`")
+    return "N/A"
+end
+
 """ Return the name of the unknown field of this problem. """
 function get_unknown_field_name{P}(problem::Problem{P})
     return get_unknown_field_name(P)

src/types.jl

@@ -72,5 +72,5 @@ end
 
 function IP(id, weight, coords::Vector)
     warn("Consider giving coords as Tuple.")
-    return IP(id, weight, [c for c in coords], Dict(), IntegrationPoint())
+    return IP(id, weight, tuple(coords...), Dict(), IntegrationPoint())
 end

test/runtests.jl

@@ -6,6 +6,13 @@ using TimerOutputs
 const to = TimerOutput()
 
 test_files = String[]
+push!(test_files, "test_add_elements.jl")
+push!(test_files, "test_common_failures.jl")
+push!(test_files, "test_elements.jl")
+push!(test_files, "test_fields.jl")
+push!(test_files, "test_integrate.jl")
+push!(test_files, "test_sparse.jl")
+push!(test_files, "test_types.jl")
 
 @testset "FEMBase.jl" begin
     for fn in test_files

test/test_add_elements.jl

@@ -0,0 +1,16 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using Base.Test
+
+type Dummy <: FEMBase.FieldProblem
+end
+
+@testset "add elements to problem" begin
+    problem = Problem(Dummy, "test", 2)
+    element = Element(Quad4, [1, 2, 3, 4])
+    elements = [element]
+    add_elements!(problem, elements)
+    @test problem.elements[1] == element
+end

test/test_common_failures.jl

@@ -0,0 +1,15 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using Base.Test
+
+type Dummy <: FieldProblem end
+
+@testset "no unknown field name or assemble!-function defined" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    pr = Problem(Dummy, "problem", 2)
+    add_elements!(pr, [el])
+    assemble!(pr)
+    @test true
+end

test/test_elements.jl

@@ -0,0 +1,226 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using Base.Test
+
+using FEMBase: Field
+using FEMBase: group_by_element_type
+using FEMBase: get_local_coordinates, inside
+using FEMBase: get_basis, get_dbasis, get_integration_order
+using FEMBase: get_reference_element_coordinates, get_reference_coordinates
+using FEMBase: get_element_type, is_element_type, get_element_id, filter_by_element_type
+
+@testset "add time dependent field to element" begin
+    el = Element(Seg2, [1, 2])
+    u1 = Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]
+    u2 = Vector{Float64}[[1.0, 1.0], [1.0, 1.0]]
+    update!(el, "displacement", 0.0 => u1)
+    update!(el, "displacement", 1.0 => u2)
+    @test length(el["displacement"]) == 2
+    @test isapprox(el("displacement", [0.0], 0.0), [0.0, 0.0])
+    @test isapprox(el("displacement", [0.0], 0.5), [0.5, 0.5])
+    @test isapprox(el("displacement", [0.0], 1.0), [1.0, 1.0])
+    el2 = Element(Poi1, [1])
+    update!(el2, "force 1", 0.0 => 1.0)
+end
+
+@testset "add CVTV field to element" begin
+    el = Element(Seg2, [1, 2])
+    f(xi, time) = xi[1]*time
+    update!(el, "my field", f)
+    v = el("my field", [1.0], 2.0)
+    @test isapprox(v, 2.0)
+end
+
+@testset "add DCTI to element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    update!(el, "displacement load", DCTI([4.0, 8.0]))
+    @test isa(el["displacement load"], DCTI)
+    @test !isa(el["displacement load"].data, DCTI)
+    update!(el, "displacement load 2", [4.0, 8.0])
+    @test isa(el["displacement load 2"], DCTI)
+    update!(el, "temperature", [1.0, 2.0, 3.0, 4.0])
+    @test isa(el["temperature"], DVTI)
+    @test isapprox(el("displacement load", [0.0, 0.0], 0.0), [4.0, 8.0])
+end
+
+@testset "interpolate DCTI from element" begin
+    el = Element(Seg2, [1, 2])
+    update!(el, "foobar", 1.0)
+    fb = el("foobar", [0.0], 0.0)
+    @test isa(fb, Float64)
+    @test isapprox(fb, 1.0)
+end
+
+@testset "add elements to elements" begin
+    el1 = Element(Seg2, [1, 2])
+    el2 = Element(Seg2, [3, 4])
+    update!(el1, "master elements", [el2])
+    lst = el1("master elements", 0.0)
+    @test isa(lst, Vector)
+end
+
+@testset "extend basis" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    expected = [
+        0.25 0.00 0.25 0.00 0.25 0.00 0.25 0.00
+        0.00 0.25 0.00 0.25 0.00 0.25 0.00 0.25]
+    @test isapprox(el([0.0, 0.0], 0.0, 2), expected)
+end
+
+@testset "group elements" begin
+    e1 = Element(Seg2, [1, 2])
+    e2 = Element(Quad4, [1, 2, 3, 4])
+    elements = [e1, e2]
+    r = group_by_element_type(elements)
+    @test length(r) == 2
+    @test first(r[Element{Seg2}]) == e1
+    @test first(r[Element{Quad4}]) == e2
+    @test get_element_type(e1) == Seg2
+    e1.id = 1
+    @test get_element_id(e1) == 1
+    @test is_element_type(e1, Seg2)
+    seg2_elements = filter_by_element_type(elements, Seg2)
+end
+
+@testset "inverse isoparametric mapping" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    X1 = el("geometry", [0.1, 0.2], time)
+    xi = get_local_coordinates(el, X1, time)
+    X2 = el("geometry", xi, time)
+    info("X1 = $X1, X2 = $X2")
+    @test isapprox(X1, X2)
+end
+
+@testset "inside of linear element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    @test inside(el, [0.5, 0.5], time) == true
+    @test inside(el, [1.0, 0.5], time) == true
+    @test inside(el, [1.0, 1.0], time) == true
+    @test inside(el, [1.01, 1.0], time) == false
+    @test inside(el, [1.0, 1.01], time) == false
+end
+
+@testset "inside of quadratic element" begin
+    el = Element(Tri6, [1, 2, 3, 4, 5, 6])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [0.0, 1.0],
+        4 => [0.5, 0.2],
+        5 => [0.8, 0.6],
+        6 => [-0.2, 0.5])
+    update!(el, "geometry", X)
+    p = [0.94, 0.3] # visually checked to be inside
+    @test inside(el, p, 0.0) == true
+    p = [-0.2, 0.8] # visually checked to be outside
+    @test inside(el, p, 0.0) == false
+end
+
+@testset "dict field" begin
+    el = Element(Seg2, [1, 2])
+    X = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [0.5, 0.5])
+    f = Field(X)
+    debug("field = $f")
+    #update!(el, "geometry", X)
+    el["geometry"] = f
+    @test isapprox(el("geometry")[1], [0.0, 0.0])
+    @test isapprox(el("geometry", 0.0)[1], [0.0, 0.0])
+    @test isapprox(el("geometry", 0.0)[3], [0.5, 0.5])
+    @test isapprox(el("geometry", [0.0], 0.0), [0.5, 0.0])
+end
+
+@testset "test nodal element" begin
+    el = Element(Poi1, [1])
+    @test get_basis(el, [0.0], 1.0) == [1]
+    @test get_dbasis(el, [0.0], 1.0) == [0]
+    @test isapprox(el([0.0], 0.0, Val{:detJ}), 1.0)
+    @test get_integration_order(el.properties) == 1
+    ips = get_integration_points(el.properties, 1)
+    @test length(ips) == 1
+    @test length(ips[1]) == 2
+    @test size(el) == (0, 1)
+    @test length(el) == 1
+    @test get_reference_element_coordinates(Poi1) == Vector{Float64}[[0.0]]
+end
+
+@testset "get reference coordinates" begin
+    el = Element(Seg2, [1, 2])
+    X = get_reference_coordinates(el)
+    @test isapprox(X[1][1], -1.0)
+    @test isapprox(X[2][1],  1.0)
+end
+
+@testset "function field" begin
+    function f(element, xi, time)
+        return xi[1]*time
+    end
+    e1 = Element(Seg2, [1, 2])
+    e1["f"] = f
+    @test isapprox(e1("f", [1.0], 1.0), 1.0)
+end
+
+@testset "index" begin
+    e1 = Element(Seg2, [1, 2])
+    update!(e1, "f", 0.0 => 1.0)
+    update!(e1, "f", 1.0 => 2.0)
+    @test isapprox(last(e1, "f").data, 2.0)
+end
+
+@testset "gradient, jacobian and determinant of jacobian" begin
+    X = Dict(
+             1 => [0.0, 0.0],
+             2 => [2.0, 0.0],
+             3 => [2.0, 2.0],
+             4 => [0.0, 2.0])
+    element = Element(Quad4, [1, 2, 3, 4])
+    element.fields["geometry"] = Field(X)
+    J = element([0.0, 0.0], 0.0, Val{:Jacobian})
+    detJ = element([0.0, 0.0], 0.0, Val{:detJ})
+    @test isapprox(det(J), detJ)
+    @test isapprox(detJ, 1.0)
+    element2 = Element(Seg2, [1, 2])
+    update!(element2, "geometry", X)
+    detJ = element2([0.0], 0.0, Val{:detJ})
+    @test isapprox(detJ, 1.0)
+    X = Dict(
+             1 => [0.0, 0.0, 0.0],
+             2 => [2.0, 0.0, 0.0],
+             3 => [2.0, 2.0, 0.0],
+             4 => [0.0, 2.0, 0.0])
+    element3 = Element(Quad4, [1, 2, 3, 4])
+    update!(element3, "geometry", X)
+    detJ = element3([0.0, 0.0], 0.0, Val{:detJ})
+    @test isapprox(detJ, 1.0)
+    grad = element([0.0, 0.0], 0.0, Val{:Grad})
+    @test isapprox(grad, 1/4*[-1.0 1.0 1.0 -1.0; -1.0 -1.0 1.0 1.0])
+    update!(element3, "youngs modulus", 1.0)
+    @test isapprox(element3("youngs modulus", [0.0], 0.0), 1.0)
+    update!(element3, "test", [1.0, 2.0, 3.0])
+    test_val = element3("test", [0.0], 1.0)
+end
+
+@testset "update data to element" begin
+    element = Element(Seg2, [1, 2])
+    X = Dict(1 => 1.0)
+    @test_throws KeyError update!(element, "fail", X)
+    @test_throws KeyError update!(element, "fail2", 0.0 => X)
+    update!(element, "success", 0.0 => [1.0, 2.0])
+    update!(element, "success2", 0.0 => [1, 2])
+    update!(element, "success3", 0.0 => Vector{Float64}[[1.0], [2.0]])
+end

test/test_elements_2.jl

@@ -0,0 +1,54 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using FEMBase: get_local_coordinates, inside
+using Base.Test
+
+@testset "inverse isoparametric mapping" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    X1 = el("geometry", [0.1, 0.2], time)
+    xi = get_local_coordinates(el, X1, time)
+    X2 = el("geometry", xi, time)
+    info("X1 = $X1, X2 = $X2")
+    @test isapprox(X1, X2)
+end
+
+@testset "inside of linear element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    @test inside(el, [0.5, 0.5], time) == true
+    @test inside(el, [1.0, 0.5], time) == true
+    @test inside(el, [1.0, 1.0], time) == true
+    @test inside(el, [1.01, 1.0], time) == false
+    @test inside(el, [1.0, 1.01], time) == false
+end
+
+@testset "inside of quadratic element" begin
+    el = Element(Tri6, [1, 2, 3, 4, 5, 6])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [0.0, 1.0],
+        4 => [0.5, 0.2],
+        5 => [0.8, 0.6],
+        6 => [-0.2, 0.5])
+    update!(el, "geometry", X)
+    p = [0.94, 0.3] # visually checked to be inside
+    @test inside(el, p, 0.0) == true
+    p = [-0.2, 0.8] # visually checked to be outside
+    @test inside(el, p, 0.0) == false
+end