diff --git a/src/assembly.jl b/src/assembly.jl index 57907bb..692a2f6 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -113,52 +113,3 @@ function assemble!(problem::Problem, time::Real, ::Type{Val{:mass_matrix}}; dens end end -# Static condensation routines - -function eliminate_interior_dofs(K::SparseMatrixCSC, f::SparseMatrixCSC, B::Vector{Int64}, I::Vector{Int64}; F=nothing, chunk_size=100000) - - dim = size(K, 1) - Kib = K[I,B] - - if F == nothing - F = cholfact(1/2*(K + K')[I,I]) - end - - if dim < chunk_size - # for small problems we don't need to care about memory usage - Kd = Kib' * (F \ Kib) - else - # for larger problems calculate schur complement in pieces - nb = length(B) - p = nb > 10 ? round(Int, nb/10) : nb - Kd = zeros(nb, nb) - for bi in 1:nb - done = round(Int, bi/nb*100) - mod(bi, p) == 0 && info("Static condensation: $done % done") - C = full(F \ Kib[:, bi]) - for bj in bi:nb - d = Kib[:, bj] - @inbounds Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d)) - end - end - Kd += tril(Kd, -1)' - end - - Kc = spzeros(dim, dim) - Kc[B,B] = K[B,B] - Kd - - fc = spzeros(dim, 1) - fc[B] = f[B] - Kib' * (F \ f[I]) - - return Kc, fc -end - -#= -function reconstruct!(ca::CAssembly, x::SparseMatrixCSC) - if isa(ca.F, Factorization) - x[ca.interior_dofs] = ca.F \ (ca.fi - ca.Kib*x[ca.boundary_dofs]) - else # normal inverse of matrix - x[ca.interior_dofs] = ca.F * (ca.fi - ca.Kib*x[ca.boundary_dofs]) - end -end -=# diff --git a/test/test_assembly.jl b/test/test_assembly.jl deleted file mode 100644 index 7fc9085..0000000 --- a/test/test_assembly.jl +++ /dev/null @@ -1,49 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "test static condensation" begin - - K = sparse([ - 4.0 -1.0 -2.0 -1.0 - -1.0 4.0 -1.0 -2.0 - -2.0 -1.0 4.0 -1.0 - -1.0 -2.0 -1.0 4.0]) - - f = sparse([6.0, 6.0, 3.0, 3.0]) - - I = [1, 2] - B = [3, 4] - - Kc, fc = eliminate_interior_dofs(K, f, B, I) - Kc = full(Kc) - fc = full(fc) - dump(Kc) - dump(fc) - - Kc_expected = [ - 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 - 0.0 0.0 2.4 -2.4 - 0.0 0.0 -2.4 2.4] - fc_expected = [0.0, 0.0, 9.0, 9.0] - - # TODO: needs to check numbers - @test isapprox(Kc, Kc_expected) - @test isapprox(fc, fc_expected) - - #= - x = sparse(zeros(4))' - la = sparse(zeros(4))' - la[3] = la[4] = 24.0 - reconstruct!(cass, x) - x = full(x) - info(la) - info(x) - @test isapprox(x[1], 1.0) - @test isapprox(x[2], 1.0) - =# -end -