From f3995ee445d4ac4f54ee11aa5c342dfdba140958 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Wed, 1 Feb 2017 12:36:26 +0200
Subject: [PATCH] solvers.jl: Another way to solve Ax = b

Conflicts:
	src/solvers.jl
---
 src/solvers.jl | 42 +++++++++++++++++++++++++++++++++++++-----
 1 file changed, 37 insertions(+), 5 deletions(-)

diff --git a/src/solvers.jl b/src/solvers.jl
index d3a991c..918e504 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -233,15 +233,47 @@ end
 """
 Solve linear system using LU factorization (UMFPACK). This version solves
 directly the saddle point problem without elimination of boundary conditions.
+It is assumed that C1 == C2 and D = 0, so problem is symmetric and zero rows
+cand be removed from total system before solution. This kind of system arises
+in e.g. mesh tie problem
 """
 function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
-    nz = ones(solver.ndofs)
-    nz[get_nonzero_rows(C2)] = 0.0
-    nz[get_nonzero_rows(D)] = 0.0
-    D += spdiagm(nz)
+
+    C1 == C2 || return false
+    length(D) == 0 || return false
+
     A = [K C1'; C2  D]
     b = [f; g]
+
+    nz1 = get_nonzero_rows(A)
+    nz2 = get_nonzero_columns(A)
+    nz1 == nz2 || return false
+    
+    x = zeros(2*solver.ndofs)
+    x[nz1] = lufact(A[nz1,nz2]) \ full(b[nz1])
+
+    u[:] = x[1:solver.ndofs]
+    la[:] = x[solver.ndofs+1:end]
+
+    return true
+end
+
+"""
+Solve linear system using LU factorization (UMFPACK). This version solves
+directly the saddle point problem without elimination of boundary conditions.
+If matrix has zero rows, diagonal term is added to that matrix is invertible.
+"""
+function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{3}})
+
+    A = [K C1'; C2  D]
+    b = [f; g]
+
+    nz = ones(2*solver.ndofs)
+    nz[get_nonzero_rows(A)] = 0.0
+    A += spdiagm(nz)
+
     x = lufact(A) \ full(b)
+
     u[:] = x[1:solver.ndofs]
     la[:] = x[solver.ndofs+1:end]
     return true
@@ -281,7 +313,7 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
     la = zeros(ndofs)
     is_solved = false
     i = 0
-    for i in [1, 2]
+    for i in [1, 2, 3]
         is_solved = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
         if is_solved
             break