From 1e68873aadf1bdc7cdfd6996346f68711860081d Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Wed, 17 Jul 2013 23:47:31 -0700 Subject: [PATCH 1/2] Added the curl and grad and changed the area,vol such that they can work with h as defined by the grid --- SimPEG/getDiffOpps.py | 150 ++++++++++++++++++++++++++++++++++++++++++ SimPEG/sputils.py | 60 ++++++++++++++++- 2 files changed, 209 insertions(+), 1 deletion(-) create mode 100644 SimPEG/getDiffOpps.py diff --git a/SimPEG/getDiffOpps.py b/SimPEG/getDiffOpps.py new file mode 100644 index 00000000..c8391ab8 --- /dev/null +++ b/SimPEG/getDiffOpps.py @@ -0,0 +1,150 @@ +import numpy as np +from scipy import sparse +from utils import mkvc +from sputils import * +#from sputils import ddx, sdiag, speye, kron3, spzeros, appendBottom3, + +def getvol(h): + + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # Compute cell volumes + V = mkvc(np.outer(mkvc(np.outer(h1,h2)),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 = mkvc(np.outer(np.ones(n1+1),np.outer(h2,h3))) + area2 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),h3)))) + area3 = mkvc(np.outer(h1,mkvc(np.outer(h2,np.ones(n3+1))))) + area = np.hstack((np.hstack((area1, area2)), area3)) + + return area + +def getLength(h): + + 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 the length of each edge + Length1 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),np.ones(n3+1))))) + Length2 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(h2,np.ones(n3+1))))) + Length3 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(np.ones(n2+1),h3)))) + + Length = np.hstack((np.hstack((Length1, Length2)), Length3)) + + return Length + + +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) + + area = getarea(h) + S = sdiag(area) + + # Compute cell volumes + 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 + + +def getCurlMatrix(h): + """Edge CURL """ + + # 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) + + d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3) + # derivatives on x-edge variables + D32 = kron3(d3, speye(n2), speye(n1+1)) + D23 = kron3(speye(n3), d2, speye(n1+1)) + D31 = kron3(d3, speye(n2+1), speye(n1)) + D13 = kron3(speye(n3), speye(n2+1), d1) + D21 = kron3(speye(n3+1), d2, speye(n1)) + D12 = kron3(speye(n3+1), speye(n2), d1) + + O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) + O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) + O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) + + CURL = appendBottom3( + appendRight3(O1, -D32, D23), + appendRight3(D31, O2, -D13), + appendRight3(-D21, D12, O3)) + + + area = getarea(h) + S = sdiag(1/area) + + # Compute edge length + lngth = getLength(h) + L = sdiag(lngth) + + return S*(CURL*L) + + +def getNodalGradient(h): + """Nodal Gradients""" + + # 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) + + + D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1)) + D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1)) + D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1)) + + # topological gradient + GRAD = appendBottom3(D1, D2, D3) + + # scale for non-uniform mesh + # Compute edge length + lngth = getLength(h) + L = sdiag(1/lngth) + + return L*GRAD diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index 204e71cc..7371e46b 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -1,5 +1,7 @@ import numpy as np from scipy import sparse +from numpy import ones + def ddx(n): """Define 1D derivatives""" @@ -15,4 +17,60 @@ def speye(n): def kron3(A, B, C): """Two kron prods""" - return sparse.kron(sparse.kron(A, B), C) \ No newline at end of file + return sparse.kron(sparse.kron(A, B), C) + +def av(n): + """Define 1D average""" + return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)) + +def spzeros(n1, n2): + """spzeros""" + return sparse.coo_matrix((n1, n2)) + +def appendBottom(A, B): + """append on bottom""" + C = sparse.vstack((A, B)) + C = C.tocsr() + return C + + +def appendBottom3(A, B, C): + """append on bottom""" + C = appendBottom(appendBottom(A, B), C) + C = C.tocsr() + return C + + +def appendRight(A, B): + """append on right""" + C = sparse.hstack((A, B)) + C = C.tocsr() + return C + + +def appendRight3(A, B, C): + """append on right""" + C = appendRight(appendRight(A, B), C) + C = C.tocsr() + return C + + +def blkDiag(A, B): + """blockdigonal""" + O12 = sparse.coo_matrix((np.shape(A)[0], np.shape(B)[1])) + O21 = sparse.coo_matrix((np.shape(B)[0], np.shape(A)[1])) + C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B)))) + C = C.tocsr() + return C + + +def blkDiag3(A, B, C): + """blockdigonal 3""" + ABC = blkDiag(blkDiag(A, B), C) + ABC = ABC.tocsr() + return ABC + + + + + \ No newline at end of file From 12b84d14ee9aefddaf72487de5a5ebd743dc7ba1 Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Thu, 18 Jul 2013 00:09:48 -0700 Subject: [PATCH 2/2] A few more changes - added mass matrices as well as name conventions --- SimPEG/getDiffOpps.py | 67 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 62 insertions(+), 5 deletions(-) diff --git a/SimPEG/getDiffOpps.py b/SimPEG/getDiffOpps.py index c8391ab8..bcfe2c8a 100644 --- a/SimPEG/getDiffOpps.py +++ b/SimPEG/getDiffOpps.py @@ -4,7 +4,7 @@ from utils import mkvc from sputils import * #from sputils import ddx, sdiag, speye, kron3, spzeros, appendBottom3, -def getvol(h): +def getVol(h): # Cell sizes in each direction h1 = h[0] @@ -16,7 +16,7 @@ def getvol(h): return V -def getarea(h): +def getArea(h): # Cell sizes in each direction h1 = h[0] @@ -69,11 +69,11 @@ def getDivMatrix(h): n2 = np.size(h2) n3 = np.size(h3) - area = getarea(h) + area = getArea(h) S = sdiag(area) # Compute cell volumes - V = getvol(h) + V = getVol(h) # Compute divergence operator on faces d1 = ddx(n1) @@ -115,7 +115,7 @@ def getCurlMatrix(h): appendRight3(-D21, D12, O3)) - area = getarea(h) + area = getArea(h) S = sdiag(1/area) # Compute edge length @@ -148,3 +148,60 @@ def getNodalGradient(h): L = sdiag(1/lngth) return L*GRAD + + +def getEdgeToCellAverge(h): + + """Average from Edge to Cell center """ + + # 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) + + a1 = av(n1); a2 = av(n2); a3 = av(n3) + # derivatives on x-edge variables + A1 = kron3(a3, a2, speye(n1)) + A2 = kron3(a3, speye(n2), a1) + A3 = kron3(speye(n3), a2, a1) + + return appendRight3(A1, A2, A3) + + +def getFaceToCellAverge(h): + + """Average from Edge to Cell center """ + + # 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) + + a1 = av(n1); a2 = av(n2); a3 = av(n3) + # derivatives on x-edge variables + A1 = kron3(speye(n3), speye(n2), a1) + A2 = kron3(speye(n3), a2, speye(n1)) + A3 = kron3(a3, speye(n2), speye(n1)) + + return appendRight3(A1, A2, A3) + + +def getEdgeMassMatrix(h,sigma): + # mass matix for products of edge functions w'*M(sigma)*e + + Av = getEdgeToCellAverge(h) + v = getVol(h) + sigma = mkvc(sigma) + + return sdiag(Av.T*(v*sigma)) + +def getFaceMassMatrix(h,sigma): + # mass matix for products of edge functions w'*M(sigma)*e + + Av = getFaceToCellAverge(h) + v = getVol(h) + sigma = mkvc(sigma) + + return sdiag(Av.T*(v*sigma)) \ No newline at end of file