From 8fb9f3683b509595d54c01bbf95878de5e65d2ac Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 17:52:40 -0800 Subject: [PATCH] Easy way to create a padded tensor mesh. --- SimPEG/mesh/TensorMesh.py | 6 ++- SimPEG/utils/__init__.py | 4 +- SimPEG/utils/lomutils.py | 23 ----------- SimPEG/utils/meshutils.py | 58 +++++++++++++++++++++++++++ docs/api_Utils.rst | 7 ++++ docs/examples/mesh/plot_TensorMesh.py | 21 ---------- 6 files changed, 73 insertions(+), 46 deletions(-) create mode 100644 SimPEG/utils/meshutils.py delete mode 100644 docs/examples/mesh/plot_TensorMesh.py diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 4c0841b8..65f94b3e 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -24,7 +24,11 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): Example of a padded tensor mesh: - .. plot:: examples/mesh/plot_TensorMesh.py + .. plot:: + + from SimPEG import mesh, utils + M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() For a quick tensor mesh on a (10x12x15) unit cube:: diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 19aecec6..4bb1e445 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,9 +3,11 @@ import sputils import lomutils import interputils import ModelBuilder +import meshutils from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap -from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird +from meshutils import exampleLomGird, meshTensors +from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube from interputils import interpmat from ipythonUtils import easyAnimate as animate import Solver diff --git a/SimPEG/utils/lomutils.py b/SimPEG/utils/lomutils.py index 9020ce12..653cfdde 100644 --- a/SimPEG/utils/lomutils.py +++ b/SimPEG/utils/lomutils.py @@ -4,29 +4,6 @@ from matutils import mkvc, ndgrid, sub2ind from sputils import sdiag -def exampleLomGird(nC, exType): - assert type(nC) == list, "nC must be a list containing the number of nodes" - assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" - exType = exType.lower() - - possibleTypes = ['rect', 'rotate'] - assert exType in possibleTypes, "Not a possible example type." - - if exType == 'rect': - return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) - elif exType == 'rotate': - if len(nC) == 2: - X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) - amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2) - amt[amt < 0] = 0 - return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt - elif len(nC) == 3: - X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) - amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2) - amt[amt < 0] = 0 - return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt - - def volTetra(xyz, A, B, C, D): """ Returns the volume for tetrahedras volume specified by the indexes A to D. diff --git a/SimPEG/utils/meshutils.py b/SimPEG/utils/meshutils.py new file mode 100644 index 00000000..5c0a753c --- /dev/null +++ b/SimPEG/utils/meshutils.py @@ -0,0 +1,58 @@ +import numpy as np +from scipy import sparse as sp +from matutils import mkvc, ndgrid, sub2ind +from sputils import sdiag + +def exampleLomGird(nC, exType): + assert type(nC) == list, "nC must be a list containing the number of nodes" + assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" + exType = exType.lower() + + possibleTypes = ['rect', 'rotate'] + assert exType in possibleTypes, "Not a possible example type." + + if exType == 'rect': + return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + elif exType == 'rotate': + if len(nC) == 2: + X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2) + amt[amt < 0] = 0 + return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt + elif len(nC) == 3: + X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2) + amt[amt < 0] = 0 + return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt + + +def meshTensors(*args): + """ + **meshTensors** takes any number of tuples that have the form:: + + h1 = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCode), (numPad, sizeStart [, increaseFactor]) ) + + .. plot:: + + from SimPEG import mesh, utils + M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() + + """ + def padding(num, start, factor=1.3, reverse=False): + pad = ((np.ones(num)*factor)**np.arange(num))*start + if reverse: pad = pad[::-1] + return pad + tensors = tuple() + for i, arg in enumerate(args): + tensors += (np.r_[padding(*arg[0],reverse=True),np.ones(arg[1][0])*arg[1][1],padding(*arg[2])],) + + return list(tensors) if len(tensors) > 1 else tensors[0] + +if __name__ == '__main__': + from SimPEG import mesh + import matplotlib.pyplot as plt + M = mesh.TensorMesh(meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() + plt.gca().axis('tight') + plt.show() diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 07903e3c..7960b3e1 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -37,6 +37,13 @@ LOM Utilities :members: :undoc-members: +Mesh Utilities +************** + +.. automodule:: SimPEG.utils.meshutils + :members: + :undoc-members: + Model Builder Utilities *********************** diff --git a/docs/examples/mesh/plot_TensorMesh.py b/docs/examples/mesh/plot_TensorMesh.py deleted file mode 100644 index e773ed13..00000000 --- a/docs/examples/mesh/plot_TensorMesh.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from SimPEG.mesh import TensorMesh - -pad = 7 -padfactor = 1.4 -xpad = (np.ones(pad)*padfactor)**np.arange(pad) -ypad = (np.ones(pad)*padfactor)**np.arange(pad) - -core = 15 -xcore = np.ones(core) -ycore = np.ones(core) - -h1 = np.r_[xpad[::-1],xcore,xpad] -h2 = np.r_[ypad[::-1],ycore,ypad] - -mesh = TensorMesh([h1, h2]) -mesh.plotGrid() -plt.axis('tight') - -plt.show()