mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
removed null space related code as obsolete
This commit is contained in:
Jukka Aho
@@ -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
|
||||
|
||||
+12
-2
@@ -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
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user