From db694dd780331f16c9d7d000ce2babce1877f5dc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 18 Oct 2013 17:22:13 -0700 Subject: [PATCH 1/3] moved to mesh folder. --- SimPEG/TensorView.py | 2 +- SimPEG/__init__.py | 3 +- SimPEG/forward/DCProblem/DCProblem.py | 2 +- SimPEG/{ => mesh}/BaseMesh.py | 2 +- SimPEG/mesh/DiffOperators.py | 336 +++++++++++++++++++ SimPEG/{ => mesh}/InnerProducts.py | 0 SimPEG/{ => mesh}/LogicallyOrthogonalMesh.py | 2 +- SimPEG/{ => mesh}/LomView.py | 2 +- SimPEG/{ => mesh}/TensorMesh.py | 2 +- SimPEG/mesh/TensorView.py | 335 ++++++++++++++++++ SimPEG/mesh/__init__.py | 8 + SimPEG/tests/TestUtils.py | 3 +- SimPEG/tests/test_LogicallyOrthogonalMesh.py | 7 +- SimPEG/tests/test_basemesh.py | 3 +- SimPEG/tests/test_forward_DCproblem.py | 2 +- SimPEG/tests/test_forward_problem.py | 2 +- SimPEG/tests/test_tensorMesh.py | 2 +- SimPEG/utils/ModelBuilder.py | 3 +- docs/api_BaseMesh.rst | 2 +- docs/api_DiffOperators.rst | 2 +- docs/api_InnerProducts.rst | 2 +- docs/api_LOMView.rst | 2 +- docs/api_LogicallyOrthogonalMesh.rst | 2 +- docs/api_TensorMesh.rst | 2 +- docs/api_TensorView.rst | 2 +- 25 files changed, 703 insertions(+), 27 deletions(-) rename SimPEG/{ => mesh}/BaseMesh.py (99%) create mode 100644 SimPEG/mesh/DiffOperators.py rename SimPEG/{ => mesh}/InnerProducts.py (100%) rename SimPEG/{ => mesh}/LogicallyOrthogonalMesh.py (99%) rename SimPEG/{ => mesh}/LomView.py (99%) rename SimPEG/{ => mesh}/TensorMesh.py (96%) create mode 100644 SimPEG/mesh/TensorView.py create mode 100644 SimPEG/mesh/__init__.py diff --git a/SimPEG/TensorView.py b/SimPEG/TensorView.py index ba4f1b63..687a520d 100644 --- a/SimPEG/TensorView.py +++ b/SimPEG/TensorView.py @@ -2,7 +2,7 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib from mpl_toolkits.mplot3d import Axes3D -from utils import mkvc +from SimPEG.utils import mkvc class TensorView(object): diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index ea9faf33..e8946f88 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,5 +1,4 @@ -from TensorMesh import TensorMesh -from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh +import mesh import utils import inverse from Solver import Solver diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py index b8115598..1df8897b 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -1,4 +1,4 @@ -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh from SimPEG.forward import Problem, SyntheticProblem from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag diff --git a/SimPEG/BaseMesh.py b/SimPEG/mesh/BaseMesh.py similarity index 99% rename from SimPEG/BaseMesh.py rename to SimPEG/mesh/BaseMesh.py index c57828e9..29c43dd1 100644 --- a/SimPEG/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -1,5 +1,5 @@ import numpy as np -from utils import mkvc +from SimPEG.utils import mkvc class BaseMesh(object): diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py new file mode 100644 index 00000000..110fe9bd --- /dev/null +++ b/SimPEG/mesh/DiffOperators.py @@ -0,0 +1,336 @@ +import numpy as np +from scipy import sparse as sp +from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros + + +def ddx(n): + """Define 1D derivatives, inner, this means we go from n+1 to n+1""" + return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") + + +def checkBC(bc): + """ Checks if boundary condition 'bc' is valid. """ + if(type(bc) is str): + bc = [bc, bc] + assert type(bc) is list, 'bc must be a list' + assert len(bc) == 2, 'bc must have two elements' + + for bc_i in bc: + assert type(bc_i) is str, "each bc must be a string" + assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" + return bc + + +def ddxCellGrad(n, bc): + """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" + bc = checkBC(bc) + + D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") + # Set the first side + if(bc[0] == 'dirichlet'): + D[0, 0] = 2 + elif(bc[0] == 'neumann'): + D[0, 0] = 0 + # Set the second side + if(bc[1] == 'dirichlet'): + D[-1, -1] = -2 + elif(bc[1] == 'neumann'): + D[-1, -1] = 0 + return D + + +def av(n): + """Define 1D averaging operator from cell-centres to nodes.""" + return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") + + +class DiffOperators(object): + """ + Class creates the differential operators that you need! + """ + def __init__(self): + raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') + + def faceDiv(): + doc = "Construct divergence operator (face-stg to cell-centres)." + + def fget(self): + if(self._faceDiv is None): + # The number of cell centers in each direction + n = self.n + # Compute faceDivergence operator on faces + if(self.dim == 1): + D = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0])) + D = sp.hstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) + D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) + D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) + D = sp.hstack((D1, D2, D3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.vol + self._faceDiv = sdiag(1/V)*D*sdiag(S) + + return self._faceDiv + return locals() + _faceDiv = None + faceDiv = property(**faceDiv()) + + def nodalGrad(): + doc = "Construct gradient operator (nodes to edges)." + + def fget(self): + if(self._nodalGrad is None): + # The number of cell centers in each direction + n = self.n + # Compute divergence operator on faces + if(self.dim == 1): + G = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]+1), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) + G = sp.vstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) + D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) + D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) + G = sp.vstack((D1, D2, D3), format="csr") + # Compute lengths of cell edges + L = self.edge + self._nodalGrad = sdiag(1/L)*G + return self._nodalGrad + return locals() + _nodalGrad = None + nodalGrad = property(**nodalGrad()) + + def setCellGradBC(self, BC): + """ + Function that sets the boundary conditions for cell-centred derivative operators. + + Examples:: + + BC = 'neumann' # Neumann in all directions + BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else + BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain, + # Dirichlet else + + """ + if(type(BC) is str): + BC = [BC for _ in self.n] # Repeat the str self.dim times + elif(type(BC) is list): + assert len(BC) == self.dim, 'BC list must be the size of your mesh' + else: + raise Exception("BC must be a str or a list.") + + for i, bc_i in enumerate(BC): + BC[i] = checkBC(bc_i) + + self._cellGrad = None # ensure we create a new gradient next time we call it + self._cellGradBC = BC + return BC + _cellGradBC = 'neumann' + + def cellGrad(): + doc = "The cell centered Gradient, takes you to cell faces." + + def fget(self): + if(self._cellGrad is None): + BC = self.setCellGradBC(self._cellGradBC) + n = self.n + if(self.dim == 1): + G = ddxCellGrad(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) + G = sp.vstack((G1, G2), format="csr") + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) + G = sp.vstack((G1, G2, G3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.vol + self._cellGrad = sdiag(S)*G*sdiag(1/V) + return self._cellGrad + return locals() + _cellGrad = None + cellGrad = property(**cellGrad()) + + def edgeCurl(): + doc = "Construct the 3D curl operator." + + def fget(self): + if(self._edgeCurl is None): + # The number of cell centers in each direction + n1 = self.nCx + n2 = self.nCy + n3 = self.nCz + + # Compute lengths of cell edges + L = self.edge + + # Compute areas of cell faces + S = self.area + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + + 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]) + + C = sp.vstack((sp.hstack((O1, -D32, D23)), + sp.hstack((D31, O2, -D13)), + sp.hstack((-D21, D12, O3))), format="csr") + + self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) + return self._edgeCurl + return locals() + _edgeCurl = None + edgeCurl = property(**edgeCurl()) + + def aveF2CC(): + doc = "Construct the averaging operator on cell faces to cell centers." + + def fget(self): + if(self._aveF2CC is None): + n = self.n + if(self.dim == 1): + self._aveF2CC = av(n[0]) + elif(self.dim == 2): + self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])), + sp.kron(av(n[1]), speye(n[0]))), format="csr") + elif(self.dim == 3): + self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") + return self._aveF2CC + return locals() + _aveF2CC = None + aveF2CC = property(**aveF2CC()) + + def aveE2CC(): + doc = "Construct the averaging operator on cell edges to cell centers." + + def fget(self): + if(self._aveE2CC is None): + # The number of cell centers in each direction + n = self.n + if(self.dim == 1): + raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') + elif(self.dim == 2): + self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + sp.kron(speye(n[1]), av(n[0]))), format="csr") + elif(self.dim == 3): + self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") + return self._aveE2CC + return locals() + _aveE2CC = None + aveE2CC = property(**aveE2CC()) + + def aveN2CC(): + doc = "Construct the averaging operator on cell nodes to cell centers." + + def fget(self): + if(self._aveN2CC is None): + # The number of cell centers in each direction + n = self.n + if(self.dim == 1): + self._aveN2CC = av(n[0]) + elif(self.dim == 2): + self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])), + sp.kron(av(n[1]), av(n[0]))), format="csr") + elif(self.dim == 3): + self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") + return self._aveN2CC + return locals() + _aveN2CC = None + aveN2CC = property(**aveN2CC()) + + def aveN2CCv(): + doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." + + def fget(self): + if(self._aveN2CCv is None): + # The number of cell centers in each direction + n = self.n + if(self.dim == 1): + self._aveN2CCv = av(n[0]) + elif(self.dim == 2): + self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])), + sp.kron(av(n[1]), av(n[0]))), format="csr") + elif(self.dim == 3): + self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") + return self._aveN2CCv + return locals() + _aveN2CCv = None + aveN2CCv = property(**aveN2CCv()) + + def getMass(self, materialProp=None, loc='e'): + """ Produces mass matricies. + + :param str loc: Average to location: 'e'-edges, 'f'-faces + :param None,float,numpy.ndarray materialProp: property to be averaged (see below) + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the mass matrix + + materialProp can be:: + + None -> takes materialProp = 1 (default) + float -> a constant value for entire domain + numpy.ndarray -> if materialProp.size == self.nC + 3D property model + if materialProp.size = self.nCz + 1D (layered eath) property model + """ + if materialProp is None: + materialProp = np.ones(self.nC) + elif type(materialProp) is float: + materialProp = np.ones(self.nC)*materialProp + elif materialProp.shape == (self.nCz,): + materialProp = materialProp.repeat(self.nCx*self.nCy) + materialProp = mkvc(materialProp) + assert materialProp.shape == (self.nC,), "materialProp incorrect shape" + + if loc=='e': + Av = self.aveE2CC + elif loc=='f': + Av = self.aveF2CC + else: + raise ValueError('Invalid loc') + + diag = Av.T * (self.vol * mkvc(materialProp)) + + return sdiag(diag) + + def getEdgeMass(self, materialProp=None): + """mass matrix for products of edge functions w'*M(materialProp)*e""" + return self.getMass(loc='e', materialProp=materialProp) + + def getFaceMass(self, materialProp=None): + """mass matrix for products of face functions w'*M(materialProp)*f""" + return self.getMass(loc='f', materialProp=materialProp) + + def getFaceMassDeriv(self): + Av = self.aveF2CC + return Av.T * sdiag(self.vol) diff --git a/SimPEG/InnerProducts.py b/SimPEG/mesh/InnerProducts.py similarity index 100% rename from SimPEG/InnerProducts.py rename to SimPEG/mesh/InnerProducts.py diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py similarity index 99% rename from SimPEG/LogicallyOrthogonalMesh.py rename to SimPEG/mesh/LogicallyOrthogonalMesh.py index 5d2b7ffa..7b1e1bca 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -3,7 +3,7 @@ from BaseMesh import BaseMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts from LomView import LomView -from utils import mkvc, ndgrid, volTetra, indexCube, faceInfo +from SimPEG.utils import mkvc, ndgrid, volTetra, indexCube, faceInfo # Some helper functions. length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 diff --git a/SimPEG/LomView.py b/SimPEG/mesh/LomView.py similarity index 99% rename from SimPEG/LomView.py rename to SimPEG/mesh/LomView.py index 4b4f36a3..2a9b242b 100644 --- a/SimPEG/LomView.py +++ b/SimPEG/mesh/LomView.py @@ -2,7 +2,7 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib from mpl_toolkits.mplot3d import Axes3D -from utils import mkvc +from SimPEG.utils import mkvc class LomView(object): diff --git a/SimPEG/TensorMesh.py b/SimPEG/mesh/TensorMesh.py similarity index 96% rename from SimPEG/TensorMesh.py rename to SimPEG/mesh/TensorMesh.py index 2015838c..a3329c14 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -3,7 +3,7 @@ from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -from utils import ndgrid, mkvc +from SimPEG.utils import ndgrid, mkvc class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py new file mode 100644 index 00000000..687a520d --- /dev/null +++ b/SimPEG/mesh/TensorView.py @@ -0,0 +1,335 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +from mpl_toolkits.mplot3d import Axes3D +from SimPEG.utils import mkvc + + +class TensorView(object): + """ + Provides viewing functions for TensorMesh + + This class is inherited by TensorMesh + """ + def __init__(self): + pass + + def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): + """ + Mesh.plotImage(I) + + Plots scalar fields on the given mesh. + + Input: + + :param numpy.array I: scalar field + + Optional Input: + + :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz + :param int figNum: number of figure to plot to + :param matplotlib.axes.Axes ax: axis to plot to + :param str direction: slice dimensions, 3D only ('x', 'y', 'z') + :param bool numbering: show numbering of slices, 3D only + :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' + :param bool showIt: call plt.show() + + .. plot:: examples/mesh/plot_image_2D.py + :include-source: + + .. plot:: examples/mesh/plot_image_3D.py + :include-source: + """ + assert type(I) == np.ndarray, "I must be a numpy array" + assert type(numbering) == bool, "numbering must be a bool" + assert direction in ["x", "y","z"], "direction must be either x,y, or z" + + + if imageType == 'CC': + assert I.size == self.nC, "Incorrect dimensions for CC." + elif imageType == 'N': + assert I.size == self.nN, "Incorrect dimensions for N." + elif imageType == 'Fx': + if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') + elif imageType == 'Fy': + if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') + elif imageType == 'Fz': + if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') + elif imageType == 'Ex': + if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') + elif imageType == 'Ey': + if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') + elif imageType == 'Ez': + if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') + elif imageType[0] == 'E': + plotAll = len(imageType) == 1 + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + fig = plt.figure(figNum) + # Determine the subplot number: 131, 121 + numPlots = 130 if plotAll else len(imageType)/2*10+100 + pltNum = 1 + ex, ey, ez = self.r(I,'E','E','M') + if plotAll or 'Ex' in imageType: + ax_x = plt.subplot(numPlots+pltNum) + self.plotImage(ex, imageType='Ex', ax=ax_x, **options) + pltNum +=1 + if plotAll or 'Ey' in imageType: + ax_y = plt.subplot(numPlots+pltNum) + self.plotImage(ey, imageType='Ey', ax=ax_y, **options) + pltNum +=1 + if plotAll or 'Ez' in imageType: + ax_z = plt.subplot(numPlots+pltNum) + self.plotImage(ez, imageType='Ez', ax=ax_z, **options) + pltNum +=1 + return + elif imageType[0] == 'F': + plotAll = len(imageType) == 1 + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + fig = plt.figure(figNum) + # Determine the subplot number: 131, 121 + numPlots = 130 if plotAll else len(imageType)/2*10+100 + pltNum = 1 + fx, fy, fz = self.r(I,'F','F','M') + if plotAll or 'Fx' in imageType: + ax_x = plt.subplot(numPlots+pltNum) + self.plotImage(fx, imageType='Fx', ax=ax_x, **options) + pltNum +=1 + if plotAll or 'Fy' in imageType: + ax_y = plt.subplot(numPlots+pltNum) + self.plotImage(fy, imageType='Fy', ax=ax_y, **options) + pltNum +=1 + if plotAll or 'Fz' in imageType: + ax_z = plt.subplot(numPlots+pltNum) + self.plotImage(fz, imageType='Fz', ax=ax_z, **options) + pltNum +=1 + return + else: + raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") + + + if ax is None: + fig = plt.figure(figNum) + fig.clf() + ax = plt.subplot(111) + else: + assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" + fig = ax.figure + + if self.dim == 1: + if imageType == 'CC': + ph = ax.plot(self.vectorCCx, I, '-ro') + elif imageType == 'N': + ph = ax.plot(self.vectorNx, I, '-bs') + ax.set_xlabel("x") + ax.axis('tight') + elif self.dim == 2: + if imageType == 'CC': + C = I[:].reshape(self.n, order='F') + elif imageType == 'N': + C = I[:].reshape(self.n+1, order='F') + C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) + elif imageType == 'Fx': + C = I[:].reshape(self.nFx, order='F') + C = 0.5*(C[:-1, :] + C[1:, :] ) + elif imageType == 'Fy': + C = I[:].reshape(self.nFy, order='F') + C = 0.5*(C[:, :-1] + C[:, 1:] ) + elif imageType == 'Ex': + C = I[:].reshape(self.nEx, order='F') + C = 0.5*(C[:,:-1] + C[:,1:] ) + elif imageType == 'Ey': + C = I[:].reshape(self.nEy, order='F') + C = 0.5*(C[:-1,:] + C[1:,:] ) + + ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) + ax.axis('tight') + ax.set_xlabel("x") + ax.set_ylabel("y") + + elif self.dim == 3: + if direction == 'z': + + # get copy of image and average to cell-centres is necessary + if imageType == 'CC': + Ic = I[:].reshape(self.n, order='F') + elif imageType == 'N': + Ic = I[:].reshape(self.n+1, order='F') + Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) + elif imageType == 'Fx': + Ic = I[:].reshape(self.nFx, order='F') + Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) + elif imageType == 'Fy': + Ic = I[:].reshape(self.nFy, order='F') + Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) + elif imageType == 'Fz': + Ic = I[:].reshape(self.nFz, order='F') + Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) + elif imageType == 'Ex': + Ic = I[:].reshape(self.nEx, order='F') + Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) + elif imageType == 'Ey': + Ic = I[:].reshape(self.nEy, order='F') + Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) + elif imageType == 'Ez': + Ic = I[:].reshape(self.nEz, order='F') + Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) + + # determine number oE slices in x and y dimension + nX = np.ceil(np.sqrt(self.nCz)) + nY = np.ceil(self.nCz/nX) + + # allocate space for montage + nCx = self.nCx + nCy = self.nCy + + C = np.zeros((nX*nCx,nY*nCy)) + + for iy in range(int(nY)): + for ix in range(int(nX)): + iz = ix + iy*nX + if iz < self.nCz: + C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] + else: + C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan + + C = np.ma.masked_where(np.isnan(C), C) + xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] + yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] + # Plot the mesh + ph = ax.pcolormesh(xx, yy, C.T) + # Plot the lines + gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) + gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) + # Repeat and seperate with NaN + gxX = np.c_[gx, gx, gx+np.nan].ravel() + gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() + gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() + gyY = np.c_[gy, gy, gy+np.nan].ravel() + ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) + ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) + ax.axis('tight') + + if numbering: + pad = np.sum(self.hx)*0.04 + for iy in range(int(nY)): + for ix in range(int(nX)): + iz = ix + iy*nX + if iz < self.nCz: + ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, + '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') + + ax.set_title(imageType) + if showIt: plt.show() + return ph + + def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. + + :param bool nodes: plot nodes + :param bool faces: plot faces + :param bool centers: plot centers + :param bool edges: plot edges + :param bool lines: plot lines connecting nodes + :param bool showIt: call plt.show() + + .. plot:: examples/mesh/plot_grid_2D.py + :include-source: + + .. plot:: examples/mesh/plot_grid_3D.py + :include-source: + """ + if self.dim == 1: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + ax.hold(True) + ax.plot(xn, np.ones(np.shape(xn)), 'bs') + ax.plot(xc, np.ones(np.shape(xc)), 'ro') + ax.plot(xn, np.ones(np.shape(xn)), 'k--') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + if showIt: plt.show() + elif self.dim == 2: + fig = plt.figure(2) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + xs1 = self.gridFx + xs2 = self.gridFy + + ax.hold(True) + if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') + if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') + if faces: + ax.plot(xs1[:, 0], xs1[:, 1], 'g>') + ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + + # Plot the grid lines + if lines: + NN = self.r(self.gridN, 'N', 'N', 'M') + X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() + X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] + plt.plot(X, Y) + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + if showIt: plt.show() + elif self.dim == 3: + fig = plt.figure(3) + fig.clf() + ax = fig.add_subplot(111, projection='3d') + xn = self.gridN + xc = self.gridCC + xfs1 = self.gridFx + xfs2 = self.gridFy + xfs3 = self.gridFz + + xes1 = self.gridEx + xes2 = self.gridEy + xes3 = self.gridEz + + ax.hold(True) + if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) + if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) + if faces: + ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) + ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) + ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) + if edges: + ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) + ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) + ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) + + # Plot the grid lines + if lines: + NN = self.r(self.gridN, 'N', 'N', 'M') + X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() + Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() + X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() + Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() + X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() + Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() + Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() + X = np.r_[X1, X2, X3] + Y = np.r_[Y1, Y2, Y3] + Z = np.r_[Z1, Z2, Z3] + plt.plot(X, Y, 'b-', zs=Z) + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + ax.set_zlabel('x3') + if showIt: plt.show() diff --git a/SimPEG/mesh/__init__.py b/SimPEG/mesh/__init__.py new file mode 100644 index 00000000..9179efca --- /dev/null +++ b/SimPEG/mesh/__init__.py @@ -0,0 +1,8 @@ +from TensorMesh import TensorMesh +from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh +from BaseMesh import BaseMesh +from TensorView import TensorView +from LomView import LomView +from InnerProducts import InnerProducts +from DiffOperators import DiffOperators + diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 79340da7..9b2158c4 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -2,7 +2,8 @@ import numpy as np import matplotlib.pyplot as plt from pylab import norm from SimPEG.utils import mkvc -from SimPEG import TensorMesh, utils, LogicallyOrthogonalMesh +from SimPEG import utils +from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh import numpy as np import unittest diff --git a/SimPEG/tests/test_LogicallyOrthogonalMesh.py b/SimPEG/tests/test_LogicallyOrthogonalMesh.py index 523b18d2..7f656aae 100644 --- a/SimPEG/tests/test_LogicallyOrthogonalMesh.py +++ b/SimPEG/tests/test_LogicallyOrthogonalMesh.py @@ -1,10 +1,7 @@ import numpy as np import unittest -import sys -sys.path.append('../') -from TensorMesh import TensorMesh -from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh -from utils import ndgrid +from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh +from SimPEG.utils import ndgrid class BasicLOMTests(unittest.TestCase): diff --git a/SimPEG/tests/test_basemesh.py b/SimPEG/tests/test_basemesh.py index 60373011..5762de57 100644 --- a/SimPEG/tests/test_basemesh.py +++ b/SimPEG/tests/test_basemesh.py @@ -1,7 +1,6 @@ import unittest import sys -sys.path.append('../') -from BaseMesh import BaseMesh +from SimPEG.mesh import BaseMesh import numpy as np diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index b5e6d844..6ffa4650 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh from SimPEG.utils import ModelBuilder, sdiag from SimPEG.forward import Problem, SyntheticProblem from SimPEG.forward.DCProblem import DCProblem, DCutils diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py index b7180dca..d913a697 100644 --- a/SimPEG/tests/test_forward_problem.py +++ b/SimPEG/tests/test_forward_problem.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh from SimPEG.forward import Problem from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index 300e6634..9544dab6 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh from TestUtils import OrderTest from scipy.sparse.linalg import dsolve diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py index 170da1bf..d0512680 100644 --- a/SimPEG/utils/ModelBuilder.py +++ b/SimPEG/utils/ModelBuilder.py @@ -15,6 +15,7 @@ def getIndecesBlock(p0,p1,ccMesh): ccMesh represents the cell-centered mesh The points p0 and p1 must live in the the same dimensional space as the mesh. + """ # Validation: p0 and p1 live in the same dimensional space @@ -131,7 +132,7 @@ def scalarConductivity(ccMesh,pFunction): if __name__ == '__main__': - from SimPEG import TensorMesh + from SimPEG.mesh import TensorMesh from matplotlib import pyplot as plt # Define the mesh diff --git a/docs/api_BaseMesh.rst b/docs/api_BaseMesh.rst index 66249953..e5f81299 100644 --- a/docs/api_BaseMesh.rst +++ b/docs/api_BaseMesh.rst @@ -3,6 +3,6 @@ Base Mesh ********* -.. automodule:: SimPEG.BaseMesh +.. automodule:: SimPEG.mesh.BaseMesh :members: :undoc-members: diff --git a/docs/api_DiffOperators.rst b/docs/api_DiffOperators.rst index 57763a33..022fe564 100644 --- a/docs/api_DiffOperators.rst +++ b/docs/api_DiffOperators.rst @@ -3,6 +3,6 @@ Differential Operators ********************** -.. automodule:: SimPEG.DiffOperators +.. automodule:: SimPEG.mesh.DiffOperators :members: :undoc-members: diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index 6dfd2a77..90642807 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -3,6 +3,6 @@ Inner Products ************** -.. automodule:: SimPEG.InnerProducts +.. automodule:: SimPEG.mesh.InnerProducts :members: :undoc-members: diff --git a/docs/api_LOMView.rst b/docs/api_LOMView.rst index 61630c26..cfbfa6fa 100644 --- a/docs/api_LOMView.rst +++ b/docs/api_LOMView.rst @@ -3,6 +3,6 @@ LOM View ******** -.. automodule:: SimPEG.LomView +.. automodule:: SimPEG.mesh.LomView :members: :undoc-members: diff --git a/docs/api_LogicallyOrthogonalMesh.rst b/docs/api_LogicallyOrthogonalMesh.rst index 3b8296d7..af9e73d0 100644 --- a/docs/api_LogicallyOrthogonalMesh.rst +++ b/docs/api_LogicallyOrthogonalMesh.rst @@ -3,6 +3,6 @@ Logically Orthogonal Mesh ************************* -.. automodule:: SimPEG.LogicallyOrthogonalMesh +.. automodule:: SimPEG.mesh.LogicallyOrthogonalMesh :members: :undoc-members: diff --git a/docs/api_TensorMesh.rst b/docs/api_TensorMesh.rst index 9eb346cb..36fb925d 100644 --- a/docs/api_TensorMesh.rst +++ b/docs/api_TensorMesh.rst @@ -3,6 +3,6 @@ Tensor Mesh *********** -.. automodule:: SimPEG.TensorMesh +.. automodule:: SimPEG.mesh.TensorMesh :members: :undoc-members: diff --git a/docs/api_TensorView.rst b/docs/api_TensorView.rst index 6ae80276..c6b5a9b8 100644 --- a/docs/api_TensorView.rst +++ b/docs/api_TensorView.rst @@ -3,6 +3,6 @@ Tensor View *********** -.. automodule:: SimPEG.TensorView +.. automodule:: SimPEG.mesh.TensorView :members: :undoc-members: From 17d52087f8a537a745abf0cb200166eb51dff7fb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 18 Oct 2013 17:36:14 -0700 Subject: [PATCH 2/3] some updates for examples. --- SimPEG/DiffOperators.py | 336 -------------------------- SimPEG/TensorView.py | 335 ------------------------- docs/examples/mesh/plot_TensorMesh.py | 2 +- docs/examples/mesh/plot_grid_2D.py | 2 +- docs/examples/mesh/plot_grid_3D.py | 2 +- docs/examples/mesh/plot_image_2D.py | 2 +- docs/examples/mesh/plot_image_3D.py | 2 +- docs/examples/mesh/plot_nodes.py | 2 +- 8 files changed, 6 insertions(+), 677 deletions(-) delete mode 100644 SimPEG/DiffOperators.py delete mode 100644 SimPEG/TensorView.py diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py deleted file mode 100644 index 110fe9bd..00000000 --- a/SimPEG/DiffOperators.py +++ /dev/null @@ -1,336 +0,0 @@ -import numpy as np -from scipy import sparse as sp -from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros - - -def ddx(n): - """Define 1D derivatives, inner, this means we go from n+1 to n+1""" - return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") - - -def checkBC(bc): - """ Checks if boundary condition 'bc' is valid. """ - if(type(bc) is str): - bc = [bc, bc] - assert type(bc) is list, 'bc must be a list' - assert len(bc) == 2, 'bc must have two elements' - - for bc_i in bc: - assert type(bc_i) is str, "each bc must be a string" - assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" - return bc - - -def ddxCellGrad(n, bc): - """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" - bc = checkBC(bc) - - D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") - # Set the first side - if(bc[0] == 'dirichlet'): - D[0, 0] = 2 - elif(bc[0] == 'neumann'): - D[0, 0] = 0 - # Set the second side - if(bc[1] == 'dirichlet'): - D[-1, -1] = -2 - elif(bc[1] == 'neumann'): - D[-1, -1] = 0 - return D - - -def av(n): - """Define 1D averaging operator from cell-centres to nodes.""" - return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") - - -class DiffOperators(object): - """ - Class creates the differential operators that you need! - """ - def __init__(self): - raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') - - def faceDiv(): - doc = "Construct divergence operator (face-stg to cell-centres)." - - def fget(self): - if(self._faceDiv is None): - # The number of cell centers in each direction - n = self.n - # Compute faceDivergence operator on faces - if(self.dim == 1): - D = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]), ddx(n[0])) - D2 = sp.kron(ddx(n[1]), speye(n[0])) - D = sp.hstack((D1, D2), format="csr") - elif(self.dim == 3): - D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) - D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) - D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) - D = sp.hstack((D1, D2, D3), format="csr") - # Compute areas of cell faces & volumes - S = self.area - V = self.vol - self._faceDiv = sdiag(1/V)*D*sdiag(S) - - return self._faceDiv - return locals() - _faceDiv = None - faceDiv = property(**faceDiv()) - - def nodalGrad(): - doc = "Construct gradient operator (nodes to edges)." - - def fget(self): - if(self._nodalGrad is None): - # The number of cell centers in each direction - n = self.n - # Compute divergence operator on faces - if(self.dim == 1): - G = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]+1), ddx(n[0])) - D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) - G = sp.vstack((D1, D2), format="csr") - elif(self.dim == 3): - D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) - D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) - D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) - G = sp.vstack((D1, D2, D3), format="csr") - # Compute lengths of cell edges - L = self.edge - self._nodalGrad = sdiag(1/L)*G - return self._nodalGrad - return locals() - _nodalGrad = None - nodalGrad = property(**nodalGrad()) - - def setCellGradBC(self, BC): - """ - Function that sets the boundary conditions for cell-centred derivative operators. - - Examples:: - - BC = 'neumann' # Neumann in all directions - BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else - BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain, - # Dirichlet else - - """ - if(type(BC) is str): - BC = [BC for _ in self.n] # Repeat the str self.dim times - elif(type(BC) is list): - assert len(BC) == self.dim, 'BC list must be the size of your mesh' - else: - raise Exception("BC must be a str or a list.") - - for i, bc_i in enumerate(BC): - BC[i] = checkBC(bc_i) - - self._cellGrad = None # ensure we create a new gradient next time we call it - self._cellGradBC = BC - return BC - _cellGradBC = 'neumann' - - def cellGrad(): - doc = "The cell centered Gradient, takes you to cell faces." - - def fget(self): - if(self._cellGrad is None): - BC = self.setCellGradBC(self._cellGradBC) - n = self.n - if(self.dim == 1): - G = ddxCellGrad(n[0], BC[0]) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) - G = sp.vstack((G1, G2), format="csr") - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) - G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) - G = sp.vstack((G1, G2, G3), format="csr") - # Compute areas of cell faces & volumes - S = self.area - V = self.vol - self._cellGrad = sdiag(S)*G*sdiag(1/V) - return self._cellGrad - return locals() - _cellGrad = None - cellGrad = property(**cellGrad()) - - def edgeCurl(): - doc = "Construct the 3D curl operator." - - def fget(self): - if(self._edgeCurl is None): - # The number of cell centers in each direction - n1 = self.nCx - n2 = self.nCy - n3 = self.nCz - - # Compute lengths of cell edges - L = self.edge - - # Compute areas of cell faces - S = self.area - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - - 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]) - - C = sp.vstack((sp.hstack((O1, -D32, D23)), - sp.hstack((D31, O2, -D13)), - sp.hstack((-D21, D12, O3))), format="csr") - - self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) - return self._edgeCurl - return locals() - _edgeCurl = None - edgeCurl = property(**edgeCurl()) - - def aveF2CC(): - doc = "Construct the averaging operator on cell faces to cell centers." - - def fget(self): - if(self._aveF2CC is None): - n = self.n - if(self.dim == 1): - self._aveF2CC = av(n[0]) - elif(self.dim == 2): - self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") - elif(self.dim == 3): - self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") - return self._aveF2CC - return locals() - _aveF2CC = None - aveF2CC = property(**aveF2CC()) - - def aveE2CC(): - doc = "Construct the averaging operator on cell edges to cell centers." - - def fget(self): - if(self._aveE2CC is None): - # The number of cell centers in each direction - n = self.n - if(self.dim == 1): - raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') - elif(self.dim == 2): - self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveE2CC - return locals() - _aveE2CC = None - aveE2CC = property(**aveE2CC()) - - def aveN2CC(): - doc = "Construct the averaging operator on cell nodes to cell centers." - - def fget(self): - if(self._aveN2CC is None): - # The number of cell centers in each direction - n = self.n - if(self.dim == 1): - self._aveN2CC = av(n[0]) - elif(self.dim == 2): - self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveN2CC - return locals() - _aveN2CC = None - aveN2CC = property(**aveN2CC()) - - def aveN2CCv(): - doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." - - def fget(self): - if(self._aveN2CCv is None): - # The number of cell centers in each direction - n = self.n - if(self.dim == 1): - self._aveN2CCv = av(n[0]) - elif(self.dim == 2): - self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveN2CCv - return locals() - _aveN2CCv = None - aveN2CCv = property(**aveN2CCv()) - - def getMass(self, materialProp=None, loc='e'): - """ Produces mass matricies. - - :param str loc: Average to location: 'e'-edges, 'f'-faces - :param None,float,numpy.ndarray materialProp: property to be averaged (see below) - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the mass matrix - - materialProp can be:: - - None -> takes materialProp = 1 (default) - float -> a constant value for entire domain - numpy.ndarray -> if materialProp.size == self.nC - 3D property model - if materialProp.size = self.nCz - 1D (layered eath) property model - """ - if materialProp is None: - materialProp = np.ones(self.nC) - elif type(materialProp) is float: - materialProp = np.ones(self.nC)*materialProp - elif materialProp.shape == (self.nCz,): - materialProp = materialProp.repeat(self.nCx*self.nCy) - materialProp = mkvc(materialProp) - assert materialProp.shape == (self.nC,), "materialProp incorrect shape" - - if loc=='e': - Av = self.aveE2CC - elif loc=='f': - Av = self.aveF2CC - else: - raise ValueError('Invalid loc') - - diag = Av.T * (self.vol * mkvc(materialProp)) - - return sdiag(diag) - - def getEdgeMass(self, materialProp=None): - """mass matrix for products of edge functions w'*M(materialProp)*e""" - return self.getMass(loc='e', materialProp=materialProp) - - def getFaceMass(self, materialProp=None): - """mass matrix for products of face functions w'*M(materialProp)*f""" - return self.getMass(loc='f', materialProp=materialProp) - - def getFaceMassDeriv(self): - Av = self.aveF2CC - return Av.T * sdiag(self.vol) diff --git a/SimPEG/TensorView.py b/SimPEG/TensorView.py deleted file mode 100644 index 687a520d..00000000 --- a/SimPEG/TensorView.py +++ /dev/null @@ -1,335 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -from SimPEG.utils import mkvc - - -class TensorView(object): - """ - Provides viewing functions for TensorMesh - - This class is inherited by TensorMesh - """ - def __init__(self): - pass - - def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): - """ - Mesh.plotImage(I) - - Plots scalar fields on the given mesh. - - Input: - - :param numpy.array I: scalar field - - Optional Input: - - :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz - :param int figNum: number of figure to plot to - :param matplotlib.axes.Axes ax: axis to plot to - :param str direction: slice dimensions, 3D only ('x', 'y', 'z') - :param bool numbering: show numbering of slices, 3D only - :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_image_2D.py - :include-source: - - .. plot:: examples/mesh/plot_image_3D.py - :include-source: - """ - assert type(I) == np.ndarray, "I must be a numpy array" - assert type(numbering) == bool, "numbering must be a bool" - assert direction in ["x", "y","z"], "direction must be either x,y, or z" - - - if imageType == 'CC': - assert I.size == self.nC, "Incorrect dimensions for CC." - elif imageType == 'N': - assert I.size == self.nN, "Incorrect dimensions for N." - elif imageType == 'Fx': - if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') - elif imageType == 'Fy': - if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') - elif imageType == 'Fz': - if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') - elif imageType == 'Ex': - if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') - elif imageType == 'Ey': - if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') - elif imageType == 'Ez': - if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') - elif imageType[0] == 'E': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - ex, ey, ez = self.r(I,'E','E','M') - if plotAll or 'Ex' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(ex, imageType='Ex', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Ey' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(ey, imageType='Ey', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Ez' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(ez, imageType='Ez', ax=ax_z, **options) - pltNum +=1 - return - elif imageType[0] == 'F': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - fx, fy, fz = self.r(I,'F','F','M') - if plotAll or 'Fx' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(fx, imageType='Fx', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Fy' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(fy, imageType='Fy', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Fz' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(fz, imageType='Fz', ax=ax_z, **options) - pltNum +=1 - return - else: - raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") - - - if ax is None: - fig = plt.figure(figNum) - fig.clf() - ax = plt.subplot(111) - else: - assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" - fig = ax.figure - - if self.dim == 1: - if imageType == 'CC': - ph = ax.plot(self.vectorCCx, I, '-ro') - elif imageType == 'N': - ph = ax.plot(self.vectorNx, I, '-bs') - ax.set_xlabel("x") - ax.axis('tight') - elif self.dim == 2: - if imageType == 'CC': - C = I[:].reshape(self.n, order='F') - elif imageType == 'N': - C = I[:].reshape(self.n+1, order='F') - C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) - elif imageType == 'Fx': - C = I[:].reshape(self.nFx, order='F') - C = 0.5*(C[:-1, :] + C[1:, :] ) - elif imageType == 'Fy': - C = I[:].reshape(self.nFy, order='F') - C = 0.5*(C[:, :-1] + C[:, 1:] ) - elif imageType == 'Ex': - C = I[:].reshape(self.nEx, order='F') - C = 0.5*(C[:,:-1] + C[:,1:] ) - elif imageType == 'Ey': - C = I[:].reshape(self.nEy, order='F') - C = 0.5*(C[:-1,:] + C[1:,:] ) - - ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) - ax.axis('tight') - ax.set_xlabel("x") - ax.set_ylabel("y") - - elif self.dim == 3: - if direction == 'z': - - # get copy of image and average to cell-centres is necessary - if imageType == 'CC': - Ic = I[:].reshape(self.n, order='F') - elif imageType == 'N': - Ic = I[:].reshape(self.n+1, order='F') - Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) - elif imageType == 'Fx': - Ic = I[:].reshape(self.nFx, order='F') - Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) - elif imageType == 'Fy': - Ic = I[:].reshape(self.nFy, order='F') - Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) - elif imageType == 'Fz': - Ic = I[:].reshape(self.nFz, order='F') - Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) - elif imageType == 'Ex': - Ic = I[:].reshape(self.nEx, order='F') - Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) - elif imageType == 'Ey': - Ic = I[:].reshape(self.nEy, order='F') - Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) - elif imageType == 'Ez': - Ic = I[:].reshape(self.nEz, order='F') - Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) - - # determine number oE slices in x and y dimension - nX = np.ceil(np.sqrt(self.nCz)) - nY = np.ceil(self.nCz/nX) - - # allocate space for montage - nCx = self.nCx - nCy = self.nCy - - C = np.zeros((nX*nCx,nY*nCy)) - - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] - else: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan - - C = np.ma.masked_where(np.isnan(C), C) - xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] - yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] - # Plot the mesh - ph = ax.pcolormesh(xx, yy, C.T) - # Plot the lines - gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) - gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) - # Repeat and seperate with NaN - gxX = np.c_[gx, gx, gx+np.nan].ravel() - gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() - gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() - gyY = np.c_[gy, gy, gy+np.nan].ravel() - ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) - ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) - ax.axis('tight') - - if numbering: - pad = np.sum(self.hx)*0.04 - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, - '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') - - ax.set_title(imageType) - if showIt: plt.show() - return ph - - def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. - - :param bool nodes: plot nodes - :param bool faces: plot faces - :param bool centers: plot centers - :param bool edges: plot edges - :param bool lines: plot lines connecting nodes - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_grid_2D.py - :include-source: - - .. plot:: examples/mesh/plot_grid_3D.py - :include-source: - """ - if self.dim == 1: - fig = plt.figure(1) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - ax.hold(True) - ax.plot(xn, np.ones(np.shape(xn)), 'bs') - ax.plot(xc, np.ones(np.shape(xc)), 'ro') - ax.plot(xn, np.ones(np.shape(xn)), 'k--') - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - if showIt: plt.show() - elif self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - xs1 = self.gridFx - xs2 = self.gridFy - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') - if faces: - ax.plot(xs1[:, 0], xs1[:, 1], 'g>') - ax.plot(xs2[:, 0], xs2[:, 1], 'g^') - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] - plt.plot(X, Y) - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - if showIt: plt.show() - elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') - xn = self.gridN - xc = self.gridCC - xfs1 = self.gridFx - xfs2 = self.gridFy - xfs3 = self.gridFz - - xes1 = self.gridEx - xes2 = self.gridEy - xes3 = self.gridEz - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) - if faces: - ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) - ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) - ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) - if edges: - ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) - ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) - ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() - Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() - Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() - X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() - X = np.r_[X1, X2, X3] - Y = np.r_[Y1, Y2, Y3] - Z = np.r_[Z1, Z2, Z3] - plt.plot(X, Y, 'b-', zs=Z) - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - ax.set_zlabel('x3') - if showIt: plt.show() diff --git a/docs/examples/mesh/plot_TensorMesh.py b/docs/examples/mesh/plot_TensorMesh.py index 94c627f8..e773ed13 100644 --- a/docs/examples/mesh/plot_TensorMesh.py +++ b/docs/examples/mesh/plot_TensorMesh.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh pad = 7 padfactor = 1.4 diff --git a/docs/examples/mesh/plot_grid_2D.py b/docs/examples/mesh/plot_grid_2D.py index 258fbdfc..4c6d5da5 100644 --- a/docs/examples/mesh/plot_grid_2D.py +++ b/docs/examples/mesh/plot_grid_2D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh h1 = np.linspace(.1,.5,3) h2 = np.linspace(.1,.5,5) diff --git a/docs/examples/mesh/plot_grid_3D.py b/docs/examples/mesh/plot_grid_3D.py index caf8dbd2..6278802c 100644 --- a/docs/examples/mesh/plot_grid_3D.py +++ b/docs/examples/mesh/plot_grid_3D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh h1 = np.linspace(.1,.5,3) h2 = np.linspace(.1,.5,5) diff --git a/docs/examples/mesh/plot_image_2D.py b/docs/examples/mesh/plot_image_2D.py index c7bcb0d5..9fef3f1b 100644 --- a/docs/examples/mesh/plot_image_2D.py +++ b/docs/examples/mesh/plot_image_2D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh n = 20 h = np.ones(n)/n diff --git a/docs/examples/mesh/plot_image_3D.py b/docs/examples/mesh/plot_image_3D.py index d00756a1..c67b0d84 100644 --- a/docs/examples/mesh/plot_image_3D.py +++ b/docs/examples/mesh/plot_image_3D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh n = 20 h = np.ones(n)/n diff --git a/docs/examples/mesh/plot_nodes.py b/docs/examples/mesh/plot_nodes.py index 5d2c912b..ef11b351 100644 --- a/docs/examples/mesh/plot_nodes.py +++ b/docs/examples/mesh/plot_nodes.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh x0 = np.zeros(2) h1 = np.linspace(.1,.5,3) From ec42871640f6e508e158f90cd6c9dfdcebc9efc6 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 18 Oct 2013 17:43:11 -0700 Subject: [PATCH 3/3] Added 1D Cylindrical mesh. --- SimPEG/mesh/Cyl1DMesh.py | 330 +++++++++++++++++++++++++++++++++++++++ SimPEG/mesh/__init__.py | 2 +- docs/api_Cyl1DMesh.rst | 8 + docs/index.rst | 1 + 4 files changed, 340 insertions(+), 1 deletion(-) create mode 100644 SimPEG/mesh/Cyl1DMesh.py create mode 100644 docs/api_Cyl1DMesh.rst diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py new file mode 100644 index 00000000..5fc4deb6 --- /dev/null +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -0,0 +1,330 @@ +import numpy as np +import scipy.sparse as sp +from scipy.constants import pi +from SimPEG.utils import mkvc, ndgrid, sdiag + +class Cyl1DMesh(object): + """ + Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems + """ + + _meshType = 'CYL1D' + + def __init__(self, h, z0=None): + assert len(h) == 2, "len(h) must equal 2" + if z0 is not None: + assert z0.size == 1, "z0.size must equal 1" + + for i, h_i in enumerate(h): + assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) + assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) + + # Ensure h contains 1D vectors + self._h = [mkvc(x) for x in h] + + if z0 is None: + z0 = 0 + self._z0 = z0 + + #################################################### + # Mesh properties + #################################################### + + def h(): + doc = "list containing the width of each cell" + def fget(self): + return self._h + return locals() + h = property(**h()) + + def z0(): + doc = "The z-origin" + def fget(self): + return self._z0 + return locals() + z0 = property(**z0()) + + def hr(): + doc = "Width of the cells in the r direction" + def fget(self): + return self._h[0] + return locals() + hr = property(**hr()) + + def hz(): + doc = "Width of the cells in the z direction" + def fget(self): + return self._h[1] + return locals() + hz = property(**hz()) + + #################################################### + # Counting + #################################################### + + def nCr(): + doc = "Number of cells in the radial direction" + fget = lambda self: self.hr.size + return locals() + nCr = property(**nCr()) + + def nCz(): + doc = "Number of cells in the z direction" + fget = lambda self: self.hz.size + return locals() + nCz = property(**nCz()) + + def nC(): + doc = "Total number of cells" + fget = lambda self: self.nCr * self.nCz + return locals() + nC = property(**nC()) + + def nNr(): + doc = "Number of nodes in the radial direction" + fget = lambda self: self.hr.size + return locals() + nNr = property(**nNr()) + + def nNz(): + doc = "Number of nodes in the radial direction" + fget = lambda self: self.hz.size + 1 + return locals() + nNz = property(**nNz()) + + def nN(): + doc = "Total number of nodes" + fget = lambda self: self.nNr * self.nNz + return locals() + nN = property(**nN()) + + def nFr(): + doc = "Number of r faces" + fget = lambda self: self.nNr * self.nCz + return locals() + nFr = property(**nFr()) + + def nFz(): + doc = "Number of z faces" + fget = lambda self: self.nNz * self.nCr + return locals() + nFz = property(**nFz()) + + def nF(): + doc = "Total number of faces" + fget = lambda self: self.nFr + self.nFz + return locals() + nF = property(**nF()) + + def nE(): + doc = "Number of edges" + fget = lambda self: self.nN + return locals() + nE = property(**nE()) + + #################################################### + # Vectors & Grids + #################################################### + + def vectorNr(): + doc = "Nodal grid vector (1D) in the r direction" + fget = lambda self: self.hr.cumsum() + return locals() + vectorNr = property(**vectorNr()) + + def vectorNz(): + doc = "Nodal grid vector (1D) in the z direction" + fget = lambda self: np.r_[0, self.hz.cumsum()] + self._z0 + return locals() + vectorNz = property(**vectorNz()) + + def vectorCCr(): + doc = "Cell centered grid vector (1D) in the r direction" + fget = lambda self: np.r_[0, self.hr.cumsum()[1:] - self.hr[1:]/2] + return locals() + vectorCCr = property(**vectorCCr()) + + def vectorCCz(): + doc = "Cell centered grid vector (1D) in the z direction" + fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0 + return locals() + vectorCCz = property(**vectorCCz()) + + def gridCC(): + doc = "Cell-centered grid" + def fget(self): + if self._gridCC is None: + self._gridCC = ndgrid([self.vectorCCr, self.vectorCCz]) + return self._gridCC + return locals() + _gridCC = None + gridCC = property(**gridCC()) + + def gridN(): + doc = "Nodal grid" + def fget(self): + if self._gridN is None: + self._gridN = ndgrid([self.vectorNr, self.vectorNz]) + return self._gridN + return locals() + _gridN = None + gridN = property(**gridN()) + + def gridFr(): + doc = "r face grid" + def fget(self): + if self._gridFr is None: + self._gridFr = ndgrid([self.vectorNr, self.vectorCCz]) + return self._gridFr + return locals() + _gridFr = None + gridFr = property(**gridFr()) + + def gridFz(): + doc = "z face grid" + def fget(self): + if self._gridFz is None: + self._gridFz = ndgrid([self.vectorCCr, self.vectorNz]) + return self._gridFz + return locals() + _gridFz = None + gridFz = property(**gridFz()) + + #################################################### + # Geometries + #################################################### + + def edge(): + doc = "Edge lengths" + def fget(self): + if self._edge is None: + self._edge = 2*pi*self.gridN[:,0] + return self._edge + return locals() + _edge = None + edge = property(**edge()) + + def area(): + doc = "Face areas" + def fget(self): + if self._area is None: + areaR = np.kron(self.hz, 2*pi*self.vectorNr) + areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2)) + self._area = np.r_[areaR, areaZ] + return self._area + return locals() + _area = None + area = property(**area()) + + def vol(): + doc = "Volume of each cell" + def fget(self): + if self._vol is None: + az = pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2) + self._vol = np.kron(self.hz,az) + return self._vol + return locals() + _vol = None + vol = property(**vol()) + + #################################################### + # Operators + #################################################### + + def edgeCurl(): + doc = "The edgeCurl property." + def fget(self): + if self._edgeCurl is None: + #1D Difference matricies + dr = sp.spdiags((np.ones((self.nCr+1, 1))*[-1, 1]).T, [-1,0], self.nCr, self.nCr, format="csr") + dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0,1], self.nCz, self.nCz+1, format="csr") + + #2D Difference matricies + Dr = sp.kron(sp.eye(self.nNz), dr) + Dz = -sp.kron(dz, sp.eye(self.nCr)) #Not sure about this negative + + #Edge curl operator + self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0) + return self._edgeCurl + return locals() + _edgeCurl = None + edgeCurl = property(**edgeCurl()) + + def aveE2CC(): + doc = "Averaging operator from cell edges to cell centres" + def fget(self): + if self._aveE2CC is None: + az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar[0,0] = 1 + self._aveE2CC = sp.kron(az, ar).T + return self._aveE2CC + return locals() + _aveE2CC = None + aveE2CC = property(**aveE2CC()) + + def aveF2CC(): + doc = "Averaging operator from cell faces to cell centres" + def fget(self): + if self._aveF2CC is None: + az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar[0,0] = 1 + Afr = sp.kron(sp.eye(self.nCz),ar) + Afz = sp.kron(az,sp.eye(self.nCr)) + self._aveF2CC = sp.vstack((Afr,Afz)).T + return self._aveF2CC + return locals() + _aveF2CC = None + aveF2CC = property(**aveF2CC()) + + #################################################### + # Methods + #################################################### + + + def getMass(self, materialProp=None, loc='e'): + """ Produces mass matricies. + + :param str loc: Average to location: 'e'-edges, 'f'-faces + :param None,float,numpy.ndarray materialProp: property to be averaged (see below) + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the mass matrix + + materialProp can be:: + + None -> takes materialProp = 1 (default) + float -> a constant value for entire domain + numpy.ndarray -> if materialProp.size == self.nC + 3D property model + if materialProp.size = self.nCz + 1D (layered eath) property model + """ + if materialProp is None: + materialProp = np.ones(self.nC) + elif type(materialProp) is float: + materialProp = np.ones(self.nC)*materialProp + elif materialProp.shape == (self.nCz,): + materialProp = materialProp.repeat(self.nCr) + materialProp = mkvc(materialProp) + assert materialProp.shape == (self.nC,), "materialProp incorrect shape" + + if loc=='e': + Av = self.aveE2CC + elif loc=='f': + Av = self.aveF2CC + else: + raise ValueError('Invalid loc') + + diag = Av.T * (self.vol * mkvc(materialProp)) + + return sdiag(diag) + + def getEdgeMass(self, materialProp=None): + """mass matrix for products of edge functions w'*M(materialProp)*e""" + return self.getMass(loc='e', materialProp=materialProp) + + def getFaceMass(self, materialProp=None): + """mass matrix for products of face functions w'*M(materialProp)*f""" + return self.getMass(loc='f', materialProp=materialProp) + diff --git a/SimPEG/mesh/__init__.py b/SimPEG/mesh/__init__.py index 9179efca..3b8e1eef 100644 --- a/SimPEG/mesh/__init__.py +++ b/SimPEG/mesh/__init__.py @@ -1,3 +1,4 @@ +from Cyl1DMesh import Cyl1DMesh from TensorMesh import TensorMesh from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh from BaseMesh import BaseMesh @@ -5,4 +6,3 @@ from TensorView import TensorView from LomView import LomView from InnerProducts import InnerProducts from DiffOperators import DiffOperators - diff --git a/docs/api_Cyl1DMesh.rst b/docs/api_Cyl1DMesh.rst new file mode 100644 index 00000000..f0109b0d --- /dev/null +++ b/docs/api_Cyl1DMesh.rst @@ -0,0 +1,8 @@ +.. _api_Cyl1DMesh: + +Cylindrical 1D Mesh +******************* + +.. automodule:: SimPEG.mesh.Cyl1DMesh + :members: + :undoc-members: \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index 4c01c4cd..692363a2 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -26,6 +26,7 @@ Meshing & Operators api_TensorMesh api_TensorView api_LogicallyOrthogonalMesh + api_Cyl1DMesh api_LOMView api_DiffOperators api_InnerProducts