Gets better all the time.

src/elements.jl

@@ -3,13 +3,16 @@
 
 abstract Element
 
+get_element(eq::Equation) = eq.element
+
 """
 Get jacobian of element evaluated at point xi
 """
 function get_jacobian(el::Element, xi)
-    dbasisdxi(xi) = get_dbasisdxi(el, xi)
+    dbasisdxi = get_dbasisdxi(el, xi)
     X = get_field(el, "coordinates")
-    J = interpolate(X, dbasisdxi, xi)'
+    #J = interpolate(X, dbasisdxi, xi)'
+    J = X*dbasisdxi
     return J
 end
 
@@ -103,6 +106,10 @@ type Point1 <: CG
     node_ids :: Array{Int, 1}
     fields :: Dict{ASCIIString, Any}
 end
+function Point1(node_ids)
+    fields = Dict{ASCIIString, Any}()
+    Point1(node_ids, fields)
+end
 
 # 1d Lagrange elements
 

src/math.jl

@@ -1,9 +1,7 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-"""
-This module contains math stuff, including interpolation, integration, linearization, ...
-"""
+## This module contains math stuff, including interpolation, integration, linearization, ...
 
 using ForwardDiff
 
@@ -185,3 +183,63 @@ function integrate!(f::Function, el::Element, target)
     end
 end
 
+get_integration_points(eq::Equation) = eq.integration_points
+
+"""
+Integrate f over element using Gaussian quadrature rules.
+
+Parameters
+----------
+el::Element
+    well defined element
+f::Function
+    Function to integrate
+"""
+function integrate(eq::Equation, f::Function)
+    target = []
+    for ip in get_integration_points(eq)
+        J = get_jacobian(eq.element, ip.xi)
+        push!(target, ip.weight*f(eq, ip)*det(J))
+    end
+    return sum(target)
+end
+
+
+"""
+Evaluate field in point xi using basis functions.
+"""
+function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1})
+    f = get_field(el, field)
+    if !isa(f, Array)
+        # This is scalar, nothing to interpolate
+        return f
+    end
+    basis = get_basis(el, xi)
+    dim, nnodes = size(f)
+    result = zeros(dim)
+    for i=1:nnodes
+        result += basis[i]*f[:,i]
+    end
+    return result
+end
+
+
+function linearize(eq::Equation, f::Function, field::ASCIIString)
+    function jacobian(eq::Equation, args...)
+        el = get_element(eq)
+        fld = get_field(el, field)
+        dim, nnodes = size(fld)
+        function helper(x::Vector)
+            orig = copy(fld)
+            set_field(el, field,  reshape(x, dim, nnodes))
+            y = f(eq, args...)
+            set_field(el, field, orig)
+            return y[:]
+        end
+        jac = ForwardDiff.jacobian(helper)
+        return jac(fld[:])
+    end
+    return jacobian
+end
+
+

src/types.jl

@@ -1,7 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-export IntegrationPoint, Assembly
+abstract Equation
+
 
 """
 Integration point