wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-01 21:32:22 +08:00
Code Issues Packages Projects Releases Wiki Activity

Bug Fixes in TensorView and test_utils

Browse Source 
Updated ipnbs
This commit is contained in:
Rowan Cockett
2013-08-30 21:34:15 -07:00
parent 9392747d1b
commit e11a3ccf65
5 changed files with 184 additions and 150 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/TensorView.py
+8 -7
View File
@@ -13,7 +13,7 @@ class TensorView(object):
def __init__(self):
pass
def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True):
def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w'):
"""
Mesh.plotImage(I)
@@ -126,11 +126,11 @@ class TensorView(object):
nY = np.ceil(self.nCz/nX)
# allocate space for montage
C = np.zeros((nX*self.nCx,nY*self.nCz))
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
@@ -142,6 +142,7 @@ class TensorView(object):
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]
@@ -151,8 +152,8 @@ class TensorView(object):
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, 'w-', linewidth=2)
ax.plot(gyX, gyY, 'w-', linewidth=2)
ax.plot(gxX, gxY, annotationColor+'-', linewidth=2)
ax.plot(gyX, gyY, annotationColor+'-', linewidth=2)
if numbering:
pad = np.sum(self.hx)*0.04
@@ -161,7 +162,7 @@ class TensorView(object):
iz = ix + iy*nX
if iz < self.nCz:
ax.text((ix+1)*self.vectorNx[-1]-pad,(iy)*self.vectorNy[-1]+pad,
'#%i'%iz,color='w',verticalalignment='bottom',horizontalalignment='right',size='x-large')
'#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
plt.show()
return ph
SimPEG/parameters.py
+64 -64
View File
@@ -5,64 +5,64 @@ import TensorView as tv
def getIndecesBlock(p0,p1,ccMesh):
"""
Creates a vector containing the block indexes in the cell centerd mesh.
Returns a tuple
Creates a vector containing the block indexes in the cell centerd mesh.
Returns a tuple
The block is defined by the points
p0 : describe the position of the left upper front corner, and
p0 : describe the position of the left upper front corner, and
p1 : describe the position of the right bottom back corner.
ccMesh represents the cell-centered mesh
The points p0 and p1 must live in the the same dimensional space as the mesh.
The points p0 and p1 must live in the the same dimensional space as the mesh.
"""
# Validation of the input
assert type(p0) == np.ndarray, "Vector must be a numpy array"
assert type(p1) == np.ndarray, "Vector must be a numpy array"
# Validation: p0 and p1 live in the same dimensional space
assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)"
# Validation: mesh and points live in the same dimensional space
dimMesh = np.size(ccMesh[0,:])
assert len(p0) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
if dimMesh == 1:
# Define the reference points
x1 = p0[0]
x2 = p1[0]
x1 = p0[0]
x2 = p1[0]
indX = (x1 <= ccMesh[:,0]) & (ccMesh[:,0] <= x2)
ind = np.where(indX)
elif dimMesh == 2:
# Define the reference points
x1 = p0[0]
y1 = p0[1]
x2 = p1[0]
y2 = p1[1]
x1 = p0[0]
y1 = p0[1]
x2 = p1[0]
y2 = p1[1]
indX = (x1 <= ccMesh[:,0]) & (ccMesh[:,0] <= x2)
indY = (y1 <= ccMesh[:,1]) & (ccMesh[:,1] <= y2)
ind = np.where(indX & indY)
else:
# Define the points
x1 = p0[0]
y1 = p0[1]
x1 = p0[0]
y1 = p0[1]
z1 = p0[2]
x2 = p1[0]
y2 = p1[1]
x2 = p1[0]
y2 = p1[1]
z2 = p1[2]
indX = (x1 <= ccMesh[:,0]) & (ccMesh[:,0] <= x2)
indY = (y1 <= ccMesh[:,1]) & (ccMesh[:,1] <= y2)
indZ = (z1 <= ccMesh[:,2]) & (ccMesh[:,2] <= z2)
ind = np.where(indX & indY & indZ)
# Return a tuple
@@ -76,28 +76,28 @@ def defineBlockConductivity(p0,p1,ccMesh,condVals):
"""
sigma = np.zeros(ccMesh.shape[0]) + condVals[1]
ind = getIndecesBlock(p0,p1,ccMesh)
sigma[ind] = condVals[0]
return sigma
def defineTwoLayeredConductivity(depth,ccMesh,condVals):
"""
Define a two layered model. Depth of the first layer must be specified.
CondVals vector with the conductivity values of the layers. Eg:
Convention to number the layers:
Convention to number the layers:
<----------------------------|------------------------------------>
0 depth zf
1st layer 2nd layer
"""
sigma = np.zeros(ccMesh.shape[0]) + condVals[1]
dim = np.size(ccMesh[0,:])
p0 = np.zeros(dim)
p1 = np.zeros(dim)
# Identify 1st cell centered reference point
p0[0] = ccMesh[0,0]
p0[1] = ccMesh[0,1]
@@ -107,30 +107,30 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
p1[0] = ccMesh[-1,0]
p1[1] = ccMesh[-1,1]
p1[2] = ccMesh[-1,2] - depth;
ind = getIndecesBlock(p0,p1,ccMesh)
sigma[ind] = condVals[0];
return sigma
def scalarConductivity(ccMesh,pFunction):
"""
Define the distribution conductivity in the mesh according to the
analytical expression given in pFunction
Define the distribution conductivity in the mesh according to the
analytical expression given in pFunction
"""
xCC = ccMesh[:,0]
yCC = ccMesh[:,1]
zCC = ccMesh[:,2]
sigma = pFunction(xCC,yCC,zCC)
return sigma
if __name__ == '__main__':
# Define the mesh
testDim = 3
h1 = 0.3*np.ones(7)
h1[0] = 0.5
@@ -149,44 +149,44 @@ if __name__ == '__main__':
h = [h1, h2, h3]
M = tm.TensorMesh(h, x0)
ccMesh = M.gridCC
# ------------------- Test conductivities! --------------------------
print('Testing 1 block conductivity')
p0 = np.array([0.5,0.5,0.5])
p1 = np.array([1.0,1.0,1.0])
condVals = np.array([100,1e-6])
sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
# Plot sigma model
#M.plotImage(sigma)
print sigma
print sigma.shape
M.plotImage(sigma)
print 'Done with block! :)'
# -----------------------------------------
print('Testing the two layered model')
print('Testing the two layered model')
condVals = np.array([100,1e-5]);
depth = 1.0;
sigma = defineTwoLayeredConductivity(depth,ccMesh,condVals)
#M.plotImage(sigma)
M.plotImage(sigma)
print sigma
print 'layer model!'
# -----------------------------------------
print('Testing scalar conductivity')
pFunction = lambda x,y,z: np.exp(x+y+z)
sigma = scalarConductivity(ccMesh,pFunction)
# Plot sigma model
M.plotImage(sigma)
print sigma
print 'Scalar conductivity defined!'
# -----------------------------------------
SimPEG/tests/test_utils.py
+2 -2
View File
@@ -2,8 +2,7 @@ import numpy as np
import unittest
import sys
sys.path.append('../')
from utils import mkvc, ndgrid, indexCube
from sputils import sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, sp
from utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
class TestSequenceFunctions(unittest.TestCase):
@@ -64,6 +63,7 @@ class TestSequenceFunctions(unittest.TestCase):
self.assertTrue(np.all(indexCube('H', nN) == np.array([10, 11, 13, 14, 19, 20, 22, 23])))
def test_invXXXBlockDiagonal(self):
import scipy.sparse as sp
a = [np.random.rand(5, 1) for i in range(4)]
notebooks/exLomPlots.ipynb
+15 -13
View File
File diff suppressed because one or more lines are too long
notebooks/exPlotImage2D.ipynb
+95 -64
View File
File diff suppressed because one or more lines are too long