From f5c950540f3035afb6da3d7bbd77ad81833225c4 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Tue, 23 Jun 2015 23:22:39 +0300
Subject: [PATCH] added sparse assembly

---
 src/elasticity_solver.jl       | 82 ++++++++++++++++++++++++++++++++++
 test/test_elasticity_solver.jl | 75 +++++++++++++++++++++++++++++--
 2 files changed, 153 insertions(+), 4 deletions(-)

diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl
index 406a527..3e4c0e9 100644
--- a/src/elasticity_solver.jl
+++ b/src/elasticity_solver.jl
@@ -1,6 +1,9 @@
 # This file is a part of JuliaFEM. License is MIT: https://github.com/ovainola/JuliaFEM/blob/master/README.md
 module elasticity_solver
 
+using Logging
+@Logging.configure(level=DEBUG)
+
 VERSION < v"0.4-" && using Docile
 
 export solve_elasticity_interface!
@@ -69,6 +72,85 @@ function calc_local_matrices!(X, u, R, Kt, N, dNdξ, λ_, μ_, ipoints, iweights
 end
 
 
+@doc """
+Assemble global stiffness matrix to I,J,V ready for sparse format
+
+Parameters
+----------
+ke : local matrix
+eldofs_ : Array
+  degrees of freedom
+I,J,V : Arrays for sparse matrix
+
+Notes
+-----
+eldofs can also be node ids for convenience. In that case dimension
+is calculated and eldofs are "extended" to problem dimension.
+""" ->
+function assemble!(ke, eldofs_, I, J, V)
+    n, m = size(ke)
+    dim = round(Int, n/length(eldofs_))
+    @debug("problem dim = ", dim)
+    if dim == 1
+        eldofs = eldofs_
+    else
+        eldofs = Int64[]
+        for i in eldofs_
+            for d in 1:dim
+                push!(eldofs, dim*(i-1)+d)
+            end
+        end
+        @debug("old eldofs", eldofs_)
+        @debug("new eldofs", eldofs)
+    end
+    for i in 1:n
+        for j in 1:m
+            push!(I, eldofs[i])
+            push!(J, eldofs[j])
+            push!(V, ke[i,j])
+        end
+    end
+end
+
+
+@doc """
+Assemble global RHS to I,V ready for sparse format
+
+Parameters
+----------
+fe : local vector
+eldofs_ : Array
+  degrees of freedom
+I,V : Arrays for sparse matrix
+
+Notes
+-----
+eldofs can also be node ids for convenience. In that case dimension
+is calculated and eldofs are "extended" to problem dimension.
+""" ->
+function assemble!(fe, eldofs_, I, V)
+    n = length(fe)
+    dim = round(Int, n/length(eldofs_))
+    if dim == 1
+        eldofs = eldofs_
+    else
+        eldofs = Int64[]
+        for i in eldofs_
+            for d in 1:dim
+                push!(eldofs, dim*(i-1)+d)
+            end
+        end
+        @debug("old eldofs", eldofs_)
+        @debug("new eldofs", eldofs)
+    end
+
+    for i in 1:n
+        push!(I, eldofs[i])
+        push!(V, fe[i])
+    end
+end
+
+
 @doc """
 Solve one increment of elasticity problem
 """ ->
diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index 3c5564b..9c508ec 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -1,11 +1,11 @@
-using JuliaFEM.elasticity_solver
-
 using FactCheck
 using Logging
 @Logging.configure(level=DEBUG)
 
 
+using JuliaFEM.elasticity_solver: solve_elasticity_increment!
 facts("test solve elasticity increment") do
+
   X = [0 0; 10 0; 10 1; 0 1]'
   elmap = [1; 2; 3; 4]
   nodalloads = [0 0; 0 0; 0 -2; 0 0]'
@@ -41,8 +41,8 @@ facts("test solve elasticity increment") do
   iweights = [1 1 1 1]
 
   for i=1:10
-    JuliaFEM.elasticity_solver.solve_elasticity_increment!(
-      X, u, du, R, Kt,elmap, nodalloads, dirichletbc, la, mu, N, dNdξ, ipoints, iweights)
+    solve_elasticity_increment!(X, u, du, R, Kt,elmap, nodalloads, dirichletbc,
+                                la, mu, N, dNdξ, ipoints, iweights)
     @debug("increment:\n",du)
     u += du
     if norm(du) < 1.0e-9
@@ -52,3 +52,70 @@ facts("test solve elasticity increment") do
   @debug("solution\n",u)
   @fact u[2, 3] => roughly(-2.222244754401764)  # Tested against Elmer solution
 end
+
+
+using JuliaFEM.elasticity_solver: assemble!
+
+facts("test assembly of global matrix for 1 dim/node case") do
+    I = Int64[]
+    J = Int64[]
+    V = Float64[]
+    ke = [3 1; 1 1]
+    eldofs = [1, 2]
+    assemble!(ke, eldofs, I, J, V)
+    ke = [4 2; 2 3]
+    eldofs = [2, 4]
+    assemble!(ke, eldofs, I, J, V)
+    S = full(sparse(I, J, V))
+    @fact S => [3.0 1.0 0.0 0.0
+                1.0 5.0 0.0 2.0
+                0.0 0.0 0.0 0.0
+                0.0 2.0 0.0 3.0]
+end
+
+facts("test assembly of global vector for 1 dim/node case") do
+    I = Int64[]
+    V = Float64[]
+    fe = [1;2]
+    eldofs = [1, 2]
+    assemble!(fe, eldofs, I, V)
+    fe = [3;1]
+    eldofs = [2, 4]
+    assemble!(fe, eldofs, I, V)
+    S = full(sparsevec(I, V))
+    @fact S => [1.0 5.0 0.0 1.0]'
+end
+
+facts("test assembly of global matrix for 2 dim/node case") do
+    # provide "convienence" function, if given only nodal connectivity
+    # automatically find out dimension and "extend" matrix to full
+    I = Int64[]
+    J = Int64[]
+    V = Float64[]
+    ke = reshape(1:16, 4, 4)
+    eldofs = [1, 2]
+    assemble!(ke, eldofs, I, J, V)
+    eldofs = [2, 3]
+    assemble!(2*ke, eldofs, I, J, V)
+    expected = zeros(6, 6)
+    expected[1:4,1:4] += ke
+    expected[3:6,3:6] += 2*ke
+    S = full(sparse(I, J, V))
+    @fact S => expected
+end
+
+facts("test assembly of global vector for 2 dim/node case") do
+    # provide "convienence" function, if given only nodal connectivity
+    # automatically find out dimension and "extend" matrix to full
+    I = Int64[]
+    J = Int64[]
+    V = Float64[]
+    fe = [1, 2, 3, 4]
+    eldofs = [1, 2]
+    assemble!(fe, eldofs, I, V)
+    eldofs = [2, 3]
+    assemble!(2*fe, eldofs, I, V)
+    S = full(sparsevec(I, V))
+    expected = [1.0 2.0 5.0 8.0 6.0 8.0]'
+    @fact S => expected
+end