mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-22 18:52:16 +00:00
feat: Consolidate FEMSparse.jl into JuliaFEM
Consolidated FEMSparse package into src/sparse/:
- sparsematrixcsc.jl: AssemblerSparsityPattern for efficient assembly
- sparsevectordok.jl: Skipped (old Julia syntax, not used)
Changes:
- Removed 'import FEMSparse' from JuliaFEM.jl
- Updated problems_elasticity.jl: FEMSparse.AssemblerSparsityPattern → AssemblerSparsityPattern
- Fixed include paths in sparse/sparse.jl (relative, not absolute)
Modernization:
- Fixed sparsevectordok.jl type keyword → mutable struct
- Fixed AbstractSparseArray type parameter syntax
- Chose to skip sparsevectordok for now (old {T,V} syntax, unused)
Result:
- ✅ JuliaFEM loads successfully
- ✅ 5 tests still passing (no regression)
- ✅ Three vendor packages now consolidated: FEMBase, FEMBasis, FEMQuad, FEMSparse
Remaining vendor packages: AbaqusReader, AsterReader, HeatTransfer, FEMBeam, Mortar packages
This commit is contained in:
Jukka Aho
+2
-2
@@ -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
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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}
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user