diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 40ef45c..46cbad8 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -116,8 +116,8 @@ using Reexport, ForwardDiff, LightXML, HDF5, Parameters using Tensors # For basis functions (Vec type) import Calculus # For symbolic differentiation in basis generation -import FEMSparse -# import FEMQuad # Consolidated into src/quadrature.jl +# import FEMSparse # Consolidated into src/sparse/ +# import FEMQuad # Consolidated into src/quadrature.jl # Note: Consolidating FEMBase and FEMBasis into JuliaFEM # Previously: @reexport using FEMBase diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 80affa6..07456b8 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -174,7 +174,7 @@ const X = ([-93.7197, -93.7197, 150.883], [-91.657, -85.8251, 157.885], [-100.52 const displacement_load_string = [string("displacement load ", i) for i in 1:3] """ Assemble 3d continuum elements in general solid mechanics problem. """ function assemble_element!(assembly::Assembly, - assembler::FEMSparse.AssemblerSparsityPattern, + assembler::AssemblerSparsityPattern, problem::Problem{Elasticity}, element::Element{El}, local_buffer::Elasticity3DLocalBuffers, @@ -387,10 +387,10 @@ function assemble_element!(assembly::Assembly, if use_csc # add contributions to K, Kg, f - @inbounds FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) + @inbounds assemble_local!(assembler, gdofs, Km, f_ext) if props.geometric_stiffness - @inbounds FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg) + @inbounds assemble_local_matrix!(assembler, gdofs, Kg) end else add!(assembly.f, gdofs, f_ext) diff --git a/src/sparse/sparse.jl b/src/sparse/sparse.jl index 148fc8c..cf0a812 100644 --- a/src/sparse/sparse.jl +++ b/src/sparse/sparse.jl @@ -1,9 +1,16 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE +# +# Sparse matrix utilities consolidated from FEMBase.jl and FEMSparse.jl using SparseArrays import SparseArrays: sparse, sparsevec +# SparseMatrixCOO from FEMBase (this file) +# AssemblerSparsityPattern from FEMSparse +include("sparsematrixcsc.jl") +# include("sparsevectordok.jl") # Old Julia syntax, not used, skipping for now + mutable struct SparseMatrixCOO{T<:Real} I :: Vector{Int} J :: Vector{Int} diff --git a/src/sparse/sparsematrixcsc.jl b/src/sparse/sparsematrixcsc.jl new file mode 100644 index 0000000..2265dde --- /dev/null +++ b/src/sparse/sparsematrixcsc.jl @@ -0,0 +1,171 @@ +struct AssemblerSparsityPattern{Tv,Ti} + K::SparseMatrixCSC{Tv,Ti} + f::Vector{Tv} + permutation::Vector{Int} + sorteddofs::Vector{Int} +end + +""" + start_assemble(K::SparseMatrixCSC, [f::Vector]) + +Create an `AssemblerSparsityPattern` from the sparsity pattern in `K` +and optionally a "force" vector. +""" +function start_assemble(K::SparseMatrixCSC, f = eltype(K)[]) + if !isempty(f) + @assert size(K, 1) == length(f) + end + AssemblerSparsityPattern(K, f, Int[], Int[]) +end + + +""" + assemble_local_vector!(f, dofs, fe) +Assembles the element residual `fe` into the global residual vector `f`. +""" +Base.@propagate_inbounds function assemble_local_vector!(f::AbstractVector, dofs::AbstractVector{Int}, fe::AbstractVector) + @boundscheck checkbounds(f, dofs) + @inbounds for i in 1:length(dofs) + f[dofs[i]] += fe[i] + end +end + + +function assemble_local!(A::AssemblerSparsityPattern, dofs::AbstractVector, Ke::AbstractMatrix, fe::AbstractVector) + @assert !isempty(A.f) + @boundscheck checkbounds(A.K, dofs, dofs) + @boundscheck checkbounds(A.f, dofs) + @inbounds assemble_local_vector!(A.f, dofs, fe) + @inbounds assemble_local_matrix!(A, dofs, Ke) +end + +""" + assemble!(A::AssemblerSparsityPattern, dofs2, Ke) + +Assemble a local dense element matrix `Ke` into the sparse matrix +wrapped by the assembler `A`. +location given by lists of indices `dofs1` and `dofs2`. + +# Example + +```julia +using SparseArrays +sparsity_pattern = sparse([1. 0 1; 1 0 1; 1 1 1]) +fill!(sparsity_pattern, 0) +assembler = FEMSparse.start_assemble(sparsity_pattern) +dofs = [1, 3] +Ke = [1.0 2.0; 3.0 4.0] +FEMSparse.FEMSparse.assemble_local_matrix!(assembler, dofs, Ke) +Matrix(sparsity_pattern) + +# output +julia> Matrix(sparsity_pattern) +3×3 Array{Float64,2}: + 1.0 0.0 2.0 + 0.0 0.0 0.0 + 3.0 0.0 4.0 +``` +""" +Base.@propagate_inbounds function assemble_local_matrix!(A::AssemblerSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix) + permutation = A.permutation + sorteddofs = A.sorteddofs + K = A.K + + @boundscheck checkbounds(K, dofs, dofs) + resize!(permutation, length(dofs)) + resize!(sorteddofs, length(dofs)) + copyto!(sorteddofs, dofs) + sortperm2!(sorteddofs, permutation) + + current_col = 1 + @inbounds for Kcol in sorteddofs + maxlookups = length(dofs) + current_idx = 1 + for r in nzrange(K, Kcol) + Kerow = permutation[current_idx] + if K.rowval[r] == dofs[Kerow] + Kecol = permutation[current_col] + K.nzval[r] += Ke[Kerow, Kecol] + current_idx += 1 + end + current_idx > maxlookups && break + end + if current_idx <= maxlookups + error("some row indices were not found") + end + current_col += 1 + end +end + + +################## +# Sort utilities # +################## +# We bundle a few sorting utilities here because the ones +# in have some unacceptable overhead. + +# Sorts B and stores the permutation in `ii` +function sortperm2!(B, ii) + @inbounds for i = 1:length(B) + ii[i] = i + end + quicksort!(B, ii) + return +end + +function quicksort!(A, order, i=1,j=length(A)) + @inbounds if j > i + if j - i <= 12 + # Insertion sort for small groups is faster than Quicksort + insertionsort!(A, order, i, j) + return A + end + + pivot = A[div(i+j,2)] + left, right = i, j + while left <= right + while A[left] < pivot + left += 1 + end + while A[right] > pivot + right -= 1 + end + if left <= right + A[left], A[right] = A[right], A[left] + order[left], order[right] = order[right], order[left] + + left += 1 + right -= 1 + end + end # left <= right + + quicksort!(A,order, i, right) + quicksort!(A,order, left,j) + end # j > i + + return A +end + +function insertionsort!(A, order, ii=1, jj=length(A)) + @inbounds for i = ii+1 : jj + j = i - 1 + temp = A[i] + itemp = order[i] + + while true + if j == ii-1 + break + end + if A[j] <= temp + break + end + A[j+1] = A[j] + order[j+1] = order[j] + j -= 1 + end + + A[j+1] = temp + order[j+1] = itemp + end # i + return +end diff --git a/src/sparse/sparsevectordok.jl b/src/sparse/sparsevectordok.jl new file mode 100644 index 0000000..01be819 --- /dev/null +++ b/src/sparse/sparsevectordok.jl @@ -0,0 +1,47 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/FEMSparse.jl/blob/master/LICENSE + +mutable struct SparseVectorDOK{Tv,Ti<:Integer} <: AbstractSparseArray{Tv,Ti,1} + data :: Dict{Ti,Tv} +end + +function SparseVectorDOK() + return SparseVectorDOK(Dict{Int64, Float64}()) +end + +function SparseVectorDOK{Tv,Ti<:Integer}(b::SparseVector{Tv,Ti}) + I, V = findnz(b) + c = SparseVectorDOK() + add!(c, I, V) + return c +end + +function SparseVectorDOK{T}(b::Vector{T}) + return SparseVectorDOK(sparsevec(b)) +end + +function add!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, i::Ti, v::Tv) + b.data[i] = get(b, i) + v + return nothing +end + +function add!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv}) + @assert length(dofs) == length(data) + z = Tv(0) + for i=1:length(dofs) + @inbounds b.data[dofs[i]] = Base.get(b.data, i, z) + data[i] + end + return nothing +end + +function get{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, i::Ti) + return Base.get(b.data, i, Tv(0)) +end + +function get!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv}) + for (i,j) in enumerate(dofs) + data[i] = get(b, i) + end + return nothing +end +