fixes

src/core.jl

@@ -69,3 +69,6 @@ include("directsolver.jl") # parallel sparse direct solver for non-linear proble
 include("mortar.jl")  # mortar projection
 
 include("abaqus_reader_old.jl")
+
+# rest of things
+include("utils.jl")

src/dirichlet.jl

@@ -2,11 +2,16 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 type Dirichlet <: BoundaryProblem
+    formulation :: Symbol
     dual_basis :: Bool
 end
 
 function Dirichlet()
-    Dirichlet(true)
+    Dirichlet(:Equality, true)
+end
+
+function get_unknown_field_name(::Type{Dirichlet})
+    return "reaction force"
 end
 
 function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time::Real)
@@ -14,8 +19,8 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     @assert problem.properties.dual_basis
 
     # get dimension and name of PARENT field
-    field_dim = problem.dimension
-    field_name = problem.parent_field_name
+    field_dim = get_unknown_field_dimension(problem)
+    field_name = get_parent_field_name(problem)
     gdofs = get_gdofs(element, field_dim)
 
     # calculate bi-orthogonal basis transformation matrix Ae

src/elasticity.jl

@@ -69,7 +69,7 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El
         for dim in 1:get_unknown_field_dimension(problem)
             if haskey(element, "displacement traction force $dim")
                 T = element("displacement traction force $dim", ip, time)
-                ldofs = gdofs[dim:problem.dim:end]
+                ldofs = gdofs[dim:get_unknown_field_dimension(problem):end]
                 L = w*T*N*norm(J)
                 add!(assembly.f, ldofs, vec(L))
             end

src/elements.jl

@@ -6,6 +6,7 @@ abstract AbstractElement
 type Element{E}
     connectivity :: Vector{Int}
     fields :: Dict{ASCIIString, Field}
+    dualbasis :: Matrix{Float64} # coefficients to construct dual basis
 end
 
 function Base.size{E}(::Element{E})
@@ -18,7 +19,7 @@ end
 
 function convert{E}(::Type{Element{E}}, connectivity::Vector{Int})
 #   return Element{E}(connectivity, get_integration_points(E), Dict())
-    return Element{E}(connectivity, Dict())
+    return Element{E}(connectivity, Dict(), Matrix())
 end
 
 function get_integration_points{E}(element::Element{E}, args...)
@@ -183,6 +184,25 @@ function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0.0)
     return get_basis(element, xi)
 end
 
+function call(element::Element, xi::VecOrIP, time::Real, ::Type{Val{:dualbasis}})
+    if length(element.dualbasis) == 0
+        nnodes = size(element, 2)
+        De = zeros(nnodes, nnodes)
+        Me = zeros(nnodes, nnodes)
+        for ip in get_integration_points(element, Val{3})
+            J = get_jacobian(element, ip, time)
+            w = ip.weight*norm(J)
+            N = element(ip, time)
+            De += w*diagm(vec(N))
+            Me += w*N'*N
+        end
+        element.dualbasis = De*inv(Me)
+    end
+    N = get_basis(element, xi)
+    Phi = element.dualbasis*N'
+    return Phi'
+end
+
 function get_basis{E}(element::Element{E})
     basis = CVTI(
         (xi::Vector) -> get_basis(E, xi),
@@ -248,6 +268,7 @@ end
 """ Calculate local normal-tangential coordinates for element. """
 function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real)
     ntcoords = Matrix[]
+    normals = Vector{Float64}[]
     refcoords = get_reference_element_coordinates(E)
     x = element("geometry", time)
     for xi in refcoords
@@ -257,6 +278,7 @@ function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::
         if m == 1 # plane case
             tangent = dN / norm(dN)
             normal = [-tangent[2] tangent[1]]'
+            push!(normals, vec(normal))
             push!(ntcoords, [normal tangent])
         elseif m == 2
             normal = cross(dN[:,1], dN[:,2])
@@ -270,16 +292,38 @@ function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::
             tangent1 = u2/norm(u2)
             tangent2 = u3/norm(u3)
             push!(ntcoords, [normal tangent1 tangent2])
+            push!(normals, vec(normal))
         else
             error("calculate_normal_tangential_coordinates!(): n=$n, m=$m")
         end
     end
     element["normal-tangential coordinates"] = ntcoords
+    element["normals"] = normals
 end
-function calculate_normal_tangential_coordinates!{E}(elements::Vector{Element{E}}, time::Real)
+
+""" Calculate normal-tangential coordinates for a set of elements. 
+
+Notes
+-----
+Average normals so that normals are unique in nodes.
+"""
+function calculate_normal_tangential_coordinates!(elements::Vector, time::Real)
     for element in elements
         calculate_normal_tangential_coordinates!(element, time)
     end
+    n = calculate_nodal_vector("normals", 2, elements, time)
+    n = reshape(n, 2, round(Int, length(n)/2))
+    t = zeros(n)
+    for node_id=1:size(n,2)
+        n[:,node_id] = n[:,node_id] / norm(n[:,node_id])
+        t[:,node_id] = [-n[2,node_id], n[1,node_id]]
+    end
+    for element in elements
+        node_ids = get_connectivity(element)
+        Q = Matrix{Float64}[ [n[:,node_id] t[:,node_id]] for node_id in node_ids]
+        element["normal-tangential coordinates"] = Q
+        element["normals"] = Vector{Float64}[n[:,node_id] for node_id in node_ids]
+    end
 end
 
 """ Update element field based on a dictionary of nodal data and connectivity information.

src/mortar.jl

@@ -643,6 +643,15 @@ end
 # Mortar assembly 2d
 
 type Mortar <: BoundaryProblem
+    formulation :: Symbol
+end
+
+function Mortar()
+    Mortar(:Equality)
+end
+
+function get_unknown_field_name(::Type{Mortar})
+    return "reaction force"
 end
 
 typealias MortarElements2D Union{Seg2, Seg3}
@@ -711,8 +720,8 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor
 
             add!(assembly.C1, slave_dofs, slave_dofs, S2)
             add!(assembly.C1, slave_dofs, master_dofs, -M2)
-            S2 = Q2*S2
-            M2 = Q2*M2
+            S2 = Q2'*S2
+            M2 = Q2'*M2
             add!(assembly.C2, slave_dofs, slave_dofs, S2)
             add!(assembly.C2, slave_dofs, master_dofs, -M2)
 

src/problems.jl

@@ -22,9 +22,14 @@ type Assembly
     D :: SparseMatrixCOO 
     g :: SparseMatrixCOO
 
-    solution :: Vector{Float64}  # full solution vector when solving problem Ax = b
-    previous_solution :: Vector{Float64}  # previous solution vector
-    solution_norm_change :: Real  # for convergence studies
+    u :: Vector{Float64}  # solution vector u
+    u_prev :: Vector{Float64}  # previous solution vector u
+    u_norm_change :: Real  # change of norm in u
+
+    la :: Vector{Float64}  # solution vector la
+    la_prev :: Vector{Float64}  # previous solution vector u
+    la_norm_change :: Real # change of norm in la
+
     prehooks :: Vector{Tuple{Symbol,Any,Any}}  # assign possible prehooks before assembly
     posthooks :: Vector{Tuple{Symbol,Any,Any}}  # assign possible posthooks after assembly
     changed :: Bool  # flag to control is reassembly needed
@@ -40,6 +45,7 @@ function Assembly()
         SparseMatrixCOO(),
         SparseMatrixCOO(),
         [], [], Inf,
+        [], [], Inf,
         [], [], true)
 end
 
@@ -101,22 +107,87 @@ function get_assembly(problem)
     return problem.assembly
 end
 
-""" Update problem solution vector.
+""" Initialize unknown field ready for nonlinear iterations, i.e.,
+    take last known value and set it as a initial quess for next
+    time increment.
 """
-function update!(problem::Problem, solution::Vector{Float64})
+function initialize!(problem::Problem, time::Real)
+    field_name = get_unknown_field_name(problem)
+    field_dim = get_unknown_field_dimension(problem)
+    for element in get_elements(problem)
+        gdofs = get_gdofs(element, problem)
+        if haskey(element, field_name)
+            # if field is found, copy last known solution to new time as initial guess
+            if !isapprox(last(element[field_name]).time, time)
+                last_data = copy(last(element[field_name]).data)
+                push!(element[field_name], time => last_data)
+            end
+        else # if field not found at all, initialize new zero field.
+            data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
+            element[field_name] = (time => data)
+        end
+    end
+end
+
+""" Update problem solution vector for assembly. """
+function update_assembly!(problem, u, la)
+
     assembly = get_assembly(problem)
-    # resize & fill with zeros solution vector if length mismatch with current solution
-    if length(solution) != length(assembly.solution)
-        resize!(assembly.solution, length(solution))
-        fill!(assembly.solution, 0.0)
+
+    # resize & fill with zeros vectors if length mismatch with current solution
+    if length(u) != length(assembly.u)
+        resize!(assembly.u, length(u))
+        fill!(assembly.u, 0.0)
     end
-    assembly.previous_solution = copy(assembly.solution)
+    if length(la) != length(assembly.la)
+        resize!(assembly.la, length(la))
+        fill!(assembly.la, 0.0)
+    end
+
+    # copy current solutions to previous ones and add/replace new solution
+    assembly.u_prev = copy(assembly.u)
+    assembly.la_prev = copy(assembly.la)
     if get_formulation_type(problem) == :incremental
-        assembly.solution += solution
+        assembly.u += u
+        assembly.la += la
     else
-        assembly.solution = solution
+        assembly.u = u
+        assembly.la = la
+    end
+
+    # calculate change of norm
+    assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
+    assembly.la_norm_change = norm(assembly.la - assembly.la_prev)
+    return assembly.u_norm_change, assembly.la_norm_change
+end
+
+""" Update solutions to elements.
+
+Notes
+-----
+This assumes that element is properly initialized so that last known field data
+is from current time. For boundary problems solution is updated from lambda vector
+and for field problems from actual solution vector.
+"""
+function update_elements!(problem, u, la)
+    field_name = get_unknown_field_name(problem)
+    field_dim = get_unknown_field_dimension(problem)
+    nnodes = round(Int, length(u)/field_dim)
+
+    solution = nothing
+    if is_field_problem(problem)
+        solution = reshape(u, field_dim, nnodes)
+    elseif is_boundary_problem(problem)
+        solution = reshape(la, field_dim, nnodes)
+    else
+        error("update_elements!(): unknown problem type $(typeof(problem))")
+    end
+
+    for element in get_elements(problem)
+        connectivity = get_connectivity(element) # node ids
+        local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity]
+        last(element[field_name]).data = local_sol
     end
-    assembly.solution_norm_change = norm(assembly.solution - assembly.previous_solution)
 end
 
 #=
@@ -144,52 +215,12 @@ function get_unknown_field_name{P}(problem::Problem{P})
     return get_unknown_field_name(P)
 end
 
+""" Return the name of the parent field of this (boundary) problem. """
+function get_parent_field_name{P<:BoundaryProblem}(problem::Problem{P})
+    return problem.parent_field_name
+end
+
 function push!(problem::Problem, element)
     push!(problem.elements, element)
 end
 
-
-
-# TODO: better place for utility functions?
-
-""" Calculate "nodal" vector from set of elements.
-
-For example element 1 with dofs [1, 2, 3, 4] has [1, 1, 1, 1] and
-element 2 with dofs [3, 4, 5, 6] has [2, 2, 2, 2] the result will
-be sparse matrix with values [1, 1, 3, 3, 2, 2].
-
-Parameters
-----------
-field_name
-    name of field, e.g. "geometry"
-field_dim
-    degrees of freedom / node
-elements
-    elements used to calculate vector
-time
-"""
-function calculate_nodal_vector(field_name::ASCIIString, field_dim::Int, elements::Vector{Element}, time::Real)
-    A = SparseMatrixCOO()
-    b = SparseMatrixCOO()
-    for element in elements
-        haskey(element, field_name) || continue
-        gdofs = get_gdofs(element, 1)
-        for ip in get_integration_points(element, Val{2})
-            J = get_jacobian(element, ip, time)
-            w = ip.weight*norm(J)
-            f = element(field_name, ip, time)
-            N = element(ip, time)
-            add!(A, gdofs, gdofs, w*kron(N', N))
-            for dim=1:field_dim
-                add!(b, gdofs, w*f[dim]*N, dim)
-            end
-        end
-    end
-    A = sparse(A)
-    b = sparse(b)
-    nz = sort(unique(rowvals(A)))
-    x = zeros(size(b)...)
-    x[nz, :] = A[nz,nz] \ b[nz, :]
-    return vec(transpose(x))
-end
-

src/solvers.jl

@@ -265,7 +265,7 @@ end
 
 """ Solve linear system using LU factorization (UMFPACK).
 """
-function solve_linear_system!(solver::Solver, ::Type{Val{:DirectLinearSolver}})
+function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
     info("solving linear system of $(length(solver.problems)) problems.")
     t0 = time()
 
@@ -290,15 +290,10 @@ function solve_linear_system!(solver::Solver, ::Type{Val{:DirectLinearSolver}})
     x = zeros(length(b))
     x[nz1] = lufact(A[nz1,nz2]) \ full(b[nz1])
 
-    # update solutions
     u = x[1:dim]
     la = x[dim+1:end]
-    for problem in solver.problems
-        is_field_problem(problem) && update!(problem, u)
-        is_boundary_problem(problem) && update!(problem, la)
-    end
-
     info("UMFPACK: solved in ", time()-t0, " seconds. norm = ", norm(u))
+    return u, la
 end
 
 """ Check convergence of problems.
@@ -307,12 +302,16 @@ Notes
 -----
 Default convergence criteria is obtained by checking each sub-problem convergence.
 """
-function has_converged(solver::Solver; print_convergence_information=true)
+function has_converged(solver::Solver)
     converged = true
+    eps = solver.nonlinear_system_convergence_tolerance
     for problem in solver.problems
-        has_converged = problem.assembly.solution_norm_change < solver.nonlinear_system_convergence_tolerance
-        if print_convergence_information
-            @printf "% 30s | %8.3f | %s\n" problem.name problem.assembly.solution_norm_change has_converged
+        has_converged = true
+        if is_field_problem(problem)
+            has_converged = problem.assembly.u_norm_change/norm(problem.assembly.u) < eps
+        end
+        if is_boundary_problem(problem)
+            has_converged = problem.assembly.la_norm_change/norm(problem.assembly.la) < eps
         end
         converged &= has_converged
     end
@@ -328,82 +327,6 @@ function Base.showerror(io::IO, exception::NonlinearConvergenceError)
     print(io, "nonlinear iteration did not converge in $max_iters iterations!")
 end
 
-""" Initialize unknown field ready for nonlinear iterations, i.e.,
-    take last known value and set it as a initial quess for next
-    time increment.
-"""
-function initialize!{P<:FieldProblem}(problem::Problem{P}, time::Real)
-    field_name = get_unknown_field_name(problem)
-    field_dim = get_unknown_field_dimension(problem)
-    for element in get_elements(problem)
-        gdofs = get_gdofs(element, problem)
-        if haskey(element, field_name)
-            if !isapprox(last(element[field_name]).time, time)
-                last_data = copy(last(element[field_name]).data)
-                push!(element[field_name], time => last_data)
-            end
-        else # if field not found at all, initialize new zero field.
-            data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
-            element[field_name] = (time => data)
-        end
-    end
-end
-
-function initialize!{P<:BoundaryProblem}(problem::Problem{P}, time::Real; initialize_primary_field=false)
-    field_name = problem.parent_field_name
-    field_dim = problem.dimension
-
-    for element in get_elements(problem)
-        gdofs = get_gdofs(element, problem)
-        data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
-
-        # add new field "reaction force" for boundary element if not found
-        if haskey(element, "reaction force")
-            if !isapprox(last(element["reaction force"]).time, time)
-                push!(element["reaction force"], time => data)
-            end
-        else
-            element["reaction force"] = (time => data)
-        end
-
-        if initialize_primary_field
-            # add new primary field for boundary element if not found
-            if haskey(element, field_name)
-                if !isapprox(last(element[field_name]).time, time)
-                    last_data = copy(last(element[field_name]).data)
-                    push!(element[field_name], time => last_data)
-                end
-            else
-                data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
-                element[field_name] = (time => data)
-            end
-        end
-    end
-end
-
-function update!{P<:FieldProblem}(problem::Problem{P}, solution::Vector, ::Type{Val{:elements}})
-    field_name = get_unknown_field_name(problem)
-    field_dim = get_unknown_field_dimension(problem)
-    for element in get_elements(problem)
-        gdofs = get_gdofs(element, problem)
-        local_sol = solution[gdofs]
-        local_sol = reshape(local_sol, field_dim, length(element))
-        local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
-        last(element[field_name]).data = local_sol
-    end
-end
-
-function update!{P<:BoundaryProblem}(problem::Problem{P}, solution::Vector, ::Type{Val{:elements}})
-    field_dim = get_unknown_field_dimension(problem)
-    for element in get_elements(problem)
-        gdofs = get_gdofs(element, field_dim)
-        local_sol = solution[gdofs]
-        local_sol = reshape(local_sol, field_dim, length(element))
-        local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
-        last(element["reaction force"]).data = local_sol
-    end
-end
-
 """ Main solver loop.
 """
 function call(solver::Solver)
@@ -421,11 +344,12 @@ function call(solver::Solver)
         end
 
         # 2.2 call solver for linearized system (default: direct lu factorization)
-        solve_linear_system!(solver, Val{solver.linear_system_solver})
+        u, la = solve_linear_system(solver, Val{solver.linear_system_solver})
 
         # 2.3 update solution back to elements
         for problem in solver.problems
-            update!(problem, problem.assembly.solution, Val{:elements})
+            update_assembly!(problem, u, la)
+            update_elements!(problem, u, la)
         end
 
         # 2.4 check convergence

src/utils.jl

@@ -0,0 +1,68 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+""" Calculate nodal vector from set of elements.
+
+For example element 1 with dofs [1, 2, 3, 4] has [1, 1, 1, 1] and
+element 2 with dofs [3, 4, 5, 6] has [2, 2, 2, 2] the result will
+be sparse matrix with values [1, 1, 3, 3, 2, 2].
+
+Parameters
+----------
+field_name
+    name of field, e.g. "geometry"
+field_dim
+    degrees of freedom / node
+elements
+    elements used to calculate vector
+time
+"""
+function calculate_nodal_vector(field_name::ASCIIString, field_dim::Int,
+                                elements::Vector{Element}, time::Real)
+    A = SparseMatrixCOO()
+    b = SparseMatrixCOO()
+    for element in elements
+        haskey(element, field_name) || continue
+        gdofs = get_gdofs(element, 1)
+        for ip in get_integration_points(element, Val{2})
+            J = get_jacobian(element, ip, time)
+            w = ip.weight*norm(J)
+            f = element(field_name, ip, time)
+            N = element(ip, time)
+            add!(A, gdofs, gdofs, w*kron(N', N))
+            for dim=1:field_dim
+                add!(b, gdofs, w*f[dim]*N, dim)
+            end
+        end
+    end
+    A = sparse(A)
+    b = sparse(b)
+    nz = sort(unique(rowvals(A)))
+    x = zeros(size(b)...)
+    x[nz, :] = A[nz,nz] \ b[nz, :]
+    return vec(transpose(x))
+end
+
+""" Collect normal-tangential coordinates to rotation matrix Q.
+"""
+function get_rotation_matrix(elements, time)
+    Q = Dict{Int64, Matrix{Float64}}()
+    for element in elements
+        node_ids = get_connectivity(element)
+        q = element("normal-tangential coordinates", time).data
+        for (qi, node_id) in zip(q, node_ids)
+            if haskey(Q, node_id)
+                @assert isapprox(Q[node_id], qi)
+            else
+                Q[node_id] = qi
+            end
+        end
+    end
+    R = SparseMatrixCOO()
+    for (k, q) in Q
+        dofs = [(k-1)*2+1, (k-1)*2+2]
+        add!(R, dofs, dofs, q)
+    end
+    return R
+end
+

test/hdf5broken_test_model.jl

@@ -1,138 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using FactCheck
-#using JuliaFEM
-using Logging
-using HDF5
-@Logging.configure(level=DEBUG)
-
-@doc """
-Create new field to model.
-
-Parameters
-----------
-field_type : Dict()
-    Target topology (model.model, model.nodes, model.elements,
-                     model.element_nodes, model.element_gauss
-field_name : str
-    Field name
-
-Returns
--------
-Dict
-    New field
-
-""" ->
-function new_field!(model, field_type, field_name; partition=1, time=0, increment=0)
-  h5write("$model.$partition.h5", "$time/$increment/$field_type/$field_name", Float64[])
-end
-
-
-@doc """Add new nodes to model.
-
-Parameters
-----------
-model : str
-  Path to model file
-nodes : Dict()
-  id => coords
-partition : int, optional
-  Partition number, defaults to 1
-time : float, optional
-  Time step, defaults to 0
-increment : float, optional
-  Increment number, defaults to 0
-
-Returns
--------
-None
-
-Notes
------
-Create new field "coords" to model if not found
-""" ->
-function add_nodes!(model, nodes; partition=1, time=0, increment=0)
-  coords = Float64[]
-  node_ids = Int64[]
-  for (nid, ncoords) in nodes
-    @debug("adding nid: ", nid, " with coords: ", ncoords)
-    @assert length(ncoords) == 3
-    for x in ncoords
-      push!(coords, x)
-    end
-    push!(node_ids, nid)
-  end
-  h5write("$model.$partition.h5", "$time/$increment/nodes/coords", coords)
-  h5write("$model.$partition.h5", "$time/$increment/nodes/node_ids", node_ids)
-end
-
-
-
-@doc """Return subset of nodes from model.
-
-Parameters
-----------
-node_ids : array, optional
-  List of node ids to return. If not given, return all nodes.
-partition : int, optional
-  Partition number, defaults to 1
-time : float, optional
-  Time step, defaults to 0
-increment : float, optional
-  Increment number, defaults to 0
-
-Returns
--------
-Dict()
-    id => coords
-""" ->
-function get_nodes(model, node_ids=[]; partition=1, time=0, increment=0)
-  subset = Dict()
-  all_node_coords = h5read("$model.$partition.h5", "$time/$increment/nodes/coords")
-  all_node_ids = h5read("$model.$partition.h5", "$time/$increment/nodes/node_ids")
-  @debug("all node coords: ", all_node_coords)
-  @debug("all node ids: ", all_node_ids)
-  if length(node_ids) == 0
-    node_ids = all_node_ids
-  end
-  dim = 3
-  for node_id in node_ids
-    @debug("fetching node ", node_id)
-    idx = findfirst(all_node_ids, node_id)
-    @debug("found node coords from idx ", idx)
-    subset[node_id] = all_node_coords[dim*(idx-1)+1:dim*(idx-1)+dim]
-  end
-  return subset
-end
-
-
-
-facts("create new field to model") do
-  model = tempname()
-  @debug("model file name", model)
-  new_field!(model, "nodes", "coords")
-  data = h5read("$model.1.h5", "0/0/nodes/coords")
-  @fact length(data) => 0
-end
-
-facts("add nodes to model") do
-  model = tempname()
-  nodes = Dict(1 => [1.0, 2.0, 3.0])
-  add_nodes!(model, nodes)
-  @fact h5read("$model.1.h5", "0/0/nodes/coords") => [1.0, 2.0, 3.0]
-  @fact h5read("$model.1.h5", "0/0/nodes/node_ids") => [1]
-end
-
-facts("get nodes from model") do
-  model = tempname()
-  nodes = Dict(1 => [1.0, 2.0, 3.0], 2 => [2.0, 3.0, 4.0])
-  add_nodes!(model, nodes)
-  subset = get_nodes(model, [1])
-  @fact subset[1] => [1.0, 2.0, 3.0]
-  subset = get_nodes(model, [2])
-  @fact subset[2] => [2.0, 3.0, 4.0]
-  subset = get_nodes(model, [1,2])
-  @fact subset[1] => [1.0, 2.0, 3.0]
-  @fact subset[2] => [2.0, 3.0, 4.0]
-end