forked from simpeg/simpeg
/
ModelBuilder.py
409 lines (296 loc) · 11.2 KB
/
ModelBuilder.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
from __future__ import print_function
import numpy as np
import scipy.ndimage as ndi
import scipy.sparse as sp
from .matutils import mkvc
import sys
if sys.version_info < (3,):
num_types = [int, long, float]
else:
num_types = [int, float]
def addBlock(gridCC, modelCC, p0, p1, blockProp):
"""
Add a block to an exsisting cell centered model, modelCC
:param numpy.array gridCC: mesh.gridCC is the cell centered grid
:param numpy.array modelCC: cell centered model
:param numpy.array p0: bottom, southwest corner of block
:param numpy.array p1: top, northeast corner of block
:blockProp float blockProp: property to assign to the model
:return numpy.array, modelBlock: model with block
"""
ind = getIndicesBlock(p0, p1, gridCC)
modelBlock = modelCC.copy()
modelBlock[ind] = blockProp
return modelBlock
def getIndicesBlock(p0, p1, ccMesh):
"""
Creates a vector containing the block indices in the cell centers 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: 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"
for ii in range(len(p0)):
p0[ii], p1[ii] = np.min([p0[ii], p1[ii]]), np.max([p0[ii], p1[ii]])
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)
elif dimMesh == 3:
# 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 defineBlock(ccMesh, p0, p1, vals=None):
"""
Build a block with the conductivity specified by condVal. Returns an array.
vals[0] conductivity of the block
vals[1] conductivity of the ground
"""
if vals is None:
vals = [0,1]
sigma = np.zeros(ccMesh.shape[0]) + vals[1]
ind = getIndicesBlock(p0,p1,ccMesh)
sigma[ind] = vals[0]
return mkvc(sigma)
def defineElipse(ccMesh, center=None, anisotropy=None, slope=10., theta=0.):
if center is None:
center = [0,0,0]
if anisotropy is None:
anisotropy = [1,1,1]
G = ccMesh.copy()
dim = ccMesh.shape[1]
for i in range(dim):
G[:, i] = G[:,i] - center[i]
theta = -theta*np.pi/180
M = np.array([[np.cos(theta),-np.sin(theta),0],[np.sin(theta),np.cos(theta),0],[0,0,1.]])
M = M[:dim,:dim]
G = M.dot(G.T).T
for i in range(dim):
G[:, i] = G[:,i]/anisotropy[i]*2.
D = np.sqrt(np.sum(G**2,axis=1))
return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5
def getIndicesSphere(center, radius, ccMesh):
"""
Creates a vector containing the sphere indices in the cell centers mesh.
Returns a tuple
The sphere is defined by the points
p0, describe the position of the center of the cell
r, describe the radius of the sphere.
ccMesh represents the cell-centered mesh
The points p0 must live in the the same dimensional space as the mesh.
"""
# Validation: mesh and point (p0) live in the same dimensional space
dimMesh = np.size(ccMesh[0,:])
assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
if dimMesh == 1:
# Define the reference points
ind = np.abs(center[0] - ccMesh[:,0]) < radius
elif dimMesh == 2:
# Define the reference points
ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius
elif dimMesh == 3:
# Define the points
ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius
# Return a tuple
return ind
def defineTwoLayers(ccMesh, depth, vals=None):
"""
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
"""
if vals is None:
vals = [0,1]
sigma = np.zeros(ccMesh.shape[0]) + vals[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]
if dim>1: p0[1] = ccMesh[0,1]
if dim>2: p0[2] = ccMesh[0,2]
# Identify the last cell-centered reference point
p1[0] = ccMesh[-1,0]
if dim>1: p1[1] = ccMesh[-1,1]
if dim>2: p1[2] = ccMesh[-1,2]
# The depth is always defined on the last one.
p1[len(p1)-1] -= depth
ind = getIndicesBlock(p0,p1,ccMesh)
sigma[ind] = vals[0];
return mkvc(sigma)
def scalarConductivity(ccMesh, pFunction):
"""
Define the distribution conductivity in the mesh according to the
analytical expression given in pFunction
"""
dim = np.size(ccMesh[0,:])
CC = [ccMesh[:,0]]
if dim>1: CC.append(ccMesh[:,1])
if dim>2: CC.append(ccMesh[:,2])
sigma = pFunction(*CC)
return mkvc(sigma)
def layeredModel(ccMesh, layerTops, layerValues):
"""
Define a layered model from layerTops (z-positive up)
:param numpy.array ccMesh: cell-centered mesh
:param numpy.array layerTops: z-locations of the tops of each layer
:param numpy.array layerValue: values of the property to assign for each layer (starting at the top)
:rtype: numpy.array
:return: M, layered model on the mesh
"""
descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
# TODO: put an error check to make sure that there is an ordering... needs to work with inf elts
# assert ascending or descending, "Layers must be listed in either ascending or descending order"
# start from bottom up
if not descending:
zprop = np.hstack([mkvc(layerTops,2),mkvc(layerValues,2)])
zprop.sort(axis=0)
layerTops, layerValues = zprop[::-1,0], zprop[::-1,1]
# put in vector form
layerTops, layerValues = mkvc(layerTops), mkvc(layerValues)
# initialize with bottom layer
dim = ccMesh.shape[1]
if dim == 3:
z = ccMesh[:,2]
elif dim == 2:
z = ccMesh[:,1]
elif dim == 1:
z = ccMesh[:,0]
model = np.zeros(ccMesh.shape[0])
for i, top in enumerate(layerTops):
zind = z <= top
model[zind] = layerValues[i]
return model
def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
"""
Create a random model by convolving a kernel with a
uniformly distributed model.
:param tuple shape: shape of the model.
:param int seed: pick which model to produce, prints the seed if you don't choose.
:param numpy.ndarray anisotropy: this is the (3 x n) blurring kernel that is used.
:param int its: number of smoothing iterations
:param list bounds: bounds on the model, len(list) == 2
:rtype: numpy.ndarray
:return: M, the model
.. plot::
import matplotlib.pyplot as plt
import SimPEG.Utils.ModelBuilder as MB
plt.colorbar(plt.imshow(MB.randomModel((50,50),bounds=[-4,0])))
plt.title('A very cool, yet completely random model.')
plt.show()
"""
if bounds is None:
bounds = [0,1]
if seed is None:
seed = np.random.randint(1e3)
print('Using a seed of: ', seed)
if type(shape) in num_types:
shape = (shape,) # make it a tuple for consistency
np.random.seed(seed)
mr = np.random.rand(*shape)
if anisotropy is None:
if len(shape) is 1:
smth = np.array([1,10.,1],dtype=float)
elif len(shape) is 2:
smth = np.array([[1,7,1],[2,10,2],[1,7,1]],dtype=float)
elif len(shape) is 3:
kernal = np.array([1,4,1], dtype=float).reshape((1,3))
smth = np.array(sp.kron(sp.kron(kernal,kernal.T).todense()[:],kernal).todense()).reshape((3,3,3))
else:
assert len(anisotropy.shape) is len(shape), 'Anisotropy must be the same shape.'
smth = np.array(anisotropy,dtype=float)
smth = smth/smth.sum() # normalize
mi = mr
for i in range(its):
mi = ndi.convolve(mi, smth)
# scale the model to live between the bounds.
mi = (mi - mi.min())/(mi.max()-mi.min()) # scaled between 0 and 1
mi = mi*(bounds[1]-bounds[0])+bounds[0]
return mi
if __name__ == '__main__':
from SimPEG.Mesh import TensorMesh
from matplotlib import pyplot as plt
# Define the mesh
testDim = 2
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])[:testDim]
p1 = np.array([1.0,1.0,1.0])[:testDim]
vals = np.array([100,1e-6])
sigma = defineBlockConductivity(ccMesh,p0,p1,vals)
# Plot sigma model
print(sigma.shape)
M.plotImage(sigma)
print('Done with block! :)')
plt.show()
# -----------------------------------------
print('Testing the two layered model')
vals = np.array([100,1e-5]);
depth = 1.0;
sigma = defineTwoLayeredConductivity(ccMesh,depth,vals)
M.plotImage(sigma)
print(sigma)
print('layer model!')
plt.show()
# -----------------------------------------
print('Testing scalar conductivity')
if testDim == 1:
pFunction = lambda x: np.exp(x)
elif testDim == 2:
pFunction = lambda x,y: np.exp(x+y)
elif testDim == 3:
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!')
plt.show()
# -----------------------------------------