Minor fixes (#9)

* export get_reference_element_coordinates
* refactor assembly.jl
* refactor fields.jl, implement == to compare fields
* use add! instead of push! to add single entry for sparse matrix
* test comparison of fields, simplified creating of DVTI
* assemble.jl: test creating new problems
* fixed deprecated syntax, remains backward compatible to JuliaFEM.

src/FEMBase.jl

@@ -22,7 +22,7 @@ end
 export info, debug
 
 using FEMBasis
-using FEMBasis: AbstractBasis, jacobian
+using FEMBasis: AbstractBasis, jacobian, get_reference_element_coordinates
 import FEMBasis: interpolate
 using FEMQuad: get_quadrature_points
 include("fields.jl")
@@ -48,6 +48,7 @@ export SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo,
 export is_field_problem, is_boundary_problem, get_elements,
        get_connectivity,
        get_parent_field_name, get_reference_coordinates,
+       get_reference_element_coordinates,
        get_assembly, get_nonzero_rows, get_nonzero_columns,
        eval_basis!, get_basis, get_dbasis, grad!,
        assemble_mass_matrix!, get_local_coordinates, inside,

src/assembly.jl

@@ -11,73 +11,98 @@ function isapprox(a1::Assembly, a2::Assembly)
     return T
 end
 
-function assemble_prehook!(::Problem, time) end
+function assemble_prehook!{T<:Number}(::Problem, ::T) end
 
-function assemble_posthook!(::Problem, time) end
+function assemble_posthook!{T<:Number}(::Problem, ::T) end
 
-function assemble!(problem::Problem, time=0.0; auto_initialize=true)
-    if !isempty(problem.assembly)
-        warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
-    end
-    if isempty(problem.elements)
-        warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
-    else
-        first_element = first(problem.elements)
-        unknown_field_name = get_unknown_field_name(problem)
-        if !haskey(first_element, unknown_field_name)
-            warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
-            if auto_initialize
-                info("Initializing problem $(problem.name) at time $time automatically.")
-                initialize!(problem, time)
-            end
-        end
-    end
-    assembly = get_assembly(problem)
-    assemble_prehook!(problem, time)
-    for (element_type, elements) in group_by_element_type(get_elements(problem))
-        assemble!(problem, assembly, elements, time)
-    end
-    assemble_posthook!(problem, time)
-    return true
-end
-
-function assemble!{P}(::Assembly, ::Problem{P}, element::Element, time::Float64)
+# will be deprecated
+function assemble!{P}(::Assembly, ::Problem{P}, ::Element, ::Any)
     warn("One must define assemble! function for problem of type $P. ",
          "Not doing anything.")
     return nothing
 end
 
+# will be deprecated
+function assemble!{P}(assembly::Assembly, problem::Problem{P},
+                      elements::Vector{Element}, time)
+    warn("This is default assemble! function. Decreased performance can be ",
+         "expected without preallocation of memory. One should implement ",
+         "`assemble_elements!(problem, assembly, elements, time)` function.")
+    for element in elements
+        assemble!(assembly, problem, element, time)
+    end
+    return nothing
+end
+
 """
-    assemble!(problem, assembly, elements, time)
+    assemble_elements!(problem, assembly, elements, time)
+
+Assemble elements for problem.
 
-This should be overridden with own custom assemble operation.
+This should be overridden with own `assemble_elements!`-implementation.
 """
-function assemble!{E}(problem::Problem, assembly::Assembly,
-                      elements::Vector{Element{E}}, time::Float64)
+function assemble_elements!{E}(problem::Problem, assembly::Assembly,
+                      elements::Vector{Element{E}}, time)
     elements2 = convert(Vector{Element}, elements)
     assemble!(assembly, problem, elements2, time)
 end
 
-function assemble!{P}(assembly::Assembly, problem::Problem{P},
-                      elements::Vector{Element}, time)
-    warn("This is default assemble! function. Decreased performance can be ",
-         "expected without preallocation of memory.")
-    for element in elements
-        assemble!(assembly, problem, element, time)
+function assemble!(problem::Problem, time)
+
+    assemble_prehook!(problem, time)
+    elements = get_elements(problem)
+    assembly = get_assembly(problem)
+
+    if !isempty(assembly)
+        warn("Assembling elements for problem $(problem.name): problem.assembly ",
+             "is not empty before assembling. This is probably causing unexpected ",
+             "results. To remove old assembly, use `empty!(problem.assembly)`")
+        assemble_posthook!(problem, time)
+        return nothing
     end
+
+    if isempty(elements)
+        warn("Assembling elements for problem $(problem.name): problem.elements ",
+             "is empty, there is no elements in problem. Before assembling ",
+             "problem, elements must be added using ",
+             "`add_elements!(problem, elements)`.")
+        assemble_posthook!(problem, time)
+        return nothing
+    end
+
+    first_element = first(elements)
+    unknown_field_name = get_unknown_field_name(problem)
+    if !haskey(first_element, unknown_field_name)
+        warn("Assembling elements for problem $(problem.name): seems that ",
+             "problem is uninitialized. To initialize problem, use ",
+             "`initialize!(problem, time)`.")
+        info("Initializing problem $(problem.name) at time $time automatically.")
+        initialize!(problem, time)
+    end
+
+    for (element_type, elements) in group_by_element_type(elements)
+        assemble_elements!(problem, assembly, elements, time)
+    end
+    assemble_posthook!(problem, time)
     return nothing
 end
 
-function assemble_mass_matrix!(problem::Problem, time)
+function assemble!(problem::Problem)
+    warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)")
+    assemble!(problem, 0.0)
+end
+
+function assemble_mass_matrix!(problem::Problem, time::Float64)
     if !isempty(problem.assembly.M)
-        info("Mass matrix for $(problem.name) is already assembled, skipping assemble routine")
-        return
+        info("Mass matrix for $(problem.name) is already assembled, ",
+             "not assembling.")
+        return nothing
     end
     elements = get_elements(problem)
     for (element_type, elements) in group_by_element_type(get_elements(problem))
         assemble_mass_matrix!(problem::Problem, elements, time)
     end
-    return
+    return nothing
 end
 
 function assemble_mass_matrix!{Basis}(problem::Problem, elements::Vector{Element{Basis}}, time)
@@ -137,7 +162,7 @@ function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{Tet10}
         -6 -4 -6 -4 16 16  8 16 32 16
         -6 -6 -4 -4  8 16 16 16 16 32]
 
-    function is_CM(element::Element{Tet10}, X; rtol=1.0e-6)
+    function is_CM(::Element{Tet10}, X; rtol=1.0e-6)
         isapprox(X[5],  1/2*(X[1]+X[2]); rtol=rtol) || return false
         isapprox(X[6],  1/2*(X[2]+X[3]); rtol=rtol) || return false
         isapprox(X[7],  1/2*(X[3]+X[1]); rtol=rtol) || return false

src/elements.jl

@@ -149,14 +149,14 @@ function interpolate(element::Element, field_name::String, time::Float64)
     end
 end
 
-function get_basis{B}(element::Element{B}, ip, time)
+function get_basis{B}(element::Element{B}, ip, ::Any)
     T = typeof(first(ip))
     N = zeros(T, 1, length(element))
     eval_basis!(B, N, tuple(ip...))
     return N
 end
 
-function get_dbasis{B}(element::Element{B}, ip, time)
+function get_dbasis{B}(element::Element{B}, ip, ::Any)
     T = typeof(first(ip))
     dN = zeros(T, size(element)...)
     eval_dbasis!(B, dN, tuple(ip...))
@@ -221,39 +221,30 @@ end
 
 function interpolate(element::Element, field_name, ip, time)
     field = element[field_name]
-    return interpolate_(element, field, ip, time)
+    return interpolate_field(element, field, ip, time)
 end
 
-@lintpragma("Ignore unused ip")
-@lintpragma("Ignore unused element")
-function interpolate_(element::Element, field::DCTI, ip, time::Float64)
-    return field.data
+function interpolate_field{T<:Union{DCTI,DCTV}}(::Element, field::T, ::Any, time)
+    return interpolate_field(field, time)
 end
 
-function interpolate_(element::Element, field::DCTV, ip, time::Float64)
-    return interpolate_(field, time)
+function interpolate_field{T<:Union{DVTV,DVTI}}(element::Element, field::T, ip, time)
+    data = interpolate_field(field, time)
+    basis = get_basis(element, ip, time)
+    N = length(basis)
+    return sum(data[i]*basis[i] for i=1:N)
 end
 
-function interpolate_(element::Element, field::CVTV, ip, time::Float64)
-    return field(ip, time)
+function interpolate_field{T<:Union{DVTVd,DVTId}}(element::Element, field::T, ip, time)
+    data = interpolate_field(field, time)
+    basis = element(ip, time)
+    N = length(element)
+    c = get_connectivity(element)
+    return sum(data[c[i]]*basis[i] for i=1:N)
 end
 
-function interpolate_{F<:AbstractField}(element::Element, field::F, ip, time::Float64)
-    field_ = interpolate(field, time)
-    basis = element(ip, time)
-    n = length(element)
-    if isa(field_, Tuple)
-        m = length(field_)
-        if n != m
-            error("Error when trying to interpolate field $field at coords $ip ",
-                  "and time $time: element length is $n and field length is $m, ",
-                  "f = Nᵢfᵢ makes no sense!")
-        end
-        return sum(field_[i]*basis[i] for i=1:n)
-    else
-        c = get_connectivity(element)
-        return sum(field_[c[i]]*basis[i] for i=1:n)
-    end
+function interpolate_field(::Element, field::CVTV, ip, time)
+    return field(ip, time)
 end
 
 function size(element::Element, dim)

src/elements_lagrange.jl

@@ -8,24 +8,23 @@ using FEMBasis
 type Poi1 <: AbstractBasis
 end
 
-function get_basis(element::Element{Poi1}, ip, time)
+function get_basis(::Element{Poi1}, ::Any, ::Any)
     return [1]
 end
 
-function get_dbasis(element::Element{Poi1}, ip, time)
+function get_dbasis(::Element{Poi1}, ::Any, ::Any)
     return [0]
 end
 
-function (element::Element{Poi1})(ip, time::Float64, ::Type{Val{:detJ}})
+function (::Element{Poi1})(::Any, ::Float64, ::Type{Val{:detJ}})
     return 1.0
 end
 
-function get_integration_order(element::Poi1)
+function get_integration_order(::Poi1)
     return 1
 end
 
-@lintpragma("Ignore unused order")
-function get_integration_points(element::Poi1, order::Int64)
+function get_integration_points(::Poi1, ::Int64)
     return [ (1.0, (0.0, )) ]
 end
 
@@ -51,7 +50,7 @@ function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet
     return all(xi .>= 0.0) && (sum(xi) <= 1.0)
 end
 
-function get_reference_coordinates{B}(element::Element{B})
+function get_reference_coordinates{B}(::Element{B})
     return get_reference_element_coordinates(B)
 end