New function element_info!

src/FEMBase.jl

@@ -22,14 +22,15 @@ end
 export info, debug
 
 using FEMBasis
-using FEMBasis: AbstractBasis
+using FEMBasis: AbstractBasis, jacobian
 import FEMBasis: interpolate
 using FEMQuad: get_quadrature_points
 include("fields.jl")
 export interpolate, update!, DCTI, DCTV, DVTI, DVTV, CVTV, DVTId, DVTVd, field
 include("types.jl")
 include("sparse.jl")
 include("elements.jl")
+export element_info!
 include("elements_lagrange.jl")
 include("integrate.jl")
 include("problems.jl")

src/elements.jl

@@ -190,7 +190,7 @@ end
 
 function (element::Element)(ip, time, ::Type{Val{:Jacobian}})
     X = element("geometry", time)
-    J = FEMBasis.jacobian(element.properties, X, ip)
+    J = jacobian(element.properties, X, ip)
     return J
 end
 
@@ -219,26 +219,26 @@ function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:G
     return grad(element.properties, u, X, ip)
 end
 
-function (element::Element)(field_name::String, ip, time::Float64)
+function interpolate(element::Element, field_name, ip, time)
     field = element[field_name]
-    return interpolate(element, field, ip, time)
+    return interpolate_(element, field, ip, time)
 end
 
 @lintpragma("Ignore unused ip")
 @lintpragma("Ignore unused element")
-function interpolate(element::Element, field::DCTI, ip, time::Float64)
+function interpolate_(element::Element, field::DCTI, ip, time::Float64)
     return field.data
 end
 
-function interpolate(element::Element, field::DCTV, ip, time::Float64)
-    return interpolate(field, time)
+function interpolate_(element::Element, field::DCTV, ip, time::Float64)
+    return interpolate_(field, time)
 end
 
-function interpolate(element::Element, field::CVTV, ip, time::Float64)
+function interpolate_(element::Element, field::CVTV, ip, time::Float64)
     return field(ip, time)
 end
 
-function interpolate{F<:AbstractField}(element::Element, field::F, ip, time::Float64)
+function interpolate_{F<:AbstractField}(element::Element, field::F, ip, time::Float64)
     field_ = interpolate(field, time)
     basis = element(ip, time)
     n = length(element)
@@ -341,3 +341,14 @@ end
 function (element::Element)(field_name::String, time::Float64)
     return interpolate(element, field_name, time)
 end
+
+# element("displacement", (0.0, 0.0), 0.0)
+function (element::Element)(field_name::String, ip, time::Float64)
+    return interpolate(element, field_name, ip, time)
+end
+
+function element_info!{E,T}(bi::BasisInfo{E,T}, element::Element{E}, ip, time)
+    X = interpolate(element, "geometry", time)
+    eval_basis!(bi, X, ip)
+    return bi.J, bi.detJ, bi.N, bi.grad
+end

test/test_problems.jl

@@ -8,7 +8,7 @@ using FEMBase: get_assembly, get_elements
 using Base.Test
 
 import FEMBase: get_unknown_field_dimension, get_unknown_field_name
-import FEMBase: get_formulation_type
+import FEMBase: get_formulation_type, assemble!
 
 type P1 <: FieldProblem
     A :: Bool
@@ -129,3 +129,55 @@ end
     @test isapprox(e2("P1", (0.0,), 0.0), [0.0, 0.0])
     @test isapprox(e2("lambda", (0.0,), 0.0), [0.0, 0.0])
 end
+
+function assemble!{E}(problem::Problem{P1},
+                      assembly::Assembly,
+                      elements::Vector{Element{E}},
+                      time::Float64)
+
+    info("Assembling elements of kind $E")
+    bi = BasisInfo(E)
+    ndofs = length(E)
+    Ke = zeros(ndofs, ndofs)
+    K = assembly.K
+
+    for element in elements
+        fill!(Ke, 0.0)
+        for ip in get_integration_points(element)
+            J, detJ, N, dN = element_info!(bi, element, ip, time)
+            c = element("coefficient", ip, time)
+            Ke += ip.weight * c*dN'*dN * detJ
+        end
+        gdofs = get_gdofs(problem, element)
+        add!(K, gdofs, gdofs, Ke)
+    end
+
+    return nothing
+
+end
+
+@testset "test assemble test problem" begin
+    el1 = Element(Quad4, [1, 2, 3, 4])
+    el2 = Element(Tri3, [3, 2, 5])
+    X = Dict(1 => [0.0, 0.0],
+             2 => [1.0, 0.0],
+             3 => [1.0, 1.0],
+             4 => [0.0, 1.0],
+             5 => [2.0, 1.0])
+    elements = [el1, el2]
+    update!(elements, "geometry", X)
+    update!(elements, "coefficient", 6.0)
+
+    problem = Problem(P1, "test problem", 1)
+    add_elements!(problem, elements)
+    time = 0.0
+    assemble!(problem, time)
+    K = full(problem.assembly.K)
+    K_expected = [
+                  4.0 -1.0 -2.0 -1.0  0.0
+                 -1.0  7.0 -4.0 -2.0  0.0
+                 -2.0 -4.0 10.0 -1.0 -3.0
+                 -1.0 -2.0 -1.0  4.0  0.0
+                  0.0  0.0 -3.0  0.0  3.0]
+    @test isapprox(K, K_expected)
+end