-
Notifications
You must be signed in to change notification settings - Fork 35
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Ref/properties #70
Ref/properties #70
Changes from all commits
fa11a0a
145bf00
19882a2
af1035d
693c5d0
809a149
44fea79
d688b85
209856a
85a4ba0
5fc99d5
d9d59cc
2606404
24accb0
d347514
7b03487
f8c9c62
393851f
e770ff9
8b061fe
f471d17
fd85204
abee619
5507fc1
2930503
d76cd60
018d585
5060219
5980f43
a72837e
23b5f45
76d78c5
d15dda3
ad2d77a
a75afa5
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,4 +1,4 @@ | ||
[bumpversion] | ||
current_version = 0.1.9 | ||
current_version = 0.1.10b0 | ||
files = setup.py discretize/__init__.py docs/conf.py | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,41 +1,102 @@ | ||
import numpy as np | ||
from discretize import utils | ||
import properties | ||
import os | ||
import json | ||
from .utils.matutils import mkvc | ||
|
||
|
||
class BaseMesh(object): | ||
"""BaseMesh does all the counting you don't want to do. | ||
class BaseMesh(properties.HasProperties): | ||
""" | ||
BaseMesh does all the counting you don't want to do. | ||
BaseMesh should be inherited by meshes with a regular structure. | ||
|
||
:param numpy.array n: (or list) number of cells in each direction (dim, ) | ||
:param numpy.array x0: (or list) Origin of the mesh (dim, ) | ||
|
||
""" | ||
|
||
def __init__(self, n, x0=None): | ||
_REGISTRY = {} | ||
|
||
# Check inputs | ||
if x0 is None: | ||
x0 = np.zeros(len(n)) | ||
# Properties | ||
_n = properties.Array( | ||
"number of cells in each direction (dim, )", | ||
dtype=int, | ||
required=True, | ||
shape=('*',) | ||
) | ||
|
||
if not len(n) == len(x0): | ||
raise Exception("Dimension mismatch. x0 != len(n)") | ||
x0 = properties.Array( | ||
"origin of the mesh (dim, )", | ||
dtype=float, | ||
shape=('*',), | ||
required=True | ||
) | ||
|
||
if len(n) > 3: | ||
raise Exception("Dimensions higher than 3 are not supported.") | ||
# Instantiate the class | ||
def __init__(self, n, x0=None): | ||
self._n = n # number of dimensions | ||
|
||
# Ensure x0 & n are 1D vectors | ||
self._n = np.array(n, dtype=int).ravel() | ||
self._x0 = np.array(x0, dtype=float).ravel() | ||
self._dim = len(self._x0) | ||
if x0 is None: | ||
self.x0 = np.zeros(len(self._n)) | ||
else: | ||
self.x0 = x0 | ||
|
||
@property | ||
def x0(self): | ||
"""Origin of the mesh | ||
# Validators | ||
@properties.validator('_n') | ||
def check_n_shape(self, change): | ||
assert ( | ||
not isinstance(change['value'], properties.utils.Sentinel) and | ||
change['value'] is not None | ||
), "Cannot delete n. Instead, create a new mesh" | ||
|
||
change['value'] = np.array(change['value'], dtype=int).ravel() | ||
if len(change['value']) > 3: | ||
raise Exception( | ||
"Dimensions of {}, which is higher than 3 are not " | ||
"supported".format(change['value']) | ||
) | ||
|
||
if change['previous'] != properties.undefined: | ||
# can't change dimension of the mesh | ||
assert len(change['previous']) == len(change['value']), ( | ||
"Cannot change dimensionality of the mesh. Expected {} " | ||
"dimensions, got {} dimensions".format( | ||
len(change['previous']), len(change['value']) | ||
) | ||
) | ||
|
||
# check that if h has been set, sizes still agree | ||
if getattr(self, 'h', None) is not None and len(self.h) > 0: | ||
for i in range(len(change['value'])): | ||
assert len(self.h[i]) == change['value'][i], ( | ||
"Mismatched shape of n. Expected {}, len(h[{}]), got " | ||
"{}".format( | ||
len(self.h[i]), i, change['value'][i] | ||
) | ||
) | ||
|
||
# check that if nodes have been set for curvi mesh, sizes still | ||
# agree | ||
if ( | ||
getattr(self, 'nodes', None) is not None and | ||
len(self.nodes) > 0 | ||
): | ||
for i in range(len(change['value'])): | ||
assert self.nodes[0].shape[i]-1 == change['value'][i], ( | ||
"Mismatched shape of n. Expected {}, len(nodes[{}]), " | ||
"got {}".format( | ||
self.nodes[0].shape[i]-1, i, change['value'][i] | ||
) | ||
) | ||
|
||
@properties.validator('x0') | ||
def check_x0(self, change): | ||
assert ( | ||
not isinstance(change['value'], properties.utils.Sentinel) and | ||
change['value'] is not None | ||
), "n must be set prior to setting x0" | ||
|
||
:rtype: numpy.array | ||
:return: x0, (dim, ) | ||
""" | ||
return self._x0 | ||
if len(self._n) != len(change['value']): | ||
raise Exception( | ||
"Dimension mismatch. x0 has length {} != len(n) which is " | ||
"{}".format(len(x0), len(n)) | ||
) | ||
|
||
@property | ||
def dim(self): | ||
|
@@ -44,7 +105,7 @@ def dim(self): | |
:rtype: int | ||
:return: dim | ||
""" | ||
return self._dim | ||
return len(self._n) | ||
|
||
@property | ||
def nC(self): | ||
|
@@ -290,11 +351,33 @@ def projectEdgeVector(self, eV): | |
), 'eV must be an ndarray of shape (nE x dim)' | ||
return np.sum(eV*self.tangents, 1) | ||
|
||
def save(self, filename='mesh.json', verbose=False): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This might want to be on a base class for simpeg in future (or in the utils of discretize)? It looks like a pretty generic function. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Agreed. Once we have injected properties elsewhere in SimPEG we can abstract this up. |
||
""" | ||
Save the mesh to json | ||
:param str file: filename for saving the casing properties | ||
:param str directory: working directory for saving the file | ||
""" | ||
|
||
f = os.path.abspath(filename) # make sure we are working with abs path | ||
with open(f, 'w') as outfile: | ||
json.dump(self.serialize(), outfile) | ||
|
||
if verbose is True: | ||
print('Saved {}'.format(f)) | ||
|
||
return f | ||
|
||
def copy(self): | ||
""" | ||
Make a copy of the current mesh | ||
""" | ||
return properties.copy(self) | ||
|
||
|
||
class BaseRectangularMesh(BaseMesh): | ||
"""BaseRectangularMesh""" | ||
def __init__(self, n, x0=None): | ||
BaseMesh.__init__(self, n, x0) | ||
BaseMesh.__init__(self, n, x0=x0) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. this isn't strictly necessary. just for clarity |
||
|
||
@property | ||
def nCx(self): | ||
|
@@ -590,44 +673,69 @@ def r(self, x, xType='CC', outType='CC', format='V'): | |
eX, eY, eZ = r(edgeVector, 'E', 'E', 'V') | ||
""" | ||
|
||
assert (type(x) == list or isinstance(x, np.ndarray)), "x must be either a list or a ndarray" | ||
assert xType in ['CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', 'Ez'], "xType must be either 'CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'" | ||
assert outType in ['CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', 'Ez'], "outType must be either 'CC', 'N', 'F', Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'" | ||
allowed_xType = [ | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. whitespace cleanup |
||
'CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', 'Ez' | ||
] | ||
assert ( | ||
type(x) == list or isinstance(x, np.ndarray) | ||
), "x must be either a list or a ndarray" | ||
assert xType in allowed_xType, ( | ||
"xType must be either " | ||
"'CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'" | ||
) | ||
assert outType in allowed_xType, ( | ||
"outType must be either " | ||
"'CC', 'N', 'F', Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'" | ||
) | ||
assert format in ['M', 'V'], "format must be either 'M' or 'V'" | ||
assert outType[:len(xType)] == xType, "You cannot change types when reshaping." | ||
assert xType in outType, 'You cannot change type of components.' | ||
assert outType[:len(xType)] == xType, ( | ||
"You cannot change types when reshaping." | ||
) | ||
assert xType in outType, "You cannot change type of components." | ||
|
||
if type(x) == list: | ||
for i, xi in enumerate(x): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It looks like this never gets hit by tests. Did it before? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The only thing I changed above is whitespace. This particular line has not changed at all in this pr |
||
assert isinstance(x, np.ndarray), "x[{0:d}] must be a numpy array".format(i) | ||
assert xi.size == x[0].size, "Number of elements in list must not change." | ||
assert isinstance(x, np.ndarray), ( | ||
"x[{0:d}] must be a numpy array".format(i) | ||
) | ||
assert xi.size == x[0].size, ( | ||
"Number of elements in list must not change." | ||
) | ||
|
||
x_array = np.ones((x.size, len(x))) | ||
# Unwrap it and put it in a np array | ||
for i, xi in enumerate(x): | ||
x_array[:, i] = utils.mkvc(xi) | ||
x_array[:, i] = mkvc(xi) | ||
x = x_array | ||
|
||
assert isinstance(x, np.ndarray), "x must be a numpy array" | ||
|
||
x = x[:] # make a copy. | ||
xTypeIsFExyz = len(xType) > 1 and xType[0] in ['F', 'E'] and xType[1] in ['x', 'y', 'z'] | ||
xTypeIsFExyz = ( | ||
len(xType) > 1 and | ||
xType[0] in ['F', 'E'] and | ||
xType[1] in ['x', 'y', 'z'] | ||
) | ||
|
||
def outKernal(xx, nn): | ||
"""Returns xx as either a matrix (shape == nn) or a vector.""" | ||
if format == 'M': | ||
return xx.reshape(nn, order='F') | ||
elif format == 'V': | ||
return utils.mkvc(xx) | ||
return mkvc(xx) | ||
|
||
def switchKernal(xx): | ||
"""Switches over the different options.""" | ||
if xType in ['CC', 'N']: | ||
nn = (self._n) if xType == 'CC' else (self._n+1) | ||
assert xx.size == np.prod(nn), "Number of elements must not change." | ||
assert xx.size == np.prod(nn), ( | ||
"Number of elements must not change." | ||
) | ||
return outKernal(xx, nn) | ||
elif xType in ['F', 'E']: | ||
# This will only deal with components of fields, not full 'F' or 'E' | ||
xx = utils.mkvc(xx) # unwrap it in case it is a matrix | ||
# This will only deal with components of fields, | ||
# not full 'F' or 'E' | ||
xx = mkvc(xx) # unwrap it in case it is a matrix | ||
nn = self.vnF if xType == 'F' else self.vnE | ||
nn = np.r_[0, nn] | ||
|
||
|
@@ -638,13 +746,20 @@ def switchKernal(xx): | |
|
||
for dim, dimName in enumerate(['x', 'y', 'z']): | ||
if dimName in outType: | ||
assert self.dim > dim, ("Dimensions of mesh not great enough for %s%s", (xType, dimName)) | ||
assert xx.size == np.sum(nn), 'Vector is not the right size.' | ||
assert self.dim > dim, ( | ||
"Dimensions of mesh not great enough for " | ||
"{}{}".format(xType, dimName) | ||
) | ||
assert xx.size == np.sum(nn), ( | ||
"Vector is not the right size." | ||
) | ||
start = np.sum(nn[:dim+1]) | ||
end = np.sum(nn[:dim+2]) | ||
return outKernal(xx[start:end], nx[dim]) | ||
|
||
elif xTypeIsFExyz: | ||
# This will deal with partial components (x, y or z) lying on edges or faces | ||
# This will deal with partial components (x, y or z) | ||
# lying on edges or faces | ||
if 'x' in xType: | ||
nn = self.vnFx if 'F' in xType else self.vnEx | ||
elif 'y' in xType: | ||
|
@@ -658,7 +773,9 @@ def switchKernal(xx): | |
isVectorQuantity = len(x.shape) == 2 and x.shape[1] == self.dim | ||
|
||
if outType in ['F', 'E']: | ||
assert ~isVectorQuantity, 'Not sure what to do with a vector vector quantity..' | ||
assert ~isVectorQuantity, ( | ||
'Not sure what to do with a vector vector quantity..' | ||
) | ||
outTypeCopy = outType | ||
out = () | ||
for ii, dirName in enumerate(['x', 'y', 'z'][:self.dim]): | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
we weren't really using this and they discontinued it