Add dirichlet node elimination

Author: j-fu

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "1.3.1"
+version = "1.4"
 
 [deps]
 AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"

docs/src/extsparse.md

@@ -13,6 +13,13 @@ LinearAlgebra.lu!
 LinearAlgebra.ldiv!
 ```
 
+## Handling of homogeneous Dirichlet BC
+```@docs
+mark_dirichlet
+eliminate_dirichlet!
+eliminate_dirichlet
+```
+
 ## Test matrix creation
 
 ```@autodocs

src/ExtendableSparse.jl

@@ -28,6 +28,8 @@ include("matrix/extendable.jl")
 export SparseMatrixLNK,
        ExtendableSparseMatrix, flush!, nnz, updateindex!, rawupdateindex!, colptrs, sparse
 
+export eliminate_dirichlet, eliminate_dirichlet!, mark_dirichlet
+
 include("factorizations/factorizations.jl")
 
 export JacobiPreconditioner,

src/matrix/extendable.jl

@@ -596,3 +596,20 @@ function pointblock(A0::ExtendableSparseMatrix{Tv,Ti},blocksize) where {Tv,Ti}
     flush!(Ab)
 end
 
+
+function mark_dirichlet(A::ExtendableSparseMatrix;penalty=1.0e20)
+    flush!(A)
+    mark_dirichlet(A.cscmatrix;penalty)
+end
+
+function eliminate_dirichlet(A::ExtendableSparseMatrix,dirichlet)
+    flush!(A)
+    ExtendableSparseMatrix(eliminate_dirichlet(A.cscmatrix,dirichlet))
+end
+
+function eliminate_dirichlet!(A::ExtendableSparseMatrix,dirichlet)
+    flush!(A)
+    eliminate_dirichlet!(A.cscmatrix,dirichlet)
+    A
+end
+

src/matrix/sparsematrixcsc.jl

@@ -88,3 +88,73 @@ end
 function pointblock(A::SparseMatrixCSC,blocksize)
     SparseMatrixCSC(pointblock(ExtendableSparseMatrix(A),blocksize))
 end
+
+"""
+   mark_dirichlet(A; penalty=1.0e20)
+
+Return boolean vector marking Dirichlet nodes, known by `A[i,i]>=penalty`
+"""
+function mark_dirichlet(A::SparseMatrixCSC; penalty=1.0e20)
+    (;colptr,rowval,nzval,n)=A
+    dirichlet=zeros(Bool,n)
+    for i=1:n
+        for j=colptr[i]:colptr[i+1]-1
+            if rowval[j]==i && nzval[j]>=penalty
+                dirichlet[i]=true
+            end
+        end
+    end
+    dirichlet
+end
+
+"""
+    eliminate_dirichlet!(A,dirichlet_marker)
+
+Eliminate dirichlet nodes in matrix by setting 
+```julia
+    A[:,i]=0; A[i,:]=0; A[i,i]=1
+```
+for a node `i` marked as Dirichlet.
+
+Returns A.
+"""
+function eliminate_dirichlet!(A::SparseMatrixCSC,dirichlet)
+    (;colptr,rowval,nzval,n)=A
+    for i=1:n
+        # set off-diagonal column indiced to zero
+        if !iszero(dirichlet[i])
+            for j=colptr[i]:colptr[i+1]-1
+                if rowval[j]==i
+                    nzval[j]=1
+                else
+                    nzval[j]=0
+                end
+            end
+        end
+        # set off-diagonal row indices to zero
+        for j=colptr[i]:colptr[i+1]-1
+            if rowval[j]!=i && !iszero(dirichlet[rowval[j]])
+                nzval[j]=0
+            end
+        end
+    end
+    A
+end
+
+"""
+    eliminate_dirichlet(A,dirichlet_marker)
+
+Create a copy B of A sharing the sparsity pattern.
+Eliminate dirichlet nodes in B by setting 
+```julia
+    B[:,i]=0; B[i,:]=0; B[i,i]=1
+```
+for a node `i` marked as Dirichlet.
+
+Returns B.
+"""
+function eliminate_dirichlet(A::SparseMatrixCSC,dirichlet)
+    (;m,n,colptr,rowval,nzval)=A
+    B=SparseMatrixCSC(m,n,colptr,rowval,copy(nzval))
+    eliminate_dirichlet!(B,dirichlet)
+end

test/runtests.jl

@@ -24,6 +24,8 @@ using ForwardDiff
 
 @testset "Backslash" begin include("test_backslash.jl") end
 
+@testset "Dirichlet" begin include("test_dirichlet.jl") end
+
 @testset "LinearSolve" begin include("test_linearsolve.jl") end
 
 @testset "Preconditioners" begin include("test_preconditioners.jl") end

test/test_dirichlet.jl

@@ -0,0 +1,30 @@
+module test_dirichlet
+using Test
+using ExtendableSparse
+using SparseArrays
+using LinearAlgebra
+
+function tdirichlet(A)
+    n=size(A,1)
+    for i=1:10:n
+        A[i,i]=1.0e30
+    end
+    f=ones(n)
+    u=A\f
+    diri=mark_dirichlet(A)
+    fD=f.* (1 .-diri)
+    AD=eliminate_dirichlet(A,diri)
+    uD=AD\fD
+    norm(uD-u,Inf) < 1.0e-20
+end
+
+@test tdirichlet(fdrand(1000; matrixtype = SparseMatrixCSC))
+@test tdirichlet(fdrand(1000; matrixtype = ExtendableSparseMatrix))
+
+@test tdirichlet(fdrand(20,20; matrixtype = SparseMatrixCSC))
+@test tdirichlet(fdrand(20,20; matrixtype = ExtendableSparseMatrix))
+
+@test tdirichlet(fdrand(10,10,10, matrixtype = SparseMatrixCSC))
+@test tdirichlet(fdrand(10,10,10; matrixtype = ExtendableSparseMatrix))
+
+end