added tests

src/JuliaFEM.jl

@@ -6,10 +6,10 @@ This is JuliaFEM -- Finite Element Package
 """
 module JuliaFEM
 
-using Lexicon
 using Logging
 @Logging.configure(level=DEBUG)
 
+using Lexicon
 using ForwardDiff
 autodiffcache = ForwardDiffCache()
 
@@ -28,7 +28,6 @@ function Base.linspace(X1, X2, n)
     [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
 end
 
-
 include("types.jl")  # type definitions
 include("interpolate.jl")  # interpolation routines
 
@@ -36,21 +35,26 @@ include("interpolate.jl")  # interpolation routines
 include("elements.jl")
 include("lagrange.jl") # Lagrange elements
 #include("hierarchical.jl") # P-elements
-include("integrate.jl")  # integration points
 
+### EQUATIONS ###
+include("integrate.jl")  # default integration points for elements
 include("equations.jl")
 include("problems.jl")
-include("solvers.jl")
-
-# pre- and postprocess
-include("xdmf.jl")
-include("abaqus_reader.jl")
-#include("interfaces.jl")
 
+### FORMULATIION ###
 include("dirichlet.jl")
 include("heat.jl")
-#include("elasticity_solver.jl")
+include("elasticity.jl")
+
+### ASSEMBLY + SOLVE ###
+include("assembly.jl")
+include("solvers.jl")
+
+# PRE AND POSTPROCESS
+include("xdmf.jl")
+include("abaqus_reader.jl")
 
 end # module
 
 FEM = JuliaFEM
+

src/assembly.jl

@@ -0,0 +1,63 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Functions to handle global assembly of problem
+
+""" Global assembly. """
+type GlobalAssembly <: Assembly
+    ndofs :: Int
+    mass_matrix :: SparseMatrixCSC
+    stiffness_matrix :: SparseMatrixCSC
+    force_vector :: SparseMatrixCSC
+end
+
+""" Initialize global assembly of size ndofs. """
+function initialize_global_assembly(ndofs::Int=1)
+    mass_matrix = spzeros(ndofs, ndofs)
+    stiffness_matrix = spzeros(ndofs, ndofs)
+    force_vector = spzeros(ndofs, 1)
+    return GlobalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector)
+end
+
+""" Initialize global assembly, get dimension from problem. """
+function initialize_global_assembly(problem::Problem)
+    dim, ndofs = size(problem)
+    return initialize_global_assembly(ndofs)
+end
+
+""" Initialize or empty workspace for global assembly. """
+function initialize_global_assembly!(assembly::GlobalAssembly, problem::Problem)
+    ndofs = prod(size(problem))
+    if ndofs != assembly.ndofs
+        # if problem size changes, automatically initialize new work space
+        assembly.ndofs = ndofs
+        assembly.mass_matrix = spzeros(ndofs, ndofs)
+        assembly.stiffness_matrix = spzeros(ndofs, ndofs)
+        assembly.force_vector = spzeros(ndofs, 1)
+        return
+    end
+    # otherwise, empty workspace ready for next iteration
+    fill!(assembly.mass_matrix, 0.0)
+    fill!(assembly.stiffness_matrix, 0.0)
+    fill!(assembly.force_vector, 0.0)
+    return
+end
+
+""" Calculate global assembly for a problem. """
+function calculate_global_assembly!(assembly::GlobalAssembly, problem::Problem, time::Number=Inf)
+
+    unknown_field_name = get_unknown_field_name(problem)
+    initialize_global_assembly!(assembly, problem) # zero all
+    dim, ndofs = size(problem)
+    Logging.info("assembling problem for $unknown_field_name")
+    Logging.info("dimension of unknown field: $dim, problem dofs: $ndofs")
+    local_assembly = initialize_local_assembly()
+    for (i, equation) in enumerate(get_equations(problem))
+        calculate_local_assembly!(local_assembly, equation, unknown_field_name, time)
+        conn = get_connectivity(get_element(equation))
+        gdofs = get_gdofs(problem, equation)
+        assembly.mass_matrix[gdofs, gdofs] += local_assembly.mass_matrix
+        assembly.stiffness_matrix[gdofs, gdofs] += local_assembly.stiffness_matrix
+        assembly.force_vector[gdofs] += local_assembly.force_vector
+    end
+end

src/dirichlet.jl

@@ -34,7 +34,9 @@ function DBC2D2(element::Seg2)
     integration_points = [
         IntegrationPoint([-sqrt(1/3)], 1.0),
         IntegrationPoint([+sqrt(1/3)], 1.0)]
-    push!(element, FieldSet("reaction force"))
+    if !haskey(element, "reaction force")
+        element["reaction force"] = FieldSet()
+    end
     DBC2D2(element, integration_points)
 end
 Base.size(equation::DBC2D2) = (1, 2)

src/elasticity.jl

@@ -0,0 +1,130 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Elasticity problems
+
+abstract ElasticityProblem <: Problem
+abstract ElasticityEquation <: Equation
+
+### Formulation ###
+
+""" Calculate internal energy for elasticity equation.
+
+Override this to define your own material model. By default we use
+Saint Venant-Kirchhoff material model, which is simply
+
+    S(E) = λtr(E) + 2μE
+"""
+function get_internal_energy(equation::Equation, ip::IntegrationPoint,
+                             time::Number, F::Matrix)
+    element = get_element(equation)
+    basis = get_basis(element)
+    dbasis = grad(basis)
+
+    # material parameters
+    young = basis("youngs modulus", ip, time)
+    poisson = basis("poissons ratio", ip, time)
+    mu = young/(2*(1+poisson))
+    lambda = young*poisson/((1+poisson)*(1-2*poisson))
+    if isa(equation, PlaneStressElasticityEquation)
+        lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d
+    end
+
+    # material model
+    E = 1/2*(F'*F - I)  # strain
+    S = lambda*trace(E)*I + 2*mu*E
+    P = F*S
+
+    return P*dbasis(ip, time)
+end
+
+""" Elasticity equations.
+
+Formulation
+-----------
+
+Field equation is:
+∂u/∂t = ∇⋅f - b
+
+Weak form is: find u∈U such that ∀v in V
+
+    δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0
+
+where
+
+    ρ₀ = density
+    b₀ = displacement load
+    t₀ = displacement traction
+
+References
+----------
+
+https://en.wikipedia.org/wiki/Linear_elasticity
+https://en.wikipedia.org/wiki/Finite_strain_theory
+https://en.wikipedia.org/wiki/Stress_measures
+https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid
+https://en.wikipedia.org/wiki/Strain_energy_density_function
+https://en.wikipedia.org/wiki/Plane_stress
+https://en.wikipedia.org/wiki/Hooke's_law
+
+"""
+function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint,
+                             time::Number; variation=nothing)
+
+    element = get_element(equation)
+    basis = get_basis(element)
+    dbasis = grad(basis)
+
+    u = basis("displacement", ip, time, variation)
+    gradu = dbasis("displacement", ip, time, variation)
+    F = I + gradu # deformation gradient
+
+    # residual vector - internal energy
+    r = get_internal_energy(equation, ip, time, F)
+
+    # external forces - volume load
+    if haskey(element, "displacement load")
+        b = basis("displacement load", ip, time)
+        r -= b*basis(ip, time)
+    end
+
+    return vec(r)
+end
+has_residual_vector(equation::ElasticityEquation) = true
+
+### Problem 1 - plane elasticity ###
+
+abstract PlaneElasticityProblem <: ElasticityProblem
+
+type PlaneStressElasticityProblem <: PlaneElasticityProblem
+    unknown_field_name :: ASCIIString
+    unknown_field_dimension :: Int
+    equations :: Array{ElasticityEquation, 1}
+    element_mapping :: Dict{DataType, DataType}
+end
+
+function PlaneStressElasticityProblem(equations=[])
+    element_mapping = Dict(
+        Quad4 => CPS4)
+    return PlaneStressElasticityProblem("displacement", 2, equations, element_mapping)
+end
+
+### Equations ###
+
+abstract PlaneElasticityEquation <: ElasticityEquation
+abstract PlaneStressElasticityEquation <: PlaneElasticityEquation
+
+type CPS4 <: PlaneStressElasticityEquation
+    element :: Quad4
+    integration_points :: Array{IntegrationPoint, 1}
+end
+function CPS4(element::Quad4)
+    integration_points = get_default_integration_points(element)
+    if !haskey(element, "displacement")
+        element["displacement"] = FieldSet()
+        push!(element["displacement"], Field(Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))
+    end
+    CPS4(element, integration_points)
+end
+Base.size(equation::CPS4) = (2, 4)
+

src/elasticity_solver.jl

@@ -1,333 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-module elasticity_solver
-
-using ForwardDiff
-
-using Logging
-@Logging.configure(level=INFO)
-
-# Below this line is internal functions related to solver. They can be used
-# directly if needed or using general interface combining data model and
-# solver.
-
-"""
-This is dummy function. Testing doctests and documentation.
-
-Parameters
-----------
-x : Array{Float64, 1}
-
-Returns
--------
-Array{float64, 1}
-  x + 1
-
-Notes
------
-This is dummy function
-
-Raises
-------
-Exception
-  if things are not going right
-
-Examples
---------
->>> a = [1.0, 2.0, 3.0]
->>> dummy(a)
-[2.0, 3.0, 4.0]
-"""
-function dummy(a)
-    # not doing anything useful.
-    return a+1
-end
-
-
-"""
-Calculate local tangent stiffness matrix and residual force vector
-R = T - F for elasticity problem.
-
-Parameters
-----------
-X : Element coordinates
-u : Displacement field
-R : Residual force vector
-K : Tangent stiffness matrix
-basis : Basis functions
-dbasis : Derivative of basis functions
-lambda : Material parameter
-mu : Material parameter
-ipoints : integration points
-iweights : integration weights
-
-Returns
--------
-None
-
-Notes
------
-If material parameters are given in list, they are interpolated to gauss
-points using shape functions.
-"""
-function calc_local_matrices!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights)
-    dim, nnodes = size(X)
-    I = eye(dim)
-    R[:,:] = 0.0
-
-    #dF = zeros(dim, dim)
-
-    function calc_R!(u, R)
-        for m = 1:length(iweights)
-            w = iweights[m]
-            xi = ipoints[m, :]
-            # calculate material parameters
-            lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi))
-            mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi))
-            Jt = X*dbasis(xi)
-            detJ = det(Jt)
-            dbasisdX = dbasis(xi)*inv(Jt)
-
-            gradu = u*dbasisdX
-            F = I + gradu  # Deformation gradient
-            E = 1/2*(gradu' + gradu + gradu'*gradu)  # Green-Lagrange strain tensor
-            S = lambda*trace(E)*I + 2*mu*E  # PK2 stress tensor
-            P = F*S  # PK1 stress tensor
-
-            R[:,:] += w*P*dbasisdX'*detJ
-        end
-    end
-
-    # herlper for tangent stiffness matrix
-    function R!(u, R)
-        R[:] = 0
-        calc_R!(reshape(u, dim, nnodes), reshape(R, dim, nnodes))
-        #calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4))
-    end
-    Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
-
-    K[:, :] = Jacobian(reshape(u, dim*nnodes))
-    R!(reshape(u, dim*nnodes), reshape(R, dim*nnodes))
-
-end
-
-
-"""
-Assemble global stiffness matrix to I,J,V ready for sparse format
-
-Parameters
-----------
-ke : local matrix
-eldofs_ : Array
-  degrees of freedom
-I,J,V : Arrays for sparse matrix
-
-Notes
------
-eldofs can also be node ids for convenience. In that case dimension
-is calculated and eldofs are "extended" to problem dimension.
-"""
-function assemble!(ke, eldofs_, I, J, V)
-    n, m = size(ke)
-    dim = round(Int, n/length(eldofs_))
-    @debug("problem dim = ", dim)
-    if dim == 1
-        eldofs = eldofs_
-    else
-        eldofs = Int64[]
-        for i in eldofs_
-            for d in 1:dim
-                push!(eldofs, dim*(i-1)+d)
-            end
-        end
-        @debug("old eldofs", eldofs_)
-        @debug("new eldofs", eldofs)
-    end
-    for i in 1:n
-        for j in 1:m
-            push!(I, eldofs[i])
-            push!(J, eldofs[j])
-            push!(V, ke[i,j])
-        end
-    end
-end
-
-
-"""
-Assemble global RHS to I,V ready for sparse format
-
-Parameters
-----------
-fe : local vector
-eldofs_ : Array
-  degrees of freedom
-I,V : Arrays for sparse matrix
-
-Notes
------
-eldofs can also be node ids for convenience. In that case dimension
-is calculated and eldofs are "extended" to problem dimension.
-"""
-function assemble!(fe, eldofs_, I, V)
-    n = length(fe)
-    dim = round(Int, n/length(eldofs_))
-    if dim == 1
-        eldofs = eldofs_
-    else
-        eldofs = Int64[]
-        for i in eldofs_
-            for d in 1:dim
-                push!(eldofs, dim*(i-1)+d)
-            end
-        end
-    end
-
-    for i in 1:n
-        push!(I, eldofs[i])
-        push!(V, fe[i])
-    end
-end
-
-
-"""
-Eliminate Dirichlet boundary conditions from matrix
-
-Parameters
-----------
-dirichletbc : array [dim x nnodes]
-I, J, V : sparse matrix arrays
-
-Returns
--------
-I, J, V : boundary conditions removed
-
-Notes
------
-pros:
-- matrix assembly remains positive definite
-cons:
-- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.)
-- if u != 0 in dirichlet boundary requires extra care
-
-Raises
-------
-Exception, if displacement boundary conditions given, i.e.
-DX=2 for some node, for example.
-
-"""
-function eliminate_boundary_conditions(dirichletbc, I, J, V)
-    if any(dirichletbc .> 0)
-        throw("displacement boundary condition not supported")
-    end
-    # dofs to remove
-    free_dofs = find(isnan(dirichletbc))
-    remove_dofs = find(!isnan(dirichletbc))
-    @debug("Removing dofs: ", remove_dofs)
-    # this can be done more clever by removing corresponging indexes from I, J, and V
-    A = sparse(I, J, V)
-    A = A[free_dofs, free_dofs]
-    return findnz(A)
-end
-
-"""
-Eliminate Dirichlet boundary conditions from vector
-
-Parameters
-----------
-dirichletbc : array [dim x nnodes]
-I, V : sparse vector arrays
-
-Returns
--------
-I, V : boundary conditions removed
-
-Notes
------
-pros:
-- matrix assembly remains positive definite
-cons:
-- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.)
-- if u != 0 in dirichlet boundary requires extra care
-
-Raises
-------
-Exception, if displacement boundary conditions given, i.e.
-DX=2 for some node, for example.
-""" 
-function eliminate_boundary_conditions(dirichletbc, I, V)
-    if any(dirichletbc .> 0)
-        throw("displacement boundary condition not supported")
-    end
-    # dofs to remove
-    free_dofs = find(isnan(dirichletbc))
-    remove_dofs = find(!isnan(dirichletbc))
-    @debug("Removing dofs: ", remove_dofs)
-    # this can be done more clever by removing corresponging indexes from I, J, and V
-    A = sparsevec(I, V)
-    A = A[free_dofs]
-    @debug("new vector: ", A)
-    i, j, v = findnz(A)
-    return i, v
-end
-
-
-
-"""
-Solve one increment of elasticity problem
-""" 
-function solve_elasticity_increment!(X, u, du, elmap, nodalloads,
-                                     dirichletbc, lambda, mu, N, dNdchi, ipoints,
-                                     iweights)
-    if length(size(elmap)) == 1
-        # quick hack for just one element
-        elmap = elmap''
-    end
-    nelnodes, nelements = size(elmap)
-    dim, nnodes = size(u)
-    dofs = dim*nelnodes
-
-    Imat = Int64[]
-    Jmat = Int64[]
-    Vmat = Float64[]
-    Ivec = Int64[]
-    Vvec = Float64[]
-
-    # FIXME: different number of nodes/element
-    R = zeros(dim, nelnodes)
-    Kt = zeros(dofs, dofs)
-
-    # this can be parallelized
-    for i in 1:nelements
-        eldofs = elmap[:,i]
-        calc_local_matrices!(X[:, eldofs], u[:, eldofs], R, Kt, N, dNdchi,
-                             lambda[eldofs], mu[eldofs], ipoints, iweights)
-        assemble!(Kt, eldofs, Imat, Jmat, Vmat)
-        assemble!(R, eldofs, Ivec, Vvec)
-    end
-
-    # add additional neumann boundary conditions to force vector
-    for (i, nodal_load) in enumerate(nodalloads)
-        if nodal_load == 0
-            continue
-        end
-        push!(Ivec, i)
-        push!(Vvec, -nodal_load)
-    end
-
-    # Create sparse matrix and vector
-    A = sparse(Imat, Jmat, Vmat)
-    b = sparsevec(Ivec, Vvec)
-
-    # Remove dirichlet boundary conditions
-    free_dofs = find(isnan(dirichletbc))
-    #Imat, Jmat, Vmat = eliminate_boundary_conditions(dirichletbc, Imat, Jmat, Vmat)
-    b = b[free_dofs]
-    A = A[free_dofs, free_dofs]
-
-    # solution
-    du[free_dofs] = lufact(A) \ -full(b)
-end
-
-
-end

src/elements.jl

@@ -1,83 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-#=
-Related notebooks
------------------
-
-2015-08-29-developing-juliafem.ipynb
-=#
-
 using FactCheck
 using ForwardDiff
 
 abstract Element
 
-""" Get FieldSet from element. """
-function Base.getindex(element::Element, field_name)
-    element.fields[field_name]
-end
-
-""" Add new FieldSet to element. """
-function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name)
-    fieldset.name = fieldset_name
-    element.fields[fieldset.name] = fieldset
-end
-function Base.push!(element::Element, fieldset::FieldSet)
-    element[fieldset.name] = fieldset
-end
-
-
-#= ELEMENT DEFINITIONS
-
-Example
--------
-
-This is example how to create new element. This is commented because I use code
-generation for simple elements like Lagrage elements. Feel free to use
-code generation but elements can be of course created manually too!
-
-abstract CG <: Element # create new element family "Continous Galerkin"
-
-type Quad4 <: CG
-    connectivity :: Array{Int, 1}
-    fields :: Dict{Any, Any}
-end
-
-""" Default contructor. """
-Quad4(connectivity) = Quad4(connectivity, Dict{Any, Any}())
-
-""" Return number of basis functions of this element. """
-get_number_of_basis_functions(el::Type{Quad4}) = 4
-
-""" Return element dimension (length of xi vector). """
-get_element_dimension(el::Type{Quad4}) = 2
-
-""" Return basis functions for this element (xi dim = 2, functions = 4). """
-function get_basis(el::Quad4, xi)
-    [(1-xi[1])*(1-xi[2])/4
-     (1+xi[1])*(1-xi[2])/4
-     (1+xi[1])*(1+xi[2])/4
-     (1-xi[1])*(1+xi[2])/4]
-end
-
-""" Return partial derivatives of basis functions. """
-function get_dbasisdxi(el::Quad4, xi)
-    [-(1-xi[2])/4.0    -(1-xi[1])/4.0
-      (1-xi[2])/4.0    -(1+xi[1])/4.0
-      (1+xi[2])/4.0     (1+xi[1])/4.0
-     -(1+xi[2])/4.0     (1-xi[1])/4.0]
-end
-
-End of example.
-
-=#
-
-# define size of your element as (dim, nbasis) tuple where first integer is spatial dimension and second is number of basis functions.
-# Base.size(element::Type{Element}) = nothing
-
-### COMMON ELEMENT ROUTINES ###
-
 """
 Test routine for element. If this passes, element interface is properly
 defined.
@@ -114,9 +42,9 @@ function test_element(element_type)
     end
 
     # try to interpolate some scalar field
-    push!(element, FieldSet("field1", [Field(0.0, collect(1:n))]))
+    element["field1"] = Field(0.0, collect(1:n))
     # TODO: how to parametrize this?
-    push!(element, FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))
+    element["geometry"] = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
 
     # evaluate basis functions at middle point of element
     basis = get_basis(element)
@@ -134,26 +62,94 @@ function test_element(element_type)
     Logging.info("Element $element_type passed tests.")
 end
 
+""" Get FieldSet from element. """
+function Base.getindex(element::Element, field_name)
+    element.fields[field_name]
+end
+
+"""Add new FieldSet to element.
+
+Examples
+--------
+>>> field = Field(0.0, [1, 2, 3, 4])
+>>> fieldset = FieldSet("geometry", Field[field])
+>>> element["geometry"] = fieldset
+JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
+"""
+function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name)
+    fieldset.name = fieldset_name
+    element.fields[fieldset.name] = fieldset
+end
+
+"""Add new FieldSet to element.
+
+Examples
+--------
+>>> field = Field(0.0, [1, 2, 3, 4])
+>>> element["geometry"] = field
+JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
+"""
+function Base.setindex!(element::Element, field::Field, fieldset_name)
+    element[fieldset_name] = FieldSet(field)
+end
+
+"""Add new FieldSet to element.
+
+Examples
+--------
+>>> element["geometry"] = [1, 2, 3, 4]
+JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
+"""
+function Base.setindex!(element::Element, field_data::Union{Number, Array}, fieldset_name)
+    element[fieldset_name] = Field(field_data)
+end
+
+"""Add new FieldSet to element.
+
+Notes
+-----
+This last version takes tuple and each cell in tuple is converted to new field.
+Time in field is 0.0, 1.0, ..., n
+
+Examples
+--------
+>>> element["load"] = (1, 2)
+JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("load"=>JuliaFEM.FieldSet("load",JuliaFEM.Field[JuliaFEM.Field{Int64}(0.0,0,1),JuliaFEM.Field{Int64}(1.0,0,2)])))
+"""
+function Base.setindex!(element::Element, field_data::Tuple, fieldset_name)
+    fields = Field[Field(Float64(i-1), field) for (i,field) in enumerate(field_data)]
+    element.fields[fieldset_name] = FieldSet(fieldset_name, fields)
+end
 
 function get_connectivity(el::Element)
     el.connectivity
 end
 
-type MixedFunctionSpace
+abstract AbstractFunctionSpace
+
+type FunctionSpace <: AbstractFunctionSpace
+    element :: Element
+end
+
+type GradientFunctionSpace <: AbstractFunctionSpace
+    element :: Element
+end
+
+type MixedFunctionSpace <: AbstractFunctionSpace
     element1 :: Element
     element2 :: Element
 end
 
-type FunctionSpace
-    element :: Element
+function get_basis(element::Element)
+    return FunctionSpace(element)
 end
 
-type GradientFunctionSpace
-    element :: Element
+function get_dbasis(element::Element)
+    return GradientFunctionSpace(element)
 end
 
 function grad(u::FunctionSpace)
-    GradientFunctionSpace(u.element)
+    return GradientFunctionSpace(u.element)
 end
 
 """ Evaluate field on element function space. """
@@ -195,8 +191,8 @@ end
 call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
 call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
 # i think these will be the most called functions.
-call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation)
-call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation)
 call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
 call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
 
@@ -215,6 +211,7 @@ function get_fieldset(u::FunctionSpace, field_name)
     return u.element[field_name]
 end
 
+""" Get a determinant of element in point ξ. """
 function LinAlg.det(u::FunctionSpace, xi::Vector, t::Number=Inf)
     X = u.element["geometry"](t)
     dN = u.element.basis.dbasisdxi(xi)
@@ -229,156 +226,13 @@ function LinAlg.det(u::FunctionSpace)
     return (args...) -> det(u, args...)
 end
 
-function get_basis(element::Element)
-    return FunctionSpace(element)
-end
-
 Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
 Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
 Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)
 Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...)
 
-
 """ Check does fieldset exist. """
 function Base.haskey(element::Element, what)
     haskey(element.fields, what)
 end
 
-
-# FIXME: These two needs integration -- maybe not in elements.jl ..?
-"""
-Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min!
-
-Parameters
-----------
-f::Function
-    Needs to take (el::Element, xi::Vector) as argument
-fixed_coeffs::Int[]
-    These coefficients are not changed during fitting -> constrained optimizatio
-"""
-function fit_field!(el::Element, field, f, fixed_coeffs=Int[])
-    w = [
-        128/225,
-        (332+13*sqrt(70))/900,
-        (332+13*sqrt(70))/900,
-        (332-13*sqrt(70))/900,
-        (332-13*sqrt(70))/900]
-    xi = Vector[
-        [0.0],
-        [ 1/3*sqrt(5 - 2*sqrt(10/7))],
-        [-1/3*sqrt(5 - 2*sqrt(10/7))],
-        [ 1/3*sqrt(5 + 2*sqrt(10/7))],
-        [-1/3*sqrt(5 + 2*sqrt(10/7))]]
-    n = get_number_of_basis_functions(el)
-    fld = get_field(el, field)
-    nfld = length(fld[1])
-    #Logging.debug("dim of field $field: $nfld")
-
-    M = zeros(n, n)
-    b = zeros(n, nfld)
-    for i=1:length(w)
-        detJ = get_detJ(el, xi[i])
-        N = get_basis(el, xi[i])
-        M += w[i]*N*N'*detJ
-        fi = f(el, xi[i])
-        for j=1:nfld
-            b[:, j] += w[i]*N*fi[j]*detJ
-        end
-    end
-
-    coeffs = zeros(n)
-    for j=1:nfld
-        for k=1:n
-            coeffs[k] = fld[k][j]
-        end
-        if length(fixed_coeffs) != 0
-            # constrained problem, some coefficients are fixed
-            N = Int[] # rest of coeffs
-            S = Int[] # fixed coeffs
-            for i = 1:n
-                if i in fixed_coeffs
-                    push!(S, i)
-                else
-                    push!(N, i)
-                end
-            end
-            lhs = M[N,N]
-            rhs = b[N,j] - M[N,S]*coeffs[S]
-            coeffs[N] = lhs \ rhs
-        else
-            coeffs[:] = M \ b[:,j]
-        end
-        for k=1:n
-            fld[k][j] = coeffs[k]
-        end
-    end
-    set_field(el, field, fld)
-    return
-end
-
-
-"""
-Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min!
-"""
-function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[])
-    w = [
-        128/225,
-        (332+13*sqrt(70))/900,
-        (332+13*sqrt(70))/900,
-        (332-13*sqrt(70))/900,
-        (332-13*sqrt(70))/900]
-    xi = Vector[
-        [0.0],
-        [ 1/3*sqrt(5 - 2*sqrt(10/7))],
-        [-1/3*sqrt(5 - 2*sqrt(10/7))],
-        [ 1/3*sqrt(5 + 2*sqrt(10/7))],
-        [-1/3*sqrt(5 + 2*sqrt(10/7))]]
-    n = get_number_of_basis_functions(el)
-    fld = get_field(el, field)
-    nfld = length(fld[1])
-    #Logging.debug("dim of field $field: $nfld")
-
-    M = zeros(n, n)
-    b = zeros(n, nfld)
-    for i=1:length(w)
-        detJ = get_detJ(el, xi[i])
-        dNdxi = get_dbasisdxi(el, xi[i])
-        dNdX = dNdxi / detJ
-        M += w[i]*dNdX*dNdX'*detJ
-        fi = f(el, xi[i])
-        for j=1:nfld
-            b[:, j] += w[i]*dNdX*fi[j]*detJ
-        end
-    end
-
-    coeffs = zeros(n)
-    for j=1:nfld
-        for k=1:n
-            coeffs[k] = fld[k][j]
-        end
-        if length(fixed_coeffs) != 0
-            #Logging.info("constrained problem, some coefficients are fixed")
-            N = Int[] # rest of coeffs
-            S = Int[] # fixed coeffs
-            for i = 1:n
-                if i in fixed_coeffs
-                    push!(S, i)
-                else
-                    push!(N, i)
-                end
-            end
-            lhs = M[N,N]
-            rhs = b[N,j] - M[N,S]*coeffs[S]
-            coeffs[N] = lhs \ rhs
-        else
-            coeffs[:] = M \ b[:,j]
-        end
-        for k=1:n
-            fld[k][j] = coeffs[k]
-        end
-    end
-    set_field(el, field, fld)
-    return
-end
-
-

src/equations.jl

@@ -1,6 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+# Functions to handle element level things -- integration, assembly, ...
+
 abstract Equation
 
 abstract Assembly
@@ -55,17 +57,61 @@ function initialize_local_assembly!(assembly::LocalAssembly, equation::Equation)
 end
 
 has_mass_matrix(equation::Equation) = false
-get_mass_matrix(equation::Equation, ip, time) = nothing
+function get_mass_matrix(equation::Equation, ip, time=Inf, problem=nothing)
+    get_mass_matrix(equation, ip, time)
+end
+function get_mass_matrix(equation::Equation, ip, time=Inf)
+    get_mass_matrix(equation, ip)
+end
+function get_mass_matrix(equation::Equation, ip)
+    nothing
+end
+
 has_stiffness_matrix(equation::Equation) = false
-get_stiffness_matrix(equation::Equation, ip, time) = nothing
+function get_stiffness_matrix(equation::Equation, ip, time=Inf, problem=nothing)
+    get_stiffness_matrix(equation, ip, time)
+end
+function get_stiffness_matrix(equation::Equation, ip, time=Inf)
+    get_stiffness_matrix(equation, ip)
+end
+function get_stiffness_matrix(equation::Equation, ip)
+    nothing
+end
+
 has_force_vector(equation::Equation) = false
-get_force_vector(equation::Equation, ip, time) = nothing
+function get_force_vector(equation::Equation, ip, time=Inf, problem=nothing)
+    get_force_vector(equation, ip, time)
+end
+function get_force_vector(equation::Equation, ip, time=Inf)
+    get_force_vector(equation, ip)
+end
+function get_force_vector(equation::Equation, ip)
+    nothing
+end
+
 has_residual_vector(equation::Equation) = false
-get_residual_vector(equation::Equation, ip, time) = nothing
+function get_residual_vector(equation::Equation, ip, time=Inf, problem=nothing)
+    get_residual_vector(equation, ip, time)
+end
+function get_residual_vector(equation::Equation, ip, time=Inf)
+    get_residual_vector(equation, ip)
+end
+function get_residual_vector(equation::Equation, ip)
+    nothing
+end
+
 has_potential_energy(equation::Equation) = false
-get_potential_energy(equation::Equation, ip, time) = nothing
+function get_potential_energy(equation::Equation, ip, time=Inf, problem=nothing)
+    get_potential_energy(equation, ip, time)
+end
+function get_potential_energy(equation::Equation, ip, time=Inf)
+    get_potential_energy(equation, ip)
+end
+function get_potential_energy(equation::Equation, ip)
+    nothing
+end
+
 get_element(equation::Equation) = equation.element
-get_number_of_dofs(equation::Equation) = nothing
 get_integration_points(equation::Equation) = equation.integration_points
 
 
@@ -85,19 +131,19 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
         for ip in get_integration_points(equation)
             s = ip.weight*detJ(ip)
             if has_mass_matrix(equation)
-                assembly.mass_matrix += s*get_mass_matrix(equation, ip, time)
+                assembly.mass_matrix += s*get_mass_matrix(equation, ip, time, problem)
             end
             if has_stiffness_matrix(equation)
-                assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time)
+                assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time, problem)
             end
             if has_force_vector(equation)
-                assembly.force_vector += s*get_force_vector(equation, ip, time)[:]
-            end
-            # external loads -- if any nodal loads is defined add to force vector
-            if haskey(element, "$unknown_field_name nodal load")
-                assembly.force_vector += element["$unknown_field_name nodal load"](time)[:]
+                assembly.force_vector += s*get_force_vector(equation, ip, time, problem)
             end
         end
+        # external loads -- if any nodal loads is defined add to force vector
+        if haskey(element, "$unknown_field_name nodal load")
+            assembly.force_vector += vec(element["$unknown_field_name nodal load"](time))
+        end
     end
 
     # 2. variational / energy form - user has defined some potential energy / variational form
@@ -116,14 +162,14 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
             # external energy -- if any nodal loads is defined, decrease from potential energy
             if haskey(element, "$unknown_field_name nodal load")
                 P = element["$unknown_field_name nodal load"](time)
-                assembly.potential_energy -= dot(P[:], df[:])
+                assembly.potential_energy -= dot(vec(P), vec(df))
             end
             if isa(assembly.potential_energy, Array)
                 return assembly.potential_energy[1]
             end
             return assembly.potential_energy
         end
-        hessian, allresults = ForwardDiff.hessian(potential_energy, field[:],
+        hessian, allresults = ForwardDiff.hessian(potential_energy, vec(field),
                                                   AllResults, cache=autodiffcache)
         assembly.stiffness_matrix += hessian
         assembly.force_vector -= ForwardDiff.gradient(allresults)  # <--- minus explained in tutorial
@@ -144,15 +190,14 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
             end
             # external loads -- if any nodal loads is defined, remove from residual
             if haskey(element, "$unknown_field_name nodal load")
-                assembly.residual_vector -= element["$unknown_field_name nodal load"](time)[:]
+                assembly.residual_vector -= vec(element["$unknown_field_name nodal load"](time))
             end
             return assembly.residual_vector
         end
-        jacobian, allresults = ForwardDiff.jacobian(residual_vector, field[:],
+        jacobian, allresults = ForwardDiff.jacobian(residual_vector, vec(field),
                                                     AllResults, cache=autodiffcache)
         assembly.stiffness_matrix += jacobian
         assembly.force_vector -= ForwardDiff.value(allresults)  # <-- minus explained in tutorial
     end
 
 end
-

src/field_utils.jl

@@ -0,0 +1,140 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# There are here for now until I figure a better place for them.
+
+
+"""
+Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min!
+
+Parameters
+----------
+f::Function
+    Needs to take (el::Element, xi::Vector) as argument
+fixed_coeffs::Int[]
+    These coefficients are not changed during fitting -> constrained optimizatio
+"""
+function fit_field!(el::Element, field, f, fixed_coeffs=Int[])
+    w = [
+        128/225,
+        (332+13*sqrt(70))/900,
+        (332+13*sqrt(70))/900,
+        (332-13*sqrt(70))/900,
+        (332-13*sqrt(70))/900]
+    xi = Vector[
+        [0.0],
+        [ 1/3*sqrt(5 - 2*sqrt(10/7))],
+        [-1/3*sqrt(5 - 2*sqrt(10/7))],
+        [ 1/3*sqrt(5 + 2*sqrt(10/7))],
+        [-1/3*sqrt(5 + 2*sqrt(10/7))]]
+    n = get_number_of_basis_functions(el)
+    fld = get_field(el, field)
+    nfld = length(fld[1])
+    #Logging.debug("dim of field $field: $nfld")
+
+    M = zeros(n, n)
+    b = zeros(n, nfld)
+    for i=1:length(w)
+        detJ = get_detJ(el, xi[i])
+        N = get_basis(el, xi[i])
+        M += w[i]*N*N'*detJ
+        fi = f(el, xi[i])
+        for j=1:nfld
+            b[:, j] += w[i]*N*fi[j]*detJ
+        end
+    end
+
+    coeffs = zeros(n)
+    for j=1:nfld
+        for k=1:n
+            coeffs[k] = fld[k][j]
+        end
+        if length(fixed_coeffs) != 0
+            # constrained problem, some coefficients are fixed
+            N = Int[] # rest of coeffs
+            S = Int[] # fixed coeffs
+            for i = 1:n
+                if i in fixed_coeffs
+                    push!(S, i)
+                else
+                    push!(N, i)
+                end
+            end
+            lhs = M[N,N]
+            rhs = b[N,j] - M[N,S]*coeffs[S]
+            coeffs[N] = lhs \ rhs
+        else
+            coeffs[:] = M \ b[:,j]
+        end
+        for k=1:n
+            fld[k][j] = coeffs[k]
+        end
+    end
+    set_field(el, field, fld)
+    return
+end
+
+
+"""
+Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min!
+"""
+function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[])
+    w = [
+        128/225,
+        (332+13*sqrt(70))/900,
+        (332+13*sqrt(70))/900,
+        (332-13*sqrt(70))/900,
+        (332-13*sqrt(70))/900]
+    xi = Vector[
+        [0.0],
+        [ 1/3*sqrt(5 - 2*sqrt(10/7))],
+        [-1/3*sqrt(5 - 2*sqrt(10/7))],
+        [ 1/3*sqrt(5 + 2*sqrt(10/7))],
+        [-1/3*sqrt(5 + 2*sqrt(10/7))]]
+    n = get_number_of_basis_functions(el)
+    fld = get_field(el, field)
+    nfld = length(fld[1])
+    #Logging.debug("dim of field $field: $nfld")
+
+    M = zeros(n, n)
+    b = zeros(n, nfld)
+    for i=1:length(w)
+        detJ = get_detJ(el, xi[i])
+        dNdxi = get_dbasisdxi(el, xi[i])
+        dNdX = dNdxi / detJ
+        M += w[i]*dNdX*dNdX'*detJ
+        fi = f(el, xi[i])
+        for j=1:nfld
+            b[:, j] += w[i]*dNdX*fi[j]*detJ
+        end
+    end
+
+    coeffs = zeros(n)
+    for j=1:nfld
+        for k=1:n
+            coeffs[k] = fld[k][j]
+        end
+        if length(fixed_coeffs) != 0
+            #Logging.info("constrained problem, some coefficients are fixed")
+            N = Int[] # rest of coeffs
+            S = Int[] # fixed coeffs
+            for i = 1:n
+                if i in fixed_coeffs
+                    push!(S, i)
+                else
+                    push!(N, i)
+                end
+            end
+            lhs = M[N,N]
+            rhs = b[N,j] - M[N,S]*coeffs[S]
+            coeffs[N] = lhs \ rhs
+        else
+            coeffs[:] = M \ b[:,j]
+        end
+        for k=1:n
+            fld[k][j] = coeffs[k]
+        end
+    end
+    set_field(el, field, fld)
+    return
+end

src/heat.jl

@@ -6,7 +6,66 @@
 abstract HeatProblem <: Problem
 abstract HeatEquation <: Equation
 
-### Plane heat problem + equations ###
+### Formulation ###
+
+""" Heat equations.
+
+Formulation
+-----------
+
+Field equation is:
+
+    ρc∂u/∂t = ∇⋅(k∇u) + f
+
+Weak form is: find u∈U such that ∀v in V
+
+    ∫k∇u∇v dx = ∫fv dx + ∫gv ds,
+
+where
+
+    k = temperature thermal conductivity    defined on volume
+    f = temperature load                    defined on volume
+    g = temperature flux                    defined on boundary
+
+References
+----------
+https://en.wikipedia.org/wiki/Heat_equation
+
+"""
+function calculate_local_assembly!(assembly::LocalAssembly, equation::HeatEquation,
+                                            unknown_field_name::ASCIIString, time::Number=Inf,
+                                            problem=nothing)
+
+    initialize_local_assembly!(assembly, equation)
+
+    element = get_element(equation)
+    basis = get_basis(element)
+    dbasis = grad(basis)
+    detJ = det(basis)
+    for ip in get_integration_points(equation)
+        w = ip.weight * detJ(ip)
+        N = basis(ip, time)
+        if haskey(element, "density")
+            ρ = basis("density", ip, time)
+            assembly.mass_matrix += w * ρ*N'*N
+        end
+        if haskey(element, "temperature thermal conductivity")
+            dN = dbasis(ip, time)
+            k = basis("temperature thermal conductivity", ip, time)
+            assembly.stiffness_matrix += w * k*dN'*dN
+        end
+        if haskey(element, "temperature load")
+            f = basis("temperature load", ip, time)
+            assembly.force_vector += w * N'*f
+        end
+        if haskey(element, "temperature flux")
+            g = basis("temperature flux", ip, time)
+            assembly.force_vector += w * N'*g
+        end
+    end
+end
+
+### Problems ###
 
 type PlaneHeatProblem <: HeatProblem
     unknown_field_name :: ASCIIString
@@ -23,7 +82,7 @@ function PlaneHeatProblem()
     return PlaneHeatProblem("temperature", 1, [], element_mapping)
 end
 
-
+### Equations ###
 
 """ Diffusive heat transfer for 4-node bilinear element. """
 type DC2D4 <: HeatEquation
@@ -32,7 +91,9 @@ type DC2D4 <: HeatEquation
 end
 function DC2D4(element::Quad4)
     integration_points = get_default_integration_points(element)
-    push!(element, FieldSet("temperature"))
+    if !haskey(element, "temperature")
+        element["temperature"] = FieldSet()
+    end
     DC2D4(element, integration_points)
 end
 Base.size(equation::DC2D4) = (1, 4)
@@ -44,34 +105,10 @@ type DC2D2 <: HeatEquation
 end
 function DC2D2(element::Seg2)
     integration_points = get_default_integration_points(element)
-    push!(element, FieldSet("temperature"))
+    if !haskey(element, "temperature")
+        element["temperature"] = FieldSet()
+    end
     DC2D2(element, integration_points)
 end
 Base.size(equation::DC2D2) = (1, 2)
 
-function calculate_local_assembly!(assembly::LocalAssembly, equation::HeatEquation,
-                                            unknown_field_name::ASCIIString, time::Number=Inf,
-                                            problem=nothing)
-
-    initialize_local_assembly!(assembly, equation)
-
-    element = get_element(equation)
-    basis = get_basis(element)
-    dbasis = grad(basis)
-    detJ = det(basis)
-    for ip in get_integration_points(equation)
-        w = ip.weight * detJ(ip)
-        # evaluate fields in integration point
-        ρ = basis("density", ip, time)
-        k = basis("temperature thermal conductivity", ip, time)
-        f = basis("temperature load", ip, time)
-        # evaluate basis functions and gradient in integration point
-        N = basis(ip, time)
-        dN = dbasis(ip, time)
-        # do assembly
-        assembly.mass_matrix += w * ρ*N'*N
-        assembly.stiffness_matrix += w * k*dN'*dN
-        assembly.force_vector += w * N'*f
-    end
-end
-

src/lagrange.jl

@@ -15,7 +15,6 @@ function calculate_lagrange_basis(P, X)
     for i=1:nbasis
         A[i,:] = P(X[:, i])
     end
-#   Logging.debug("Calculating inverse of A")
     invA = inv(A)'
     basis(xi) = (invA*P(xi))'
     dbasisdxi(xi) = (ForwardDiff.jacobian((xi) -> invA*P(xi), xi, cache=autodiffcache))'
@@ -34,9 +33,6 @@ macro create_lagrange_element(element_name, element_description, X, P)
     eltype = esc(element_name)
     quote
         global get_element_description
-        #global get_number_of_basis_functions, get_element_dimension
-        #dim = size($X, 1)
-        #nbasis = size($X, 2)
         basis, dbasisdxi = calculate_lagrange_basis($P, $X)
         type $eltype <: CG
             connectivity :: Array{Int, 1}
@@ -47,8 +43,6 @@ macro create_lagrange_element(element_name, element_description, X, P)
             $eltype(connectivity, Basis(basis, dbasisdxi), Dict())
         end
         get_element_description(el::Type{$eltype}) = $element_description
-        #get_number_of_basis_functions(el::Type{$eltype}) = nbasis
-        #get_element_dimension(el::Type{$eltype}) = dim
         Base.size(el::Type{$eltype}) = Base.size($X)
     end
 end

src/problems.jl

@@ -5,22 +5,72 @@ abstract Problem
 abstract BoundaryProblem <: Problem
 abstract FieldProblem <: Problem
 
+""" Return all equations beloging to this problem. """
 function get_equations(problem::Problem)
     problem.equations
 end
 
+""" Return the dimension of the unknown field of this problem. """
 function get_unknown_field_dimension(problem::Problem)
     problem.unknown_field_dimension
 end
 
+""" Return the name of the unknown field of this problem. """
 function get_unknown_field_name(problem::Problem)
     problem.unknown_field_name
 end
 
-""" Add new element to problem. """
+"""
+Add new equation to problem.
+
+Parameters
+----------
+problem
+element
+
+Notes
+-----
+Equation is automatically created during process based on problem
+element -> equation mapping and element type.
+"""
 function Base.push!(problem::Problem, element::Element)
     element_type = typeof(element)
     equation_type = problem.element_mapping[element_type]
     push!(problem.equations, equation_type(element))
 end
 
+"""
+Return the size of the problem, i.e.
+maximum number of connectivity × unknown field dimension.
+"""
+function Base.size(problem::Problem)
+    mc = 0
+    for equation in get_equations(problem)
+        element = get_element(equation)
+        mc = max(mc, get_connectivity(element)...)
+    end
+    dim = get_unknown_field_dimension(problem)
+    return (dim, dim*mc)
+end
+
+""" Assign new equation mapping to problem for some element.
+
+Examples
+--------
+>>> p = PlaneHeatProblem()
+>>> p[Seg2] = DC2D2
+
+"""
+function Base.setindex!(problem::Problem, equation, element)
+    problem.element_mapping[element] = equation
+end
+
+""" Return global degrees of freedom of element in matrix level. """
+function get_gdofs(problem::Problem, equation::Equation)
+    dim = get_unknown_field_dimension(problem)
+    element = get_element(equation)
+    conn = get_connectivity(element)
+    gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
+    return gdofs
+end
+

src/shape_functions.jl

@@ -1,3 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-

src/solvers.jl

@@ -5,17 +5,89 @@
 
 abstract Solver
 
+"""
+Solve field equations for single element with some dofs fixed. This can be used
+to test nonlinear element formulations.
+"""
+function solve!(equation::Equation, unknown_field_name::ASCIIString,
+                 free_dofs::Array{Int, 1}, time::Number=Inf; 
+                 max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+    element = get_element(equation)
+    x0 = element[unknown_field_name](-Inf)
+    x = zeros(prod(size(equation)))
+    dx = fill!(similar(x), 0.0)
+    la = initialize_local_assembly()
+    for i=1:max_iterations
+        calculate_local_assembly!(la, equation, unknown_field_name)
+        A = la.stiffness_matrix[free_dofs, free_dofs]
+        b = la.force_vector[free_dofs]
+        if dump_matrices
+            dump(full(A))
+            dump(full(b)')
+        end
+        dx[free_dofs] = A \ b
+        x += dx
+        new_field = similar(x0, x)
+        new_field.time = time
+        new_field.increment = i
+        push!(element[unknown_field_name], new_field)
+        if norm(dx) < tolerance
+            return
+        end
+    end
+    Logging.err("Did not converge in $max_iterations iterations")
+end
+
+"""
+Solve field equations for a single problem with some dofs fixed. This can be used
+to test nonlinear element formulations. Dirichlet boundary is assumed to be homogeneous
+and degrees of freedom are eliminated. So if boundary condition is known in nodal
+points and everything is zero this should be quite good.
+"""
+function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=Inf;
+                max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+    ga = initialize_global_assembly(problem)
+    x = zeros(ga.ndofs)
+    dx = fill!(similar(x), 0.0)
+    field_name = get_unknown_field_name(problem)
+    dim = get_unknown_field_dimension(problem)
+    for i=1:max_iterations
+        calculate_global_assembly!(ga, problem)
+        A = ga.stiffness_matrix[free_dofs, free_dofs]
+        b = ga.force_vector[free_dofs]
+        if dump_matrices
+            dump(full(A))
+            dump(full(b)')
+        end
+        dx[free_dofs] = lufact(A) \ full(b)
+        x += dx
+        for equation in get_equations(problem)
+            element = get_element(equation)
+            conn = get_connectivity(element)
+            gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
+            old_field = element[field_name](Inf)
+            new_field = similar(old_field, full(x[gdofs]))
+            new_field.time = time
+            new_field.increment = i
+            push!(element[field_name], new_field)
+        end
+        if norm(dx) < tolerance
+            return
+        end
+    end
+    Logging.err("Did not converge in $max_iterations iterations")
+end
+
 """ Add new problem to solver. """
 function add_problem!(solver::Solver, problem::Problem)
     push!(solver.problems, problem)
 end
+
 function Base.push!(solver::Solver, problem::Problem)
     push!(solver.problems, problem)
 end
 
-"""
-Get all problems assigned to solver
-"""
+""" Get all problems assigned to solver. """
 function get_problems(s::Solver)
     return s.problems
 end
@@ -25,6 +97,7 @@ end
 type SimpleSolver <: Solver
     problems
 end
+
 """ Default initializer. """
 function SimpleSolver()
     SimpleSolver(Problem[])
@@ -33,34 +106,32 @@ end
 """
 Call solver to solve a set of problems.
 
-This is simple serial solver for demonstration purposes. It handles the most
+This is a simple direct solver for demonstration purposes. It handles the
 common situation, i.e., some main field problem and it's Dirichlet boundary.
+
+    Au + C'λ = f
+    Cu       = g
+
 """
-function call(solver::SimpleSolver, t)
-    problems = get_problems(solver)
-    problem1 = problems[1]
-    problem2 = problems[2]
+function call(solver::SimpleSolver, time::Number=Inf)
+    p1, p2 = get_problems(solver)
 
-    # calculate order of degrees of freedom in global matrix
-    # and set the ordering to problems
-    dofmap = calculate_global_dofs(problem1)
-    assign_global_dofs!(problem1, dofmap)
-    assign_global_dofs!(problem2, dofmap)
+    ga1 = initialize_global_assembly(p1)
+    calculate_global_assembly!(ga1, p1)
+    ga2 = initialize_global_assembly(p2)
+    calculate_global_assembly!(ga2, p2)
 
-    # assemble problem 1
-    A1 = sparse(get_lhs(problem1, t)...)
-    b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1))
-
-    # assemble problem 2
-    A2 = sparse(get_lhs(problem2, t)...)
-    b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1))
-
-    # make one monolithic assembly
+    A1 = ga1.stiffness_matrix
+    b1 = ga1.force_vector
+    A2 = ga2.stiffness_matrix
+    b2 = ga2.force_vector
+    
+    # create a saddle point problem
     A = [A1 A2; A2' zeros(A2)]
     b = [b1; b2]
 
     # solve problem
-    nz = unique(rowvals(A))
+    nz = unique(rowvals(A))  # here we remove any zero rows
     x = zeros(b)
     x[nz] = lufact(A[nz,nz]) \ full(b[nz])
 
@@ -68,26 +139,23 @@ function call(solver::SimpleSolver, t)
     x1 = x[1:length(b1)]
     x2 = x[length(b1)+1:end]
 
-    # check residual
-    R1 = A1*x1 - b1
-    R2 = A2*x2 - b2
-    println("Residual norm: $(norm(R1+R2))")
-
     # update field for elements in problem 1
-    for equation in get_equations(problem1)
-        gdofs = get_global_dofs(equation)
+    for equation in get_equations(p1)
         element = get_element(equation)
-        field_name = get_unknown_field_name(equation)  # field we are solving
-        field = Field(t, full(x1[gdofs])[:])
+        field_name = get_unknown_field_name(p1)
+        gdofs = get_gdofs(p1, equation)
+        element_solution = full(x1[gdofs])
+        field = Field(time, element_solution)
         push!(element[field_name], field)
     end
 
-    # update field for elements in problem 2
-    for equation in get_equations(problem2)
-        gdofs = get_global_dofs(equation)
+    # update field for elements in problem 2 (Dirichlet boundary)
+    for equation in get_equations(p2)
         element = get_element(equation)
-        field_name = get_unknown_field_name(equation)
-        field = Field(t, full(x2[gdofs]))
+        field_name = get_unknown_field_name(p2)
+        gdofs = get_gdofs(p2, equation)
+        element_solution = full(x2[gdofs])
+        field = Field(time, element_solution)
         push!(element[field_name], field)
     end
 end

src/types.jl

@@ -7,7 +7,7 @@ using ForwardDiff
 
 """ Field is a fundamental data type which holds some values in some time t """
 type Field{T}
-    time :: Float64
+    time :: Number
     increment :: Int64
     values :: T
 end
@@ -106,12 +106,15 @@ type FieldSet
     fields :: Array{Field, 1}
 end
 """ Initializer for FieldSet. """
-function FieldSet(field_name)
+function FieldSet(field_name::ASCIIString)
     FieldSet(field_name, [])
 end
 function FieldSet()
     FieldSet("unknown field", [])
 end
+function FieldSet(fields::Array{Field, 1})
+    FieldSet("unknown field", fields)
+end
 """ Add new field to fieldset. """
 function Base.push!(fs::FieldSet, field::Field)
     push!(fs.fields, field)
@@ -130,6 +133,9 @@ end
 function Base.endof(fieldset::FieldSet)
     length(fieldset)
 end
+function Base.convert(fieldset::Type{FieldSet}, field::Field)
+    FieldSet(Field[field])
+end
 
 
 """ Basis function. """

test/test_elasticity.jl

@@ -0,0 +1,26 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# unit tests for heat equations
+
+using FactCheck
+using JuliaFEM: Quad4, Field, FieldSet, CPS4, get_basis, solve!, PlaneStressElasticityProblem
+
+
+facts("test plane elasticity on single element, volume load") do
+    element = Quad4([1, 2, 3, 4])
+    element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]))
+    element["youngs modulus"] = FieldSet(Field(500.0))
+    element["poissons ratio"] = FieldSet(Field(0.3))
+    element["displacement load"] = FieldSet(Field(0.0, Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]]))
+    equation = CPS4(element)
+    free_dofs = [3, 4, 5, 6]
+    problem = PlaneStressElasticityProblem([equation])
+    solve!(problem, free_dofs; max_iterations=10)
+    #solve!(equation, "displacement", free_dofs; max_iterations=10)
+    disp = get_basis(element)("displacement", [1.0, 1.0])[2]
+    Logging.info("displacement at tip: $disp")
+    # verified using Code Aster.
+    @fact disp --> roughly(-8.77303119819776E+00)
+end
+

test/test_elements.jl

@@ -34,12 +34,14 @@ end
 
 facts("test adding fieldsets and fields to element") do
     el = MockElement([1, 2, 3, 4])
+ 
     fieldset = JuliaFEM.FieldSet("geometry")
     field1 = JuliaFEM.Field(0.0, [0.0, 0.0, 0.0, 0.0])
     push!(fieldset, field1)
     field2 = JuliaFEM.Field(1.0, [1.0, 1.0, 1.0, 1.0])
     push!(fieldset, field2)
-    push!(el, fieldset)
+
+    el["geometry"] = fieldset
     fields = el["geometry"]
     @fact length(fields) --> 2
     @fact fields[1] --> field1
@@ -57,13 +59,13 @@ facts("interpolation of fields in some function space") do
     fieldset6 = FieldSet("vector field 3", [Field(0.0, Vector[[1.0, 5.0, 9.0], [2.0, 6.0, 10.0], [3.0, 7.0, 11.0], [4.0, 8.0, 12.0]])])
     fieldset7 = FieldSet("tensor field 1", [Field(0.0, Matrix[[1.0 5.0; 9.0 13.0], [2.0 6.0; 10.0 14.0], [3.0 7.0; 11.0 15.0], [4.0 8.0; 12.0 16.0]])])
 
-    push!(element, fieldset1)
-    push!(element, fieldset2)
-    push!(element, fieldset3)
-    push!(element, fieldset4)
-    push!(element, fieldset5)
-    push!(element, fieldset6)
-    push!(element, fieldset7)
+    element["geometry"] = fieldset1
+    element["constant scalar field"] = fieldset2
+    element["scalar field"] = fieldset3
+    element["vector field 1"] = fieldset4
+    element["vector field 2"] = fieldset5
+    element["vector field 3"] = fieldset6
+    element["tensor field 1"] = fieldset7
 
     xi = [0.0, 0.0]
     t = 0.0

test/test_global_assembly.jl

@@ -0,0 +1,32 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM: Quad4, Seg2, FieldSet, Field, PlaneHeatProblem
+using JuliaFEM: initialize_global_assembly, calculate_global_assembly!
+using FactCheck
+
+facts("assemble a simple two element problem and solve") do 
+    el1 = Quad4([1, 2, 3, 4])
+    el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    el1["temperature thermal conductivity"] = 6.0
+    el1["temperature load"] = [12.0, 12.0, 12.0, 12.0]
+    el1["density"] = 10
+
+    el2 = Seg2([1, 2])
+    el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
+
+    # Boundary load, linear ramp 0 -> 600 at time 0 -> 1
+    el2["temperature flux"] = FieldSet(Field[Field(0.0, 0.0), Field(1.0, 600.0)])
+
+    problem = PlaneHeatProblem()
+    push!(problem, el1)
+    push!(problem, el2)
+
+    global_assembly = initialize_global_assembly(problem)
+    calculate_global_assembly!(global_assembly, problem)
+    free_dofs = [1, 2]
+    A = lufact(global_assembly.stiffness_matrix[free_dofs, free_dofs])
+    b = full(global_assembly.force_vector)[free_dofs]
+    u = A \ b
+    @fact u --> roughly([101.0, 101.0])
+end

test/test_heat.jl

@@ -0,0 +1,43 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# unit tests for heat equations
+
+using FactCheck
+using JuliaFEM: Seg2, Quad4, Field, FieldSet, DC2D4, initialize_local_assembly, calculate_local_assembly!, DC2D2
+
+facts("tests on [0x1]x[0x1] domain") do
+
+    # volume element
+    element = Quad4([1, 2, 3, 4])
+    element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))
+    element["temperature thermal conductivity"] = FieldSet(Field(0.0, 6.0))
+    element["temperature load"] = FieldSet(Field(0.0, [12.0, 12.0, 12.0, 12.0]))
+    element["density"] = FieldSet(Field(0.0, 36.0))
+
+    # boundary element 
+    boundary_element = Seg2([1, 2])
+    boundary_element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [1.0, 0.0]]))
+    # linear ramp from 1 to 6 in time 0 to 1
+    boundary_element["temperature flux"] = FieldSet(Field[Field(0.0, 0.0), Field(1.0, 6.0)])
+
+    # Set constant source f=12 with k=6. Accurate solution is
+    # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
+    equation = DC2D4(element)
+    la = initialize_local_assembly()
+    calculate_local_assembly!(la, equation, "temperature")
+    fdofs = [1, 2]
+    A = la.stiffness_matrix
+    b = la.force_vector
+    @fact A[fdofs, fdofs] \ b[fdofs] --> roughly([1.0, 1.0])
+
+    # Set constant flux g=6 on boundary. Accurate solution is
+    # u(x,y) = x which equals T=1 on boundary.
+    boundary_equation = DC2D2(boundary_element);
+
+    calculate_local_assembly!(la, boundary_equation, "temperature")
+    b = la.force_vector
+    @fact A[fdofs, fdofs] \ b[fdofs] --> roughly([1.0, 1.0])
+
+end
+

test/test_interpolation.jl

@@ -4,28 +4,17 @@
 using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4
 using FactCheck
 
-element = Quad4([1, 2, 3, 4])
-
-geometry_field = Field(0.0, Vector[])  # Create empty field at time t=0.0
-push!(geometry_field, [ 0.0, 0.0])  # push some values for field
-push!(geometry_field, [ 1.0, 0.0])
-push!(geometry_field, [ 1.0, 1.0])
-push!(geometry_field, [ 0.0, 1.0])
-geometry_fieldset = FieldSet("geometry")  # create fieldset "geometry"
-push!(geometry_fieldset, geometry_field)  # add field to fieldset
-push!(element, geometry_fieldset) # add fieldset to element
-
-temperature_fieldset = FieldSet("temperature")
-push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))
-push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))
-push!(element, temperature_fieldset)
-
-displacement_fieldset = FieldSet("displacement")
-push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))
-push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))
-push!(element, displacement_fieldset)
 
 facts("basic continuum interpolations") do
+
+    element = Quad4([1, 2, 3, 4])
+
+    element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0])
+    element["displacement"] = (
+        Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]],
+        Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
+
     # from my old home works
     basis = get_basis(element)
     dbasis = grad(basis)

test/test_solver.jl

@@ -0,0 +1,71 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# test SimpleSolver
+
+using FactCheck
+using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis
+
+""" Define Problem 1:
+
+- Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}
+- Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁
+"""
+function get_heatproblem()
+    el1 = Quad4([1, 2, 3, 4])
+    # these might look like normal values but believe me, they
+    # are fields with temporal and spatial dimension
+    el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    el1["temperature thermal conductivity"] = 6.0
+    el1["density"] = 36.0
+
+    el2 = Seg2([1, 2])
+    el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
+    # Boundary load, linear ramp 0 -> 600 at time 0 -> 1
+    # yet another simplification, if field is given as a tuple,
+    # multiple fields are created. there is 1 second time step between
+    # each field. So the following is basically same as
+    # fieldset = FieldSet("temperature flux")
+    # field1 = Field(0.0, 0.0)
+    # field2 = Field(1.0, 600.0)
+    # push!(fieldset, field1)
+    # push!(fieldset, field2)
+    # element["temperature flux"] = fieldset
+    el2["temperature flux"] = (0.0, 600.0)
+
+    problem1 = PlaneHeatProblem()
+    push!(problem1, el1)
+    push!(problem1, el2)
+    return problem1
+end
+
+""" Define Problem 2:
+ - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂
+"""
+function get_boundaryproblem()
+    el3 = Seg2([3, 4])
+    el3["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
+    problem2 = DirichletProblem(1)
+    push!(problem2, el3)
+    return problem2
+end
+
+facts("test simplesolver") do
+    problem1 = get_heatproblem()
+    problem2 = get_boundaryproblem()
+    # Create a solver for a set of problems
+    solver = SimpleSolver()
+    push!(solver, problem1)
+    push!(solver, problem2)
+    # Solve problem at time t=1.0 and update fields
+    call(solver, 1.0)
+    # Postprocess.
+    # Interpolate temperature field along boundary of Γ₁ at time t=1.0
+    xi = [0.0, -1.0]
+    el2 = get_element(problem1.equations[2])
+    basis = get_basis(el2)
+    X = basis("geometry", xi, 1.0)
+    T = basis("temperature", xi, 1.0)
+    Logging.info("Temperature at point X = $X is T = $T")
+    @fact T --> roughly(100.0)
+end