WIP: Improve sparse utils

src/JuliaFEM.jl

@@ -3,11 +3,19 @@
 
 # __precompile__()
 
+push!(LOAD_PATH, Pkg.dir("JuliaFEM") * "/subpackages")
+
 """
 This is JuliaFEM -- Finite Element Package
 """
 module JuliaFEM
 
+using FEMSparse
+import FEMSparse: get_nonzero_rows, get_nonzero_columns
+
+# Magic stuff to make life with SparseMatrixCSC easier
+include("sparse.jl")
+
 using TimerOutputs
 const to = TimerOutput()
 function print_statistics()
@@ -17,7 +25,7 @@ export print_statistics
 
 import Base: getindex, setindex!, convert, length, size, isapprox, similar,
              start, first, next, done, last, endof, vec, ==, +, -, *, /, haskey, copy,
-             push!, isempty, empty!, append!, sparse, full, read
+             push!, isempty, empty!, append!, sparse, full, read, resize!
 
 using Logging
 
@@ -38,7 +46,7 @@ module Testing
 end
 
 include("fields.jl")
-export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
+export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment, isconstant
 
 include("types.jl")  # data types: Point, IntegrationPoint, ...
 export AbstractPoint, Point, IntegrationPoint, IP, Node
@@ -47,7 +55,7 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 include("elements.jl") # common element routines
 export Node, AbstractElement, Element, update!, get_connectivity, get_basis,
        get_dbasis, inside, get_local_coordinates, get_element_type,
-       filter_by_element_type, get_element_id
+       filter_by_element_type, get_element_id, get_nodal_value
 
 include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
 export get_reference_coordinates,
@@ -69,9 +77,6 @@ export NSeg, NSurf, NSolid, is_nurbs
 include("integrate.jl")  # default integration points for elements
 export get_integration_points
 
-include("sparse.jl")
-export add!, SparseMatrixCOO, SparseVectorCOO, get_nonzero_rows, get_nonzero_columns, optimize!, resize_sparse, resize_sparsevec
-
 include("problems.jl") # common problem routines
 export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
        get_unknown_field_dimension, get_gdofs, Assembly,
@@ -115,7 +120,7 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
        get_field_problems, get_boundary_problems,
        get_field_assembly, get_boundary_assembly,
        initialize!, create_projection, eliminate_interior_dofs,
-       is_field_problem, is_boundary_problem
+       is_field_problem, is_boundary_problem, get_solution_vector
 include("solvers_modal.jl")
 export Modal
 
@@ -148,16 +153,27 @@ export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
 end
 
 module Postprocess
+
+using JuliaFEM
+using FEMSparse
+using DataFrames
+using HDF5
+using LightXML
+using Formatting
+
 include("postprocess_utils.jl")
 export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!,
        calculate_area, calculate_center_of_mass,
        calculate_second_moment_of_mass, extract
+
 end
 
 module Abaqus
 include("abaqus.jl")
 export abaqus_read_model, abaqus_run_model, abaqus_open_results, create_surface_elements
 end
 
+include("deprecated.jl")
+
 end
 

src/abaqus.jl

@@ -648,10 +648,8 @@ function (model::Model)()
         step_problems = [create_boundary_problem(model, bc) for bc in step.boundary_conditions]
         # 3.2 create solver and solve set of problems
         solver_type = determine_solver_type(model, step)
-        solver_description = "$solver_type solver"
-        solver = Solver(solver_type, solver_description)
         all_problems = [model.problems; boundary_problems; step_problems]
-        push!(solver, all_problems...)
+        solver = Solver(solver_type, all_problems...)
         solver()
         info(repeat("-", 80))
         info("Simulation ready, processing output requests")

src/deprecated.jl

@@ -0,0 +1,14 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Old ways to do things, we deprecate these in some schedule
+
+function Solver{S}(::Type{S}, name::String)
+    solver = Solver(S)
+    solver.name = name
+    return solver
+end
+
+function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
+    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P())
+end

src/elements.jl

@@ -62,6 +62,43 @@ function setindex!(element::Element, data, field_name)
     element.fields[field_name] = Field(data)
 end
 
+function hasfield(element::Element, field_name::String)
+    return haskey(element.fields, field_name)
+end
+
+function hasconnectivity(element::Element, node_id::Int64)
+    return node_id in element.connectivity
+end
+
+""" Extract a field from element. """
+function getfield(element::Element, field_name::String)
+    return element.fields[field_name]
+end
+
+""" Get field nodal values from element.
+
+Returns
+-------
+Dict, id => value
+
+"""
+function get_nodal_values(element::Element, field_name::String, time::Float64)
+    hasfield(element, field_name) || error("get_nodal_values(): field $field_name not defined in element")
+    field = getfield(element, field_name)
+    field_data = field(time)
+    connectivity = get_connectivity(element)
+    if isconstant(field)
+        return Dict(c => field_data for c in connectivity)
+    else
+        return Dict(c => field_data[i] for (i, c) in enumerate(connectivity))
+    end
+end
+
+function get_nodal_value(element::Element, field_name::String, time::Float64, node_id::Int64)
+    hasconnectivity(element, node_id) || error("get_nodal_value(): node $node_id not in element")
+    return get_nodal_values(element, field_name, time)[node_id]
+end
+
 """ Return a Field object from element.
 
 Examples

src/fields.jl

@@ -27,6 +27,14 @@ typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to inter
 typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time
 typealias CVTV Field{Continuous, Variable, TimeVariant}
 
+function isconstant(f::Field)
+    return false
+end
+
+function isconstant{A,B<:Constant,C}(f::Field{A,B,C})
+    return true
+end
+
 # Discrete fields
 
 

src/postprocess_utils.jl

@@ -1,17 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using JuliaFEM
-using DataFrames
-using HDF5
-using LightXML
-using Formatting
-
 """
 Calculate field values to nodal points from Gauss points using least-squares fitting.
 """
-function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
-                            F=nothing, nz=nothing, b=nothing, return_F_and_nz=false)
+function calc_nodal_values!(elements::Vector, field_name::String, time::Float64;
+                            F=nothing, nz=nothing, b=nothing)
 
     if F == nothing
         A = SparseMatrixCOO()
@@ -24,8 +18,8 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
                 add!(A, gdofs, gdofs, w*kron(N', N))
             end
         end
-        A = sparse(A)
         nz = get_nonzero_rows(A)
+        A = sparse(A)
         A = 1/2*(A + A')
         F = ldltfact(A[nz,nz])
     end
@@ -36,31 +30,24 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
             gdofs = get_connectivity(element)
             for ip in get_integration_points(element)
                 if !haskey(ip, field_name)
-                    info("warning: integration point does not have field $field_name")
+                    warn("integration point does not have field $field_name")
                     continue
                 end
                 detJ = element(ip, time, Val{:detJ})
                 w = ip.weight*detJ
                 f = ip(field_name, time)
                 N = element(ip, time)
-                for dim=1:field_dim
-                    add!(b, gdofs, w*f[dim]*N, dim)
-                end
+                add!(b, gdofs, collect(1:length(f)), w*(f*N)')
             end
         end
         b = sparse(b)
     end
 
     x = zeros(size(b)...)
     x[nz, :] = F \ b[nz, :]
-    nodal_values = Dict()
-    for i=1:size(x,1)
-        nodal_values[i] = vec(x[i,:])
-    end
+    nodal_values = Dict(i => vec(x[i,:]) for i=1:size(x,1))
     update!(elements, field_name, time => nodal_values)
-    if return_F_and_nz
-        return F, nz
-    end
+    return F, nz
 end
 
 """

src/problems.jl

@@ -6,85 +6,10 @@ abstract FieldProblem <: AbstractProblem
 abstract BoundaryProblem <: AbstractProblem
 abstract MixedProblem <: AbstractProblem
 
-"""
-General linearized problem to solve
-    (K₁+K₂)Δu  +   C1'*Δλ = f₁+f₂
-         C2Δu  +     D*Δλ = g
-"""
-type Assembly
-
-    M :: SparseMatrixCOO  # mass matrix
-
-    # for field assembly
-    K :: SparseMatrixCOO   # stiffness matrix
-    Kg :: SparseMatrixCOO  # geometric stiffness matrix
-    f :: SparseMatrixCOO   # force vector
-    fg :: SparseMatrixCOO  #
-
-    # for boundary assembly
-    C1 :: SparseMatrixCOO
-    C2 :: SparseMatrixCOO
-    D :: SparseMatrixCOO
-    g :: SparseMatrixCOO
-    c :: SparseMatrixCOO
-
-    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
-
-    removed_dofs :: Vector{Int64} # manually remove dofs from assembly
-end
-
-function Assembly()
-    return Assembly(
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        SparseMatrixCOO(),
-        [], [], Inf,
-        [], [], Inf,
-        [])
-end
-
-function empty!(assembly::Assembly)
-    empty!(assembly.K)
-    empty!(assembly.Kg)
-    empty!(assembly.f)
-    empty!(assembly.fg)
-    empty!(assembly.C1)
-    empty!(assembly.C2)
-    empty!(assembly.D)
-    empty!(assembly.g)
-    empty!(assembly.c)
-end
-
-function isempty(assembly::Assembly)
-    T = isempty(assembly.K)
-    T &= isempty(assembly.Kg)
-    T &= isempty(assembly.f)
-    T &= isempty(assembly.fg)
-    T &= isempty(assembly.C1)
-    T &= isempty(assembly.C2)
-    T &= isempty(assembly.D)
-    T &= isempty(assembly.g)
-    T &= isempty(assembly.c)
-    return T
-end
-
 type Problem{P<:AbstractProblem}
-    name :: AbstractString           # descriptive name for problem
+    name :: String                   # descriptive name for problem
     dimension :: Int                 # degrees of freedom per node
-    parent_field_name :: AbstractString # (optional) name of parent field e.g. "displacement"
+    parent_field_name :: String      # (optional) name of parent field e.g. "displacement"
     elements :: Vector{Element}
     dofmap :: Dict{Element, Vector{Int64}} # connects element local dofs to global dofs
     assembly :: Assembly
@@ -99,15 +24,20 @@ Examples
 --------
 Create vector-valued (dim=3) elasticity problem:
 
-julia> prob1 = Problem(Elasticity, "this is my problem", 3)
-julia> prob2 = Problem(Elasticity, 3)
+julia> prob1 = Problem(Elasticity, 3)
 
 """
-function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
-    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P())
-end
 function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64)
-    return Problem(P, "$P problem", dimension)
+    name = "$P problem"
+    properties = P()
+    parent_field_name = "none"
+    elements = Vector()
+    dofmap = Dict()
+    assembly = Assembly()
+    fields = Dict()
+    postprocess_fields = Vector()
+    return Problem{P}(name, dimension, parent_field_name, elements, dofmap,
+                      assembly, fields, postprocess_fields, properties)
 end
 
 """ Construct a new boundary problem.

src/problems_dirichlet.jl

@@ -73,9 +73,9 @@ function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
             end
         end
         for (k, v) in field_vals
-            push!(problem.assembly.C1, k, k, 1.0)
-            push!(problem.assembly.C2, k, k, 1.0)
-            push!(problem.assembly.g, k, 1, v)
+            add!(problem.assembly.C1, k, k, 1.0)
+            add!(problem.assembly.C2, k, k, 1.0)
+            add!(problem.assembly.g, k, v)
         end
     end