diff --git a/src/solvers.jl b/src/solvers.jl index 490e4ec..3e1e03e 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -184,55 +184,6 @@ function get_boundary_assembly(solver::Solver) return K, C1, C2, D, f, g end -function resize!(A::SparseMatrixCSC, m::Int64, n::Int64) - (n == A.n) && (m == A.m) && return - @assert n >= A.n - @assert m >= A.m - append!(A.colptr, A.colptr[end]*ones(Int, m-A.m)) - A.n = n - A.m = m -end - -""" -Given C and g, construct new basis such that v = P*u + g - -Parameters ----------- -S set of linearly independent dofs. -""" -function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12) - n, m = size(C) - @assert n == m - if S == nothing - S = get_nonzero_rows(C) - end - # FIXME: this creates dense matrices - # efficiency / memory usage is a question - M = get_nonzero_columns(C) - F = qrfact(C[S,:]) - P = spzeros(n,m) - P[:,M] = sparse(F \ full(C[S,M])) - h = sparse(F \ full(g[S])) - resize!(P, n, m) - resize!(h, n, 1) - P = speye(n) - P - droptol!(P, tol) - return P, h -end - -""" Assume C is invertible. """ -function create_projection(C, g, ::Type{Val{:invertible}}) - nz1, nz2 = get_nonzeros(C) - P = spzeros(size(C)...) - for j=1:size(C,1) - j in nz1 && continue - P[j,j] = 1.0 - end - v = lufact(C[nz1,nz2]) \ full(g[nz1]) - return P, v -end - - """ Solve linear system using LDLt factorization (SuiteSparse). This version diff --git a/src/sparse.jl b/src/sparse.jl index d2639be..48c466d 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -178,11 +178,21 @@ function size(A::SparseMatrixCOO, idx::Int) return size(A)[idx] end +""" Resize sparse matrix A to (higher) dimension n x m. """ +function resize_sparse(A, n, m) + return sparse(findnz(A)..., n, m) +end + +""" Resize sparse vector b to (higher) dimension n. """ +function resize_sparsevec(b, n) + return sparsevec(findnz(b)..., n) +end + """ Matrix norm. Automatically convert to dense when asking for 2-norm for small matrices. """ function norm(A::SparseMatrixCOO, p=Inf; maxdim=1000) dim = size(A, 1) if p == 2 && dim > maxdim - info("Assembly norm: dim = $dim > $maxdim and p=$p, not making dense matrices for operation.") + warn("Assembly norm: dim = $dim > $maxdim and p=$p, not making dense matrices for operation.") return 0.0 end if p == 2 @@ -192,9 +202,9 @@ function norm(A::SparseMatrixCOO, p=Inf; maxdim=1000) end end +""" Approximative comparison of two matricse A and B. """ function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO) A2 = sparse(A) B2 = sparse(B, size(A2)...) return isapprox(A2, B2) end - diff --git a/test/test_projection.jl b/test/test_projection.jl deleted file mode 100644 index d7af575..0000000 --- a/test/test_projection.jl +++ /dev/null @@ -1,94 +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 -using JuliaFEM.Preprocess - -@testset "test projection" begin - C = [ - 2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 2.0 1.0 -1.0 -2.0 0.0 0.0 - 0.0 0.0 1.0 2.0 -2.0 -1.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] - g = [3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - P, h = create_projection(sparse(C), g) - P_expected = [ - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 - 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0] - h_expected = [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - @test isapprox(full(P), P_expected) - @test isapprox(full(h), h_expected) - nz = [5, 6, 7, 8] - info("is P'P positive definite? ", isposdef(P[:,nz]'P[:,nz])) - @test isposdef(P[:,nz]'P[:,nz]) -end - -@testset "test creating projection matrix from invertible problem" begin - # simple 3 element poisson problem - k = [1.0 -1.0; -1.0 1.0] - K = zeros(4, 4) - K[1:2,1:2] += k - K[2:3,2:3] += k - K[3:4,3:4] += k - # first dof homogeneous bc, last dof uâ‚„ = 1 - C = zeros(4, 4) - C[1,1] = 1.0 - C[4,4] = 1.0 - g = zeros(4) - g[1] = 0.0 - g[4] = 1.0 - P, h = create_projection(sparse(C), g, Val{:invertible}) - info("P = ") - dump(full(P)) - P_expected = zeros(4, 4) - P_expected[2,2] = P_expected[3,3] = 1.0 - @test isapprox(full(P), P_expected) - #h_expected = [0.5, 1.0] - #@test isapprox(h, h_expected) -end - -@testset "projection between surfaces" begin - # FIXME: creating mortar projection takes very long time. - meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/joint.med" - isfile(meshfile) || return - mesh = aster_read_mesh(meshfile, "JOINT") - - # top - bc1 = Problem(Dirichlet, "fixed1", 3, "displacement") - bc1.elements = create_elements(mesh, "FIXED1") - update!(bc1.elements, "displacement 1", 0.0) - update!(bc1.elements, "displacement 2", 0.0) - update!(bc1.elements, "displacement 3", 0.0) - - # bottom - bc2 = Problem(Dirichlet, "fixed2", 3, "displacement") - bc2.elements = create_elements(mesh, "FIXED2") - update!(bc2.elements, "displacement 1", 0.0) - update!(bc2.elements, "displacement 2", 0.0) - update!(bc2.elements, "displacement 3", 0.0) - - # joint - contact = Problem(Mortar, "joint", 3, "displacement") - master_elements = create_elements(mesh, "BODY1_TO_BODY2") - slave_elements = create_elements(mesh, "BODY2_TO_BODY1") - update!(slave_elements, "master elements", master_elements) - contact.elements = [master_elements; slave_elements] - - solver = Solver(Modal) - push!(solver, bc1, bc2, contact) - assemble!(solver) - Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) - P, h = create_projection(C1, g) -end -