From aa7d775ec138d1003d946a6f620870c46eb35e41 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 3 Feb 2016 20:39:03 +0200 Subject: [PATCH] fixes --- src/core.jl | 3 + src/dirichlet.jl | 11 ++- src/elasticity.jl | 2 +- src/elements.jl | 48 ++++++++++- src/mortar.jl | 13 ++- src/problems.jl | 147 ++++++++++++++++++++-------------- src/solvers.jl | 102 +++-------------------- src/utils.jl | 68 ++++++++++++++++ test/hdf5broken_test_model.jl | 138 ------------------------------- 9 files changed, 239 insertions(+), 293 deletions(-) create mode 100644 src/utils.jl delete mode 100644 test/hdf5broken_test_model.jl diff --git a/src/core.jl b/src/core.jl index 61dd385..50ceb89 100644 --- a/src/core.jl +++ b/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") diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 4207e72..e916bed 100644 --- a/src/dirichlet.jl +++ b/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 diff --git a/src/elasticity.jl b/src/elasticity.jl index d037abc..2f2f5bd 100644 --- a/src/elasticity.jl +++ b/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 diff --git a/src/elements.jl b/src/elements.jl index 802b4bd..c69e79d 100644 --- a/src/elements.jl +++ b/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. diff --git a/src/mortar.jl b/src/mortar.jl index c1ca1a2..646a13a 100644 --- a/src/mortar.jl +++ b/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) diff --git a/src/problems.jl b/src/problems.jl index a1c8bb1..dd3d275 100644 --- a/src/problems.jl +++ b/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 - diff --git a/src/solvers.jl b/src/solvers.jl index 59ec8cf..58b9de8 100644 --- a/src/solvers.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..bc9fc1e --- /dev/null +++ b/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 + diff --git a/test/hdf5broken_test_model.jl b/test/hdf5broken_test_model.jl deleted file mode 100644 index 4052fc9..0000000 --- a/test/hdf5broken_test_model.jl +++ /dev/null @@ -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