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

Merged in condModels (pull request #10)

Browse Source 
Cond Models
This commit is contained in:
rowanc1
2013-08-30 22:25:11 -07:00
parent d651a2bcbc c0efb972a9
commit 82e8df226c
11 changed files with 338 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/BaseMesh.py
+24 -17
View File
@@ -59,22 +59,22 @@ class BaseMesh(object):
For example, you have a face variable, and you want the x component of it reshaped to a 3D matrix.
Mesh.r can fulfil your dreams...
Mesh.r can fulfil your dreams::
mesh.r(V, 'F', 'Fx', 'M')
| | | { How: 'M' or ['V'] for a matrix (ndgrid style) or a vector (n x dim) }
| | { What you want: ['CC'], 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez' }
| { What is it: ['CC'], 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez' }
{ The input: as a list or ndarray }
mesh.r(V, 'F', 'Fx', 'M')
| | | { How: 'M' or ['V'] for a matrix (ndgrid style) or a vector (n x dim) }
| | { What you want: ['CC'], 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez' }
| { What is it: ['CC'], 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez' }
{ The input: as a list or ndarray }
For example:
For example::
Xex, Yex, Zex = r(mesh.gridEx, 'Ex', 'Ex', 'M') # Separates each component of the Ex grid into 3 matrices
Xex, Yex, Zex = r(mesh.gridEx, 'Ex', 'Ex', 'M') # Separates each component of the Ex grid into 3 matrices
XedgeVector = r(edgeVector, 'E', 'Ex', 'V') # Given an edge vector, this will return just the part on the x edges as a vector
XedgeVector = r(edgeVector, 'E', 'Ex', 'V') # Given an edge vector, this will return just the part on the x edges as a vector
eX, eY, eZ = r(edgeVector, 'E', 'E', 'V') # Separates each component of the edgeVector into 3 vectors
eX, eY, eZ = r(edgeVector, 'E', 'E', 'V') # Separates each component of the edgeVector into 3 vectors
"""
assert (type(x) == list or type(x) == np.ndarray), "x must be either a list or a ndarray"
@@ -341,18 +341,25 @@ class BaseMesh(object):
tangents = property(**tangents())
def projectFaceVector(self, fV):
"""Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals"""
"""
Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
:param numpy.array fV: face vector with shape (nF, dim)
:rtype: numpy.array with shape (nF, )
:return: projected face vector
"""
assert type(fV) == np.ndarray, 'fV must be an ndarray'
assert len(fV.shape) == 2 and fV.shape[0] == np.sum(self.nF) and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
return np.sum(fV*self.normals, 1)
def projectEdgeVector(self, eV):
"""Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents"""
"""
Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
:param numpy.array eV: edge vector with shape (nE, dim)
:rtype: numpy.array with shape (nE, )
:return: projected edge vector
"""
assert type(eV) == np.ndarray, 'eV must be an ndarray'
assert len(eV.shape) == 2 and eV.shape[0] == np.sum(self.nE) and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
return np.sum(eV*self.tangents, 1)
if __name__ == '__main__':
m = BaseMesh([3, 2, 4])
print m.n
SimPEG/TensorMesh.py
+5 -2
View File
@@ -13,7 +13,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
Any Mesh that has a constant width along the entire axis
such that it can defined by a single width vector, called 'h'.
e.g.
::
hx = np.array([1,1,1])
hy = np.array([1,2])
@@ -21,6 +21,9 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
mesh = TensorMesh([hx, hy, hz])
.. math::
x^2 = 5
"""
_meshType = 'TENSOR'
@@ -61,7 +64,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h)
return outStr[:-1]
if self.dim == 1:
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
SimPEG/TensorView.py
+15 -13
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)
@@ -124,12 +124,13 @@ class TensorView(object):
# 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
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
@@ -141,17 +142,18 @@ 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.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), np.sum(self.hy)).ravel())]
gy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), np.sum(self.hx)).ravel())]
gx = np.arange(nX+1)*self.vectorNx[-1]
gy = np.arange(nY+1)*self.vectorNy[-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, '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
@@ -160,9 +162,9 @@ 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')
fig.show()
plt.show()
return ph
def plotGrid(self):
@@ -180,7 +182,7 @@ class TensorView(object):
ax.grid(True)
ax.hold(False)
ax.set_xlabel('x1')
fig.show()
plt.show()
elif self.dim == 2:
fig = plt.figure(2)
fig.clf()
@@ -199,7 +201,7 @@ class TensorView(object):
ax.hold(False)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
fig.show()
plt.show()
elif self.dim == 3:
fig = plt.figure(3)
fig.clf()
@@ -228,4 +230,4 @@ class TensorView(object):
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('x3')
fig.show()
plt.show()
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)]
SimPEG/utils/ModelBuilder.py
+194
View File
@@ -0,0 +1,194 @@
import numpy as np
def getIndecesBlock(p0,p1,ccMesh):
"""
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
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.
"""
# 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]
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]
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]
z1 = p0[2]
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
return ind
def defineBlockConductivity(p0,p1,ccMesh,condVals):
"""
Build a block with the conductivity specified by condVal. Returns an array.
condVals[0] conductivity of the block
condVals[1] conductivity of the ground
"""
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:
<----------------------------|------------------------------------>
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]
p0[2] = ccMesh[0,2]
# Identify the last cell-centered reference point
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
"""
xCC = ccMesh[:,0]
yCC = ccMesh[:,1]
zCC = ccMesh[:,2]
sigma = pFunction(xCC,yCC,zCC)
return sigma
if __name__ == '__main__':
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
# Define the mesh
testDim = 3
h1 = 0.3*np.ones(7)
h1[0] = 0.5
h1[-1] = 0.6
h2 = .5 * np.ones(4)
h3 = .4 * np.ones(6)
x0 = np.zeros(3)
if testDim == 1:
h = [h1]
x0 = x0[0]
elif testDim == 2:
h = [h1, h2]
x0 = x0[0:2]
else:
h = [h1, h2, h3]
M = 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)
# Plot sigma model
print sigma.shape
M.plotImage(sigma)
print 'Done with block! :)'
# -----------------------------------------
print('Testing the two layered model')
condVals = np.array([100,1e-5]);
depth = 1.0;
sigma = defineTwoLayeredConductivity(depth,ccMesh,condVals)
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/utils/__init__.py
+2 -1
View File
@@ -1,3 +1,4 @@
from utils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
import ModelBuilder
SimPEG/utils/lomutils.py
+1 -1
View File
@@ -1,6 +1,6 @@
import numpy as np
from scipy import sparse as sp
from utils import mkvc, ndgrid, sub2ind
from matutils import mkvc, ndgrid, sub2ind
from sputils import sdiag
SimPEG/utils/utils.py → SimPEG/utils/matutils.py
View File
SimPEG/utils/sputils.py
+1 -1
View File
@@ -1,5 +1,5 @@
from scipy import sparse as sp
from utils import mkvc
from matutils import mkvc
def sdiag(h):
notebooks/exLomPlots.ipynb
+15 -13
View File
File diff suppressed because one or more lines are too long
notebooks/exPlotImage2D.ipynb
+79 -57
View File
File diff suppressed because one or more lines are too long