From 18fd590effcc776d7b1078b88d9bbcd18c87c98a Mon Sep 17 00:00:00 2001
From: SEOGI KANG <sgkang09@gmail.com>
Date: Wed, 10 Jul 2013 17:10:54 -0700
Subject: [PATCH] Generating 3D divergence matrix: getDIV.py Some funtions for
 sparse matrix: sputils.py -> Updated we've discussed before Test example for
 divergence

---
 code/getDIV.py         | 75 ++++++++++++++++++++++++++++++++++++++++++
 code/sputils.py        | 18 ++++++++++
 code/tests/test_div.py | 45 +++++++++++++++++++++++++
 3 files changed, 138 insertions(+)
 create mode 100644 code/getDIV.py
 create mode 100644 code/sputils.py
 create mode 100644 code/tests/test_div.py

diff --git a/code/getDIV.py b/code/getDIV.py
new file mode 100644
index 00000000..6c45c57c
--- /dev/null
+++ b/code/getDIV.py
@@ -0,0 +1,75 @@
+import numpy as np
+from scipy import sparse
+from utils import mkvc
+from sputils import ddx, sdiag, speye, kron3
+
+def getvol(h):   
+    
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # Compute cell volumes                        
+    v12 = h1.T*h2       
+    V = mkvc(v12.reshape(-1,1)*h3)
+    
+    return V
+        
+def getarea(h):   
+    
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)       
+    # Compute areas of cell faces
+    area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) 
+    area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
+    area3 = h1.T*mkvc(h2.T*np.ones(n3+1))               
+    area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
+    
+    return area
+
+def getDivMatrix(h):    
+    """Consturct the 3D divergence operator on Faces."""
+           
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)  
+        
+    # Compute areas of cell faces
+    #area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) 
+    #area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
+    #area3 = h1.T*mkvc(h2.T*np.ones(n3+1))               
+    #area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
+    area = getarea(h)
+            
+    S = sdiag(area)                   
+        
+    # Compute cell volumes                        
+    #v12 = h1.T*h2       
+    #V = mkvc(v12.reshape(-1,1)*h3)
+    V = getvol(h)
+        
+    # Compute divergence operator on faces                       
+    d1 = ddx(n1)
+    d2 = ddx(n2)
+    d3 = ddx(n3)
+    D1 = kron3(speye(n3), speye(n2), d1)
+    D2 = kron3(speye(n3), d2, speye(n1))
+    D3 = kron3(d3, speye(n2), speye(n1))       
+        
+    D = sparse.hstack((sparse.hstack((D1, D2)), D3))
+    return sdiag(1/V)*D*S
+
diff --git a/code/sputils.py b/code/sputils.py
new file mode 100644
index 00000000..204e71cc
--- /dev/null
+++ b/code/sputils.py
@@ -0,0 +1,18 @@
+import numpy as np
+from scipy import sparse
+
+def ddx(n):
+    """Define 1D derivatives"""
+    return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1)
+    
+def sdiag(h):
+    """Sparse diagonal matrix"""
+    return sparse.spdiags(h, 0, np.size(h), np.size(h))    
+
+def speye(n):
+    """Sparse identity"""
+    return sparse.identity(n)
+
+def kron3(A, B, C):
+    """Two kron prods"""
+    return sparse.kron(sparse.kron(A, B), C)  
\ No newline at end of file
diff --git a/code/tests/test_div.py b/code/tests/test_div.py
new file mode 100644
index 00000000..bc0c1a9d
--- /dev/null
+++ b/code/tests/test_div.py
@@ -0,0 +1,45 @@
+import numpy as np
+
+import sys
+sys.path.append('../')
+from TensorMesh import TensorMesh
+from getDIV import getDivMatrix, getarea, getvol
+
+# Define the mesh
+err=0.
+for i in range(4):
+    icount=i+1;
+    nc = 2*icount;
+    h1 = np.pi/nc*np.ones((1,nc))
+    h2 = np.pi/nc*np.ones((1,nc))
+    h3 = np.pi/nc*np.ones((1,nc))
+    h = [h1, h2, h3]
+    x0 = -np.pi/2*np.ones((3, 1))
+    M = TensorMesh(h, x0)
+    #n = M.plotGrid()
+
+    # Generate DIV matrix
+    DIV = getDivMatrix(h)
+    
+    #Test function
+    fun = lambda x: np.sin(x)
+    Fx = fun(M.gridFx[:,0])
+    Fy = fun(M.gridFy[:,1])
+    Fz = fun(M.gridFz[:,2])
+    
+    F = np.concatenate((Fx,Fy,Fz))
+    divF = DIV*F
+    sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
+    divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2])
+     
+    area = getarea(h)
+    vol = getvol(h)
+    err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
+    if icount == 1:
+        err1 = err
+        print 'h       |   2 norm   | error ratio'        
+        print '---------------------------------------'                
+        print '%6.4f  |  %8.2e  |'% (h1[0,0], err)
+    else:
+        print '%6.4f  |  %8.2e  |  %6.4f' % (h1[0,0], err, err1/err)
+