From d606bde168aeb52654157844fd976685b2ec734f Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Wed, 25 May 2016 22:33:47 +0300
Subject: [PATCH] integration points

---
 src/JuliaFEM.jl                               |  13 +-
 src/dirichlet.jl                              |  36 ++--
 src/elasticity.jl                             | 103 ++++++-----
 src/elements.jl                               |  75 +++++---
 src/fields.jl                                 |  29 ++++
 src/integrate.jl                              |  36 ++--
 src/lagrange_macro.jl                         |   4 +-
 src/test.jl                                   | 162 +-----------------
 src/types.jl                                  |  88 ++++------
 test/test_define_new_element.jl               |  65 +++++++
 test/test_dirichlet.jl                        |  69 ++++++--
 ..._elasticity_2d_linear_with_surface_load.jl |  16 ++
 ...asticity_2d_nonlinear_with_surface_load.jl |  55 ++++--
 ...ticity_2d_plane_stress_stiffness_matrix.jl |   4 +-
 test/test_fields.jl                           |  25 ++-
 test/test_integration_points.jl               |  14 ++
 16 files changed, 424 insertions(+), 370 deletions(-)
 create mode 100644 test/test_define_new_element.jl
 create mode 100644 test/test_integration_points.jl

diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index b7d25ea..a085e25 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -15,16 +15,17 @@ autodiffcache = ForwardDiffCache()
 include("common.jl")
 
 include("fields.jl")
-export DCTI, Field
+export Field, DCTI, DVTI, DCTV, DVTV
+include("types.jl")  # data types: Point, IntegrationPoint, ...
+export AbstractPoint, Point, IntegrationPoint, IP, Node
 #include("basis.jl")  # interpolation of discrete fields
 #include("symbolic.jl") # a thin symbolic layer for fields
-#include("types.jl")  # type definitions
 
 ### ELEMENTS ###
 include("elements.jl") # common element routines
 export Node, AbstractElement, Element, update!, get_connectivity
 include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
-export Seg2, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10
+export Seg2, Seg3, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10
 
 #include("hierarchical.jl") # P-elements
 #include("mortar_elements.jl") # Mortar elements
@@ -48,7 +49,11 @@ export Dirichlet
 include("heat.jl")
 export Heat
 
-export assemble
+export assemble, assemble!
+
+function assemble!(problem::Problem, element::Element, time=0.0)
+    assemble!(problem.assembly, problem, element, time)
+end
 
 ### ASSEMBLY + SOLVE ###
 include("assembly.jl")
diff --git a/src/dirichlet.jl b/src/dirichlet.jl
index 956d56e..6ffb74b 100644
--- a/src/dirichlet.jl
+++ b/src/dirichlet.jl
@@ -30,20 +30,17 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     field_name = get_parent_field_name(problem)
     gdofs = get_gdofs(element, field_dim)
 
-  
-#   if problem.properties.formulation == :dual_basis
+    if problem.properties.dual_basis
         De, Me, Ae = get_dualbasis(element, time)
-#   else
-#       Ae = eye(nnodes)
-#       De = zeros(nnodes, nnodes)
-#       for (w, xi) in get_integration_points(element, Val{3})
-#           N = element(xi, time)
-#           detJ = element(xi, time, Val{:detJ})
-#           De += w*N'*N*detJ
-#       end
-#   end
-  
-#   De = Ae = eye(nnodes)
+    else
+        Ae = eye(nnodes)
+        De = zeros(nnodes, nnodes)
+        for ip in get_integration_points(element)
+            N = element(ip, time)
+            detJ = element(ip, time, Val{:detJ})
+            De += ip.weight*N'*N*detJ
+        end
+    end
 
     # left hand side
     for i=1:field_dim
@@ -55,20 +52,21 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     end
 
     # right hand side
-    for (w, xi) in get_integration_points(element, Val{3})
-        detJ = element(xi, time, Val{:detJ})
-        N = element(xi, time)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
 
         for i=1:field_dim
             ldofs = gdofs[i:field_dim:end]
             if haskey(element, field_name*" $i")
-                g = element(field_name*" $i", xi, time)
+                g = element(field_name*" $i", ip, time)
                 if true
                     haskey(element, "displacement") || continue
-                    g_prev = element(field_name, xi, time)
+                    g_prev = element(field_name, ip, time)
                     g -= g_prev[i]
                 end
-                add!(assembly.g, ldofs, w*g*Ae*N'*detJ)
+                add!(assembly.g, ldofs, w*g*Ae*N')
             end
         end
 
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 63a8052..e354190 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -1,55 +1,40 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-""" Concrete Elasticity type. """
+""" Elasticity problem
+
+
+"""
 type Elasticity <: FieldProblem
     # these are found from problem.properties for type Problem{Elasticity}
     formulation :: Symbol
     finite_strain :: Bool
-    use_forwarddiff :: Bool
 end
 function Elasticity()
     # formulations: plane_stress, plane_strain, continuum
-    return Elasticity(:continuum, true, false)
+    return Elasticity(:continuum, true)
 end
 
-# in case of experimenting new things;
-# 1. import JuliaFEM.Core: assemble!
-# 2. copy/paste assemble! code to notebook
-# 3. change to last argument, i.e. ::Type{Val{:plane_stress}} to ::Type{Val{:my_formulation}}
-# 4. when running code: set problem.properties.formulation = :my_formulation
-# 5. let multiple dispatch do the magic for you
-
-function get_unknown_field_name(::Type{Elasticity})
+function get_unknown_field_name(problem::Problem{Elasticity})
     return "displacement"
 end
 
 function get_formulation_type(problem::Problem{Elasticity})
     # we are solving residual and add increment to previous solution vector
     return :incremental
-    #return :total
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real)
+function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time=0.0)
     props = problem.properties
     gdofs = get_gdofs(problem, element)
-    if props.use_forwarddiff
-        Kt, f = assemble(problem, element, time, Val{:forwarddiff})
-    elseif props.formulation == :continuum
-        Kt, f = assemble(problem, element, time, Val{:continuum})
-    elseif (props.formulation == :plane_stress) || (props.formulation == :plane_strain)
+    if problem.properties.formulation in [:plane_stress, :plane_strain]
         Kt, f = assemble(problem, element, time, Val{:plane})
+    else
+        Kt, f = assemble(problem, element, time, Val{problem.properties.formulation})
     end
     add!(assembly.K, gdofs, gdofs, Kt)
     add!(assembly.f, gdofs, f)
-end
-
-function assemble(problem::Problem{Elasticity}, element::Element, time=0.0)
-    problem.properties
-    if problem.properties.formulation in [:plane_stress, :plane_strain]
-        return assemble(problem, element, time, Val{:plane})
-    end
-    return assemble(problem, element, time, Val{problem.properties.formulation})
+    return Kt, f
 end
 
 """ Elasticity equations for 2d cases. """
@@ -63,16 +48,17 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
     Kt = zeros(dim*nnodes, dim*nnodes)
     f = zeros(dim*nnodes)
 
-    for (w, xi) in get_integration_points(element)
+    for ip in get_integration_points(element)
 
-        detJ = element(xi, time, Val{:detJ})
-        N = element(xi, time)
-        dN = element(xi, time, Val{:Grad})
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        dN = element(ip, time, Val{:Grad})
 
         # kinematics; calculate deformation gradient and strain
         gradu = zeros(dim, dim)
         if haskey(element, "displacement")
-            gradu += element("displacement", xi, time, Val{:Grad})
+            gradu += element("displacement", ip, time, Val{:Grad})
         end
         strain = zeros(dim , dim)
         strain += 1/2*(gradu' + gradu)
@@ -84,8 +70,8 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
 
         # constitutive equations; material model (isotropic linear material here)
         # get_material(problem, element, ...)
-        E = element("youngs modulus", xi, time)
-        nu = element("poissons ratio", xi, time)
+        E = element("youngs modulus", ip, time)
+        nu = element("poissons ratio", ip, time)
         if props.formulation == :plane_stress
             D = E/(1.0 - nu^2) .* [
                 1.0  nu 0.0
@@ -99,9 +85,15 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
         else
             error("unknown plane formulation: $(props.formulation)")
         end
-
         # calculate stress
-        S = D*[strain[1,1]; strain[2,2]; 2*strain[1,2]]
+        strain_vec = [strain[1,1]; strain[2,2]; 2*strain[1,2]]
+        stress_vec = D*strain_vec
+        stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]]
+        cauchy_stress = F'*stress*F/det(F)
+        cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
+
+        update!(ip, "strain", time => strain_vec)
+        update!(ip, "stress", time => cauchy_stress)
 
         # add contributions: material and geometric stiffness + internal forces
         fill!(BL, 0.0)
@@ -121,29 +113,29 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
             BNL[4, 2*(i-1)+2] = dN[2,i]
         end
         S2 = zeros(2*dim, 2*dim)
-        S2[1,1] = S[1]
-        S2[2,2] = S[2]
-        S2[1,2] = S2[2,1] = S[3]
+        S2[1,1] = stress_vec[1]
+        S2[2,2] = stress_vec[2]
+        S2[1,2] = S2[2,1] = stress_vec[3]
         S2[3:4,3:4] = S2[1:2,1:2]
 
-        Kt += w*BL'*D*BL*detJ # material stiffness
+        Kt += w*BL'*D*BL # material stiffness
         if props.finite_strain # add geometric stiffness
-            Kt += w*BNL'*S2*BNL*detJ # geometric stiffness
+            Kt += w*BNL'*S2*BNL # geometric stiffness
         end
 
         if get_formulation_type(problem) == :incremental
-            f -= w*BL'*S*detJ # internal force
+            f -= w*BL'*stress_vec # internal force
         end
 
         # volume load
         if haskey(element, "displacement load")
-            b = element("displacement load", xi, time)
-            f += w*vec(N'*b)*detJ
+            b = element("displacement load", ip, time)
+            f += w*vec(N'*b)
         end
         for i=1:dim
             if haskey(element, "displacement load $i")
-                b = element("displacement load $i", xi, time)
-                f[i:dim:end] += w*vec(b*N)*detJ
+                b = element("displacement load $i", ip, time)
+                f[i:dim:end] += w*vec(b*N)
             end
         end
 
@@ -160,29 +152,30 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E
     Kt = zeros(dim*nnodes, dim*nnodes)
     f = zeros(dim*nnodes)
 
-    for (w, xi) in get_integration_points(element)
+    for ip in get_integration_points(element)
 
-        detJ = element(xi, time, Val{:detJ})
-        N = element(xi, time)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
 
         if haskey(element, "displacement traction force")
-            T = element("displacement traction force", xi, time)
-            f += w*vec(T*N)*detJ
+            T = element("displacement traction force", ip, time)
+            f += w*vec(T*N)
         end
 
         for i=1:dim
             # traction force for ith component
             if haskey(element, "displacement traction force $i")
-                T = element("displacement traction force $i", xi, time)
-                f[i:dim:end] += w*vec(T*N)*detJ
+                T = element("displacement traction force $i", ip, time)
+                f[i:dim:end] += w*vec(T*N)
             end
         end
 
         if haskey(element, "nt displacement traction force")
             # traction force given in normal-tangential direction
-            T = element("nt displacement traction force", xi, time)
-            Q = element("normal-tangential coordinates", xi, time)
-            f += w*vec(Q'*T*N)*detJ
+            T = element("nt displacement traction force", ip, time)
+            Q = element("normal-tangential coordinates", ip, time)
+            f += w*vec(Q'*T*N)
         end
 
     end
diff --git a/src/elements.jl b/src/elements.jl
index 676a0ca..1142470 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -3,17 +3,19 @@
 
 abstract AbstractElement
 
-typealias Node Vector{Float64}
-
 type Element{E<:AbstractElement}
     id :: Int
     connectivity :: Vector{Int}
+    integration_points :: Vector{IP}
     fields :: Dict{ASCIIString, Field}
     properties :: E
 end
 
 function Element{E<:AbstractElement}(::Type{E}, connectivity=[], id=-1, fields=Dict(), properties...)
-    Element{E}(id, connectivity, fields, E(properties...))
+    variant = E(properties...)
+    ips = get_integration_points(variant)
+    integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
+    Element{E}(id, connectivity, integration_points, fields, variant)
 end
 
 function getindex(element::Element, field_name::ASCIIString)
@@ -28,26 +30,19 @@ function call(element::Element, field_name::ASCIIString, time=0.0)
     return element[field_name](time)
 end
 
-function call(element::Element, xi::Vector, time=0.0)
-    get_basis(element, xi, time)
+function call(element::Element, ip, time=0.0)
+    get_basis(element, ip, time)
 end
 
-function call(element::Element, field_name::ASCIIString, xi::Vector, time=0.0)
-    field = element[field_name](time)
-    isa(field, DCTI) && return field.data
-    basis = element(xi, time)
-    return basis*field
-end
-
-function call(element::Element, xi::Vector, time, ::Type{Val{:Jacobian}})
+function call(element::Element, ip, time, ::Type{Val{:Jacobian}})
     X = element["geometry"](time)
-    dN = get_dbasis(element, xi, time)
+    dN = get_dbasis(element, ip, time)
     J = sum([kron(dN[:,i], X[i]') for i=1:length(X)])
     return J
 end
 
-function call(element::Element, xi::Vector, time, ::Type{Val{:detJ}})
-    J = element(xi, time, Val{:Jacobian})
+function call(element::Element, ip, time, ::Type{Val{:detJ}})
+    J = element(ip, time, Val{:Jacobian})
     n, m = size(J)
     if n == m  # volume element
         return det(J)
@@ -60,13 +55,23 @@ function call(element::Element, xi::Vector, time, ::Type{Val{:detJ}})
     end
 end
 
-function call(element::Element, xi::Vector, time, ::Type{Val{:Grad}})
-    J = element(xi, time, Val{:Jacobian})
-    return inv(J)*get_dbasis(element, xi, time)
+function call(element::Element, ip, time, ::Type{Val{:Grad}})
+    J = element(ip, time, Val{:Jacobian})
+    return inv(J)*get_dbasis(element, ip, time)
 end
 
-function call(element::Element, field_name, xi::Vector, time, ::Type{Val{:Grad}})
-    element(xi, time, Val{:Grad})*element[field_name](time)
+function call(element::Element, field_name::ASCIIString, ip, time, ::Type{Val{:Grad}})
+    element(ip, time, Val{:Grad})*element[field_name](time)
+end
+
+function call(element::Element, field_name::ASCIIString, ip, time=0.0)
+    field = element[field_name](time)
+    isa(field, DCTI) && return field.data
+    basis = element(ip, time)
+    n = length(element)
+    m = length(field)
+    @assert n == m
+    return sum([field[i]*basis[i] for i=1:n])
 end
 
 #function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0)
@@ -122,6 +127,9 @@ function get_dbasis(element::Element, xi::Vector, time)
     basis(xi) = vec(get_basis(element, xi, time))
     return ForwardDiff.jacobian(basis, xi, cache=autodiffcache)'
 end
+function get_dbasis(element::Element, ip::IP, time)
+    get_dbasis(element, ip.coords, time)
+end
 
 """ Check existence of field. """
 function haskey(element::Element, field_name)
@@ -132,6 +140,20 @@ function get_connectivity(element::Element)
     return element.connectivity
 end
 
+function get_integration_points(element::Element)
+    return element.integration_points
+end
+
+""" This is a special case, temporarily change order
+of integration scheme mainly for mass matrix.
+"""
+function get_integration_points(element::Element, change_order::Int)
+    order = get_integration_order(element.properties)
+    order += change_order
+    ips = get_integration_points(element.properties, Val{order})
+    return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
+end
+
 function get_gdofs(element::Element)
     return get_gdofs(element, 1)
 end
@@ -141,11 +163,12 @@ function get_dualbasis(element::Element, time)
     nnodes = length(element)
     De = zeros(nnodes, nnodes)
     Me = zeros(nnodes, nnodes)
-    for (w, xi) in get_integration_points(element, Val{3})
-        detJ = element(xi, time, Val{:detJ})
-        N = element(xi, time)
-        De += w*diagm(vec(N))*detJ
-        Me += w*N'*N*detJ
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        De += w*diagm(vec(N))
+        Me += w*N'*N
     end
     return De, Me, De*inv(Me)
 end
diff --git a/src/fields.jl b/src/fields.jl
index 3401b26..38fb666 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -270,6 +270,35 @@ function Base.done(f::DVTI, s)
     return s > length(f.data)
 end
 
+""" Update time-dependent fields with new values.
+
+Examples
+--------
+
+julia> f = Field(0.0 => 1.0)
+julia> update!(f, 1.0 => 2.0)
+
+Now field has two (time, value) pairs: (0.0, 1.0) and (1.0, 2.0)
+
+Notes
+-----
+Time vector is assumed to be ordered t_i-1 < t_i < t_i+1. If updating
+field with already existing time the old value is replaced with new one.
+
+"""
+function update!{T}(field::Union{DCTV, DVTV}, val::Pair{Float64, T})
+    time, data = val
+    if isapprox(last(field).time, time)
+        last(field).data = data
+    else
+        push!(field.data, Increment(val...))
+    end
+end
+
+function update!{T}(field::Union{DCTI, DVTI}, val::T)
+    field.data = val
+end
+
 ### Accessing continuous fields
 
 function Base.call(field::CVTI, xi::Vector)
diff --git a/src/integrate.jl b/src/integrate.jl
index 345498b..56e3b60 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -51,9 +51,9 @@ end
 
 ### 1d elements
 
-typealias CartesianLineElement Union{Element{Seg2}, Element{Seg3}}
-typealias CartesianSurfaceElement Union{Element{Quad4}}
-typealias CartesianVolumeElement Union{Element{Hex8}}
+typealias CartesianLineElement Union{Seg2, Seg3}
+typealias CartesianSurfaceElement Union{Quad4}
+typealias CartesianVolumeElement Union{Hex8}
 
 function get_integration_points(element::CartesianLineElement, ::Type{Val{1}})
     w, xi = get_integration_points(Val{1})
@@ -100,7 +100,7 @@ end
 # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
 # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
 
-typealias TriangularElement Union{Element{Tri3}, Element{Tri6}}
+typealias TriangularElement Union{Tri3, Tri6}
 
 function get_integration_points(element::TriangularElement, ::Type{Val{1}})
     weights = [0.5]
@@ -152,7 +152,7 @@ end
 
 ### 3d elements
 
-typealias TetrahedralElement Union{Element{Tet4}, Element{Tet10}}
+typealias TetrahedralElement Union{Tet4, Tet10}
 
 function get_integration_points(element::TetrahedralElement, ::Type{Val{1}})
     weights = 1.0/6.0*[1.0]
@@ -191,23 +191,25 @@ end
 ### default number of integration points for each element
 ### 2 for linear elements, 3 for quadratic
 
-typealias LinearElement Union{
-    Element{Seg2},
-    Element{Tri3},
-    Element{Quad4},
-    Element{Tet4},
-    Element{Hex8}}
+typealias LinearElement Union{Seg2, Tri3, Quad4, Tet4, Hex8}
 
-typealias QuadraticElement Union{
-    Element{Seg3},
-    Element{Tri6},
-    Element{Tet10}}
+typealias QuadraticElement Union{Seg3, Tri6, Tet10}
 
-function get_integration_points(element::LinearElement; order=2)
+function get_integration_order(element::LinearElement)
+    return 2
+end
+
+function get_integration_order(element::QuadraticElement)
+    return 3
+end
+
+function get_integration_points(element::LinearElement)
+    order = get_integration_order(element)
     get_integration_points(element, Val{order})
 end
 
-function get_integration_points(element::QuadraticElement; order=2)
+function get_integration_points(element::QuadraticElement)
+    order= get_integration_order(element)
     get_integration_points(element, Val{order})
 end
 
diff --git a/src/lagrange_macro.jl b/src/lagrange_macro.jl
index 32bddf9..da4b4be 100644
--- a/src/lagrange_macro.jl
+++ b/src/lagrange_macro.jl
@@ -42,8 +42,8 @@ macro create_lagrange_element(element_name, element_description, X, P)
         A = calculate_lagrange_basis_coefficients($P, $X)
         #basis(xi) = C*$P(xi)
 
-        function get_basis(element::Element{$eltype}, xi::Vector, time)
-            return transpose($P(xi))*A
+        function get_basis(element::Element{$eltype}, ip, time)
+            return transpose($P(ip))*A
         end
 
         function size(element::Element{$eltype})
diff --git a/src/test.jl b/src/test.jl
index 5712592..2dfd955 100644
--- a/src/test.jl
+++ b/src/test.jl
@@ -1,167 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-if VERSION >= v"0.5-"
+if VERSION >= v"0.5-dev+7720"
     using Base.Test
 else
     using BaseTestNext
+    const Test = BaseTestNext
 end
 
-
-abstract TestResult
-
-type NormalTestResult <: TestResult
-    test_function :: Function
-    result
-end
-
-type CriticalTestResult <: TestResult
-    filename :: ASCIIString
-    message
-end
-
-function TestResult(test_function::Function, result)
-    NormalTestResult(test_function, result)
-end
-
-function TestResult(filename::ASCIIString, message)
-    CriticalTestResult(filename, message)
-end
-
-global test_results = []
-
-    
-function get_test_functions(func)
-    return Function[]
-end
-
-""" Return all functions from module with name starting test """
-function get_test_functions(mod::Module)
-    test_function_names = filter((k) -> startswith(string(k), "test_"), names(mod, true))
-    if haskey(ENV, "JULIAFEM_TEST_SLOW")
-        info("JULIAFEM_TEST_SLOW set, testing also tests that are taking a long time")
-        slow_test_functions = filter((k) -> startswith(string(k), "slow_test_"), names(mod, true))
-        append!(test_function_names, slow_test_functions)
-    end
-    test_function_expressions = map((k) -> :($mod.$k), test_function_names)
-    test_functions = map(eval, test_function_expressions)
-    return test_functions
-end
-
-""" Run tests from some file. """
-function run_test(filename::ASCIIString, test_function=nothing)
-    info("running tests from $filename")
-    test_module = nothing
-    try
-        test_module = include(filename)
-    catch error
-        warn("Unable to include file $filename for testing.")
-        err = Base.showerror(Base.STDOUT, error)
-        push!(test_results, TestResult(filename, "Unable to include file: $error"))
-        return
-    end
-    if isa(test_function, Void)
-        test_functions = get_test_functions(test_module)
-    else
-        test_functions = [eval( :($test_module.$test_function) )]
-    end
-    if length(test_functions) == 0
-        warn("Unable to get test functions for file $filename. Define test functions inside module, look for test_heat.jl for concrete example how to do that.")
-        push!(test_results, TestResult(filename, "Unable to find test functions"))
-        return
-    end
-    for test_function in test_functions
-        run_test(test_function)
-    end
-end
-
-""" Run single test set. """
-function run_test(test_function::Function)
-
-    allok = true
-
-    function test_handler(r::Base.Test.Success)
-        print(".")
-        result = TestResult(test_function, r)
-        push!(test_results, result)
-    end 
-
-    function test_handler(r::Base.Test.Failure)
-        allok = false
-        print(" [\x1b[31mFAIL\x1b[0m]")
-        push!(test_results, TestResult(test_function, r))
-        println()
-        println("Test failed: $(r.expr)")
-        #warn("partially evaluated expression: $(r.resultexpr)")
-    end 
-
-    function test_handler(r::Base.Test.Error)
-        allok = false
-        println(" [\x1b[31mERROR\x1b[0m]")
-        push!(test_results, TestResult(test_function, r))
-        println("Error when testing: $(r.expr)")
-    end 
-
-    Base.Test.with_handler(test_handler) do
-        print("TEST: $test_function() ")
-        # TODO: print docstring of test function if defined.
-        #@doc($test_function)
-        try
-            test_function()
-        catch error
-            allok = false
-            println("[\x1b[31mCRITICAL\x1b[0m]")
-            println("Testing $test_function() stopped for critical error:\n$error")
-            Base.showerror(Base.STDOUT, error)
-            println()
-            println("Cannot continue to test function $test_function()")
-            push!(test_results, TestResult("$test_function", error))
-        end
-
-        if allok
-            println(" [\x1b[32mPASS\x1b[0m]")
-        else
-            filename, linenum = Base.functionloc(test_function)
-            println("Test $test_function() failed: file $filename, line $linenum.")
-            println()
-        end
-    end
-
-end
-
-function print_test_statistics()
-    info("################# TEST RESULTS #####################")
-    passed = 0
-    failed = 0
-    errors = 0
-    critical = 0
-    for result in test_results
-        if isa(result, NormalTestResult)
-            if isa(result.result, Base.Test.Success)
-                passed += 1
-                continue
-            elseif isa(result.result, Base.Test.Failure)
-                failed += 1
-                filename, linenum = Base.functionloc(result.test_function)
-                functionname = Base.function_name(result.test_function)
-                warn("test failed: $(result.result.expr), partially evaluated expression: $(result.result.expr)")
-                warn("in function $functionname, file $filename, line $linenum")
-                continue
-            elseif isa(result.result, Base.Test.Error)
-                errors += 1
-                filename, linenum = Base.functionloc(result.test_function)
-                functionname = Base.function_name(result.test_function)
-                warn("error in function: $(result.result.expr)")
-                warn("in function $functionname, file $filename, line $linenum")
-                continue
-            end
-        elseif isa(result, CriticalTestResult)
-            critical += 1
-            warn("Critical error on $(result.filename): $(result.message)")
-        end
-    end
-    info("$passed test passed, $failed test failed, $errors errors, $critical critical failures")
-    return passed, failed, errors, critical
-end
-
-export @test, @testset, @test_throws, run_test, print_test_statistics
+export @test, @testset, @test_throws
diff --git a/src/types.jl b/src/types.jl
index 4b4683c..658dc8d 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -3,77 +3,63 @@
 
 typealias Node Vector{Float64}
 
+abstract AbstractPoint
 
-"""
-Integration point
-
-xi
-    (dimensionless) coordinates of integration point
-weight
-    integration weight
-fields
-    FieldSet what can be used to store internal variables, stress, strain, ...
-"""
-immutable IntegrationPoint
-    xi :: Vector
+type Point{P<:AbstractPoint}
+    id :: Int
     weight :: Float64
+    coords :: Vector{Float64}
     fields :: Dict{ASCIIString, Field}
-    changed :: Bool
+    properties :: P
 end
 
-function IntegrationPoint(xi, weight)
-    return IntegrationPoint(xi, weight, FieldSet(), false)
+function setindex!{T}(point::Point, val::Pair{Float64, T}, field_name::ASCIIString)
+    point.fields[field_name] = Field(val)
 end
 
-function setindex!{T<:ForwardDiff.ForwardDiffNumber}(ip::IntegrationPoint, data::Array{T,2}, field_name::ASCIIString)
-    data = ForwardDiff.get_value(data)
-    setindex!(ip, data, field_name)
-end
-function setindex!(ip::IntegrationPoint, data, field_name)
-    ip.fields[field_name] = Field(data)
-    ip.changed = true
+function getindex(point::Point, field_name::ASCIIString)
+    return point.fields[field_name]
 end
 
-function getindex(ip::IntegrationPoint, field_name::ASCIIString)
-    ip.fields[field_name]
+function getindex(point::Point, idx::Int)
+    return point.coords[idx]
 end
 
-function convert(::Type{Number}, ip::IntegrationPoint)
-    return ip.xi
+function haskey(point::Point, field_name::ASCIIString)
+    return haskey(point.fields, field_name)
 end
 
-function call(field::CVTI, ip::IntegrationPoint)
-    return call(field, ip.xi)
+function call(point::Point, field_name::ASCIIString, time::Float64=0.0)
+    point.fields[field_name](time).data
 end
 
-function call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
-    call(basis, field, ip.xi)
+function update!{T}(point::Point, field_name, val::Pair{Float64, T})
+    if haskey(point, field_name)
+        update!(point[field_name], val)
+    else
+        point[field_name] = val
+    end
 end
 
-function call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
-    call(basis, field, ip.xi, Val{:grad})
+#= TODO: in future
+type Node <: AbstractPoint
 end
 
-function call(basis::CVTI, field::DVTI, ip::IntegrationPoint)
-    call(basis, field, ip.xi)
+type MaterialPoint <: AbstractPoint
+end
+=#
+
+type IntegrationPoint <: AbstractPoint
 end
 
-function call(basis::CVTI, geometry::DVTI, field::Union{DCTI, DVTI}, ip::IntegrationPoint, ::Type{Val{:grad}})
-    call(basis, geometry, field, ip.xi, Val{:grad})
+typealias IP Point{IntegrationPoint}
+
+function IP(id, weight, coords)
+    return IP(id, weight, coords, Dict(), IntegrationPoint())
+end
+
+function convert(::Type{IP}, data::Tuple{Float64, Vector{Float64}})
+    weight, coords = data
+    return IP(-1, weight, coords)
 end
 
-#function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint)
-#    return call(basis, increment, ip.xi)
-#end
-#function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
-#    return call(basis, increment, ip.xi, Val{:grad})
-#end
-#function Base.call(basis::Basis, field::Field, ip::IntegrationPoint, ::Type{Val{:grad}})
-#    return call(basis, field, ip.xi, Val{:grad})
-#end
-#function Base.call(basis::Basis, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
-#    return call(basis, geometry, field, ip.xi, Val{:grad})
-#end
-#function Base.call(basis::Basis, field::Field, ip::IntegrationPoint)
-#    return call(basis, field, ip.xi)
-#end
diff --git a/test/test_define_new_element.jl b/test/test_define_new_element.jl
new file mode 100644
index 0000000..39bde9d
--- /dev/null
+++ b/test/test_define_new_element.jl
@@ -0,0 +1,65 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+importall Base
+import JuliaFEM: get_basis, get_dbasis, get_integration_points
+
+type MyQuad4 <: AbstractElement
+end
+
+function get_basis(element::Element{MyQuad4}, xi, time)
+    1/4*[(1-xi[1])*(1-xi[2]) (1+xi[1])*(1-xi[2]) (1+xi[1])*(1+xi[2]) (1-xi[1])*(1+xi[2])]
+end
+
+function get_dbasis(element::Element{MyQuad4}, xi, time)
+    1/4*[-(1-xi[2])    (1-xi[2])   (1+xi[2])  -(1+xi[2])
+         -(1-xi[1])   -(1+xi[1])   (1+xi[1])   (1-xi[1])]
+end
+
+function get_integration_points(element::MyQuad4)
+    [
+        (1.0, 1.0/sqrt(3.0)*[-1, -1]),
+        (1.0, 1.0/sqrt(3.0)*[ 1, -1]),
+        (1.0, 1.0/sqrt(3.0)*[ 1,  1]),
+        (1.0, 1.0/sqrt(3.0)*[-1,  1])
+    ]
+end
+
+function length(element::Element{MyQuad4})
+    return 4
+end
+
+function size(element::Element{MyQuad4})
+    return (2, 4)
+end
+
+@testset "test new element" begin
+    el = Element(MyQuad4)
+    el["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
+    el["displacement"] = Vector{Float64}[[0.0,0.0], [0.0,0.0], [1.0,0.0], [0.0,0.0]]
+    @test isapprox(el("geometry", [0.0, 0.0]), [0.5, 0.5])
+    @test isapprox(el("displacement", [0.0, 0.0], 0.0), [0.25, 0.0])
+    el["temperature thermal conductivity"] = 6.0
+    dim = length(el)
+    K = zeros(dim, dim)
+    A = 0.0
+    time = 0.0
+    for ip in get_integration_points(el)
+        dN = el(ip, time, Val{:Grad})
+        detJ = el(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        c = el("temperature thermal conductivity", ip, time)
+        K += w*c*dN'*dN
+        A += w
+    end
+    @test isapprox(A, 1.0)
+    K_expected = [
+        4.0  -1.0  -2.0  -1.0
+       -1.0   4.0  -1.0  -2.0
+       -2.0  -1.0   4.0  -1.0
+       -1.0  -2.0  -1.0   4.0]
+    @test isapprox(K, K_expected)
+end
+
diff --git a/test/test_dirichlet.jl b/test/test_dirichlet.jl
index d3a7de4..59ddd42 100644
--- a/test/test_dirichlet.jl
+++ b/test/test_dirichlet.jl
@@ -1,27 +1,66 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+using JuliaFEM
 using JuliaFEM.Test
-using JuliaFEM.Core: Tri3, Seg2, Dirichlet, Assembly, assemble!, Node, Problem
 
 #=
+In [36]: C = Matrix([[0], [30], [15]]) # node coordinates
+In [37]: A = Matrix([P.subs({x: C[i,0]}).T for i in range(len(P))])
+In [38]: N = P.T*A.inv()
+In [39]: Me = integrate(N.T*N, (x, 0, 30))
+In [40]: De = diag(*integrate(N, (x, 0, 30)))
+In [41]: Me
+Out[41]: 
+Matrix([
+[ 4, -1,  2],
+[-1,  4,  2],
+[ 2,  2, 16]])
+In [42]: De
+Out[42]: 
+Matrix([
+[5, 0,  0],
+[0, 5,  0],
+[0, 0, 20]])
+=#
+
 @testset "dirichlet problem in 1 dimension" begin
-    element = Seg2([1, 2])
-    element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
-    element["temperature"] = 0.0
-    problem = Problem(Dirichlet, "test problem", 1, "temperature")
-    push!(problem, element)
-    assemble!(problem, 0.0)
-    C1 = full(problem.assembly.C1)
-    info("C1")
-    dump(C1)
-    C2 = full(problem.assembly.C2)
-    g = full(problem.assembly.g)
+    element = Element(Seg2, [1, 2])
+    element["geometry"] = Vector{Float64}[[0.0, 0.0], [6.0, 0.0]]
+    element["temperature 1"] = 0.0
+    p1 = Problem(Dirichlet, "test problem 1", 1, "temperature")
+    p1.properties.dual_basis = false
+    p2 = Problem(Dirichlet, "test problem 2", 1, "temperature")
+    assemble!(p1, element)
+    assemble!(p2, element)
+    C1 = full(p1.assembly.C1)
+    C2 = full(p1.assembly.C2)
     @test isapprox(C1, C2)
-    @test isapprox(C1, 1/6*[2 1; 1 2])
-    @test isapprox(g, [0.0, 0.0])
+    @test isapprox(C1, [2.0 1.0; 1.0 2.0])
+    C1 = full(p2.assembly.C1)
+    C2 = full(p2.assembly.C2)
+    @test isapprox(C1, C2)
+    @test isapprox(C1, [3.0 0.0; 0.0 3.0])
+
+    element = Element(Seg3, [1, 2, 3])
+    element["geometry"] = Vector{Float64}[[0.0, 0.0], [30.0, 0.0], [15.0, 0.0]]
+    element["temperature 1"] = 0.0
+    p1 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
+    p1.properties.dual_basis = false
+    p2 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
+    assemble!(p1, element)
+    assemble!(p2, element)
+    C1 = full(p1.assembly.C1)
+    C2 = full(p1.assembly.C2)
+    @test isapprox(C1, C2)
+    @test isapprox(C1, [4.0 -1.0 2.0; -1.0 4.0 2.0; 2.0 2.0 16.0])
+    C1 = full(p2.assembly.C1)
+    C2 = full(p2.assembly.C2)
+    @test isapprox(C1, C2)
+    @test isapprox(C1, [5.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 20.0])
 end
 
+#=
 @testset "dirichlet problem using tri3 surface element" begin
     element = Tri3([1, 2, 3])
     element["geometry"] = Node[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
@@ -34,7 +73,6 @@ end
     @test isapprox(C1, C2)
     @test isapprox(C1, 1/24*[2 1 1; 1 2 1; 1 1 2])
 end
-=#
 
 @testset "dirichlet problem in 2 dimensions" begin
     element = Seg2([1, 2])
@@ -72,4 +110,5 @@ end
     @test isapprox(C1, C1_expected)
     @test isapprox(g, [0.0, 0.0, 0.0, 0.0])
 end
+=#
 
diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl
index d1e2e03..d7f056e 100644
--- a/test/test_elasticity_2d_linear_with_surface_load.jl
+++ b/test/test_elasticity_2d_linear_with_surface_load.jl
@@ -35,6 +35,7 @@ using JuliaFEM.Test
     solver = Solver("solve block problem")
     push!(solver, block, bc_sym)
     call(solver)
+
     f = 288.0
     g = 576.0
     E = 288.0
@@ -43,4 +44,19 @@ using JuliaFEM.Test
     u3 = reshape(block.assembly.u, 2, 4)[:,3]
     info("u3 = $u3")
     @test isapprox(u3, u3_expected)
+
+    info("strain")
+    for ip in get_integration_points(elements[1])
+        eps = ip("strain")
+        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3]
+        @test isapprox(eps, [u3; 0.0])
+    end
+
+    info("stress")
+    for ip in get_integration_points(elements[1])
+        sig = ip("stress")
+        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
+        @test isapprox(sig, [0.0; g; 0.0])
+    end
+
 end
diff --git a/test/test_elasticity_2d_nonlinear_with_surface_load.jl b/test/test_elasticity_2d_nonlinear_with_surface_load.jl
index 2370740..b6126f1 100644
--- a/test/test_elasticity_2d_nonlinear_with_surface_load.jl
+++ b/test/test_elasticity_2d_nonlinear_with_surface_load.jl
@@ -8,15 +8,21 @@ using JuliaFEM.Test
 @testset "test 2d nonlinear elasticity with surface load" begin
     meshfile = "/geometry/2d_block/BLOCK_1elem.med"
     mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile)
+
     # field problem
-    body = Problem(Elasticity, "BLOCK", 2)
-    body.properties.formulation = :plane_stress
-    body_elements = aster_create_elements(mesh, :BLOCK, :QU4)
-    update!(body_elements, "youngs modulus", 900.0)
-    update!(body_elements, "poissons ratio", 0.25)
-    trac_elements = aster_create_elements(mesh, :TOP, :SE2)
-    update!(trac_elements, "displacement traction force 2", -100.0)
-    push!(body, body_elements..., trac_elements...)
+    block = Problem(Elasticity, "BLOCK", 2)
+    block.properties.formulation = :plane_stress
+
+    elements = aster_create_elements(mesh, :BLOCK, :QU4)
+    update!(elements, "youngs modulus", 288.0)
+    update!(elements, "poissons ratio", 1/3)
+    update!(elements, "displacement load 2", 576.0)
+    push!(block, elements...)
+
+    traction = aster_create_elements(mesh, :TOP, :SE2)
+    update!(traction, "displacement traction force 2", 288.0)
+    push!(block, traction...)
+
     # boundary conditions
     bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement")
     bc_elements_left = aster_create_elements(mesh, :LEFT, :SE2)
@@ -24,11 +30,34 @@ using JuliaFEM.Test
     update!(bc_elements_left, "displacement 1", 0.0)
     update!(bc_elements_bottom, "displacement 2", 0.0)
     push!(bc_sym, bc_elements_left..., bc_elements_bottom...)
+
     solver = Solver("solve block problem")
-    push!(solver, body, bc_sym)
+    push!(solver, block, bc_sym)
     call(solver)
-    # result is verified using code aster
-    u3_expected = [3.17431158889468E-02, -1.38591518927826E-01]
-    u3 = reshape(body.assembly.u, 2, 4)[:,3]
-    @test isapprox(u3, u3_expected)
+
+    # from code aster
+    u3_expected = [-4.92316106779943E-01, 7.96321884292103E-01]
+    eps_zz = -3.71128811855451E-01
+    eps_expected = [-3.71128532282463E-01, 1.11338615599337E+00, 0.0]
+    sig_expected = [ 3.36174888827909E-05, 2.23478729403118E+03, 0.0]
+
+    u3 = reshape(block.assembly.u, 2, 4)[:, 3]
+    info("u3 = $u3")
+    @test isapprox(u3, u3_expected, atol=1.0e-5)
+
+    info("strain")
+    for ip in get_integration_points(elements[1])
+        eps = ip("strain")
+        #eps = [eps[1,1]; eps[2,2]; eps[1,2]]
+        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3]
+        @test isapprox(eps, eps_expected)
+    end
+
+    info("stress")
+    for ip in get_integration_points(elements[1])
+        sig = ip("stress")
+        #sig = [sig[1,1]; sig[2,2]; sig[1,2]]
+        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
+        @test isapprox(sig, sig_expected)
+    end
 end
diff --git a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl
index e42cd1d..1f5827b 100644
--- a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl
+++ b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl
@@ -7,7 +7,7 @@ using JuliaFEM
 using JuliaFEM.Test
 
 @testset "test 2d linear elasticity local matrices" begin
-    element = Element(Quad4)
+    element = Element(Quad4, [1, 2, 3, 4])
     element["geometry"] = Vector{Float64}[
         [0.0, 0.0],
         [1.0, 0.0],
@@ -19,7 +19,7 @@ using JuliaFEM.Test
 
     problem = Problem(Elasticity, "[0x1] x [0x1] block", 2)
     problem.properties.formulation = :plane_stress
-    K, f = assemble(problem, element)
+    K, f = assemble!(problem, element)
 
     K_expected = [
         144  54 -90   0 -72 -54  18   0
diff --git a/test/test_fields.jl b/test/test_fields.jl
index e0b93ee..ea83b14 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -1,15 +1,26 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-
-module FieldTests
-
+using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM.Core: Field
+@testset "test updating time dependent fields" begin
+    f = Field(0.0 => 1.0)
+    @test last(f).time == 0.0
+    @test last(f).data == 1.0
+    update!(f, 0.0 => 2.0)
+    @test last(f).time == 0.0
+    @test last(f).data == 2.0
+    @test length(f) == 1
+    update!(f, 1.0 => 3.0)
+    @test last(f).time == 1.0
+    @test last(f).data == 3.0
+    @test length(f) == 2
+end
 
-function test_create_field()
+@testset "test updating time invariant fields" begin
     f = Field(1.0)
-end
-
+    @test f.data == 1.0
+    update!(f, 2.0)
+    @test f.data == 2.0
 end
diff --git a/test/test_integration_points.jl b/test/test_integration_points.jl
new file mode 100644
index 0000000..d7d9f84
--- /dev/null
+++ b/test/test_integration_points.jl
@@ -0,0 +1,14 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "test integration point" begin
+    ip = IP(1, 1.0, sqrt(1.0/3.0)*[-1.0, -1.0])
+    strain = [1.0 2.0; 3.0 4.0]
+    update!(ip, "strain", 0.0 => strain)
+    @test isapprox(ip("strain", 0.0), strain)
+    @test isapprox(ip("strain"), strain)
+end
+